95 lines
2.2 KiB
Fortran
95 lines
2.2 KiB
Fortran
SUBROUTINE INTERP(X,Y,XX,YY,NX,NXX,NPOL,ILOGX,ILOGY)
|
|
C ====================================================
|
|
C
|
|
C General interpolation procedure of the (NPOL-1)-th order
|
|
C
|
|
C for ILOGX = 1 logarithmic interpolation in X
|
|
C for ILOGY = 1 logarithmic interpolation in Y
|
|
C
|
|
C Input:
|
|
C X - array of original x-coordinates
|
|
C Y - array of corresponding functional values Y=y(X)
|
|
C NX - number of elements in arrays X or Y
|
|
C XX - array of new x-coordinates (to which is to be
|
|
C interpolated
|
|
C NXX - number of elements in array XX
|
|
C Output:
|
|
C YY - interpolated functional values YY=y(XX)
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
DIMENSION X(MDEPTH),Y(MDEPTH),XX(MDEPTH),YY(MDEPTH)
|
|
C
|
|
C no interpolation for NPOL.LE.0 or NX.le.0
|
|
C
|
|
IF(NPOL.LE.0.OR.NX.LE.0) THEN
|
|
N=NX
|
|
IF(NXX.GE.NX) N=NXX
|
|
DO I=1,N
|
|
XX(I)=X(I)
|
|
YY(I)=Y(I)
|
|
END DO
|
|
RETURN
|
|
END IF
|
|
C
|
|
C interpolation
|
|
C
|
|
C if required, compute logarithms to be interpolated
|
|
C
|
|
IF(ILOGX.GT.0) THEN
|
|
DO I=1,NX
|
|
X(I)=LOG10(X(I))
|
|
END DO
|
|
DO I=1,NXX
|
|
XX(I)=LOG10(XX(I))
|
|
END DO
|
|
END IF
|
|
IF(ILOGY.GT.0) THEN
|
|
DO I=1,NX
|
|
Y(I)=LOG10(Y(I))
|
|
END DO
|
|
END IF
|
|
C
|
|
NM=(NPOL+1)/2
|
|
NM1=NM+1
|
|
NUP=NX+NM1-NPOL
|
|
DO ID=1,NXX
|
|
XXX=XX(ID)
|
|
DO I=NM1,NUP
|
|
IF(XXX.LE.X(I)) GO TO 10
|
|
END DO
|
|
I=NUP
|
|
10 J=I-NM
|
|
JJ=J+NPOL-1
|
|
YYY=0.
|
|
DO K=J,JJ
|
|
T=1.
|
|
DO 20 M=J,JJ
|
|
IF(K.EQ.M) GO TO 20
|
|
T=T*(XXX-X(M))/(X(K)-X(M))
|
|
20 CONTINUE
|
|
YYY=Y(K)*T+YYY
|
|
END DO
|
|
YY(ID)=YYY
|
|
END DO
|
|
C
|
|
IF(ILOGX.GT.0) THEN
|
|
DO I=1,NX
|
|
X(I)=EXP(X(I)*2.30258509299405)
|
|
END DO
|
|
DO I=1,NXX
|
|
XX(I)=EXP(XX(I)*2.30258509299405)
|
|
END DO
|
|
END IF
|
|
IF(ILOGY.GT.0) THEN
|
|
DO I=1,NX
|
|
Y(I)=EXP(Y(I)*2.30258509299405)
|
|
END DO
|
|
DO I=1,NXX
|
|
YY(I)=EXP(YY(I)*2.30258509299405)
|
|
END DO
|
|
END IF
|
|
C
|
|
RETURN
|
|
END
|