32 lines
763 B
Fortran
32 lines
763 B
Fortran
FUNCTION YLINTP (XINT,X,Y,N,NTOT)
|
|
C =================================
|
|
C
|
|
C linear interpolation routine. Determines YINT = Y(XINT) from
|
|
C grid Y(X) with N points and dimension NTOT.
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
DIMENSION X(NTOT),Y(NTOT)
|
|
C
|
|
C bisection (see Numerical Recipes par 3.4 page 90)
|
|
JL = 0
|
|
JU = N+1
|
|
10 IF (JU-JL.GT.1) THEN
|
|
JM = (JU+JL)/2
|
|
IF ((X(N).GT.X(1)).EQV.(XINT.GT.X(JM))) THEN
|
|
JL = JM
|
|
ELSE
|
|
JU = JM
|
|
END IF
|
|
GO TO 10
|
|
END IF
|
|
J = JL
|
|
IF (J.EQ.N) J = J-1
|
|
IF (J.EQ.0) J = J+1
|
|
CM42 version: enable extrapolation
|
|
C
|
|
RC = (Y(J+1)-Y(J))/(X(J+1)-X(J))
|
|
YLINTP = RC*(XINT-X(J))+Y(J)
|
|
C
|
|
RETURN
|
|
END
|