c
C插值! This subroutine is the Lagrange the interpolation routine progrmae
C
SUBROUTINE INTERPOLATION(XA,YA,N,X,Y,DY)
C INPUT:
C N: 存放插值节点数
C XA(N):存放插值节点值
C YA(N):存放函数值
C X:输入参数,插值自变量
C OUTPUT:
C Y:输出参数,所求的插值点的函数值
C DY:输出参数,误差估计
REAL XA(N),YA(N),X,Y,C(11),D(11)
DOUBLE PRECISION DY,DIF,DIFT,HO,HP,W,DEN
NS=1
DIF=ABS(X-XA(1))
DO I=1,N
DIFT=ABS(X-XA(I))
IF(DIFT.LT.DIF) THEN
NS=I
DIF=DIFT
END IF
C(I)=YA(I)
D(I)=YA(I)
END DO
Y=YA(NS)
NS=NS-1
DO M=1,N-1
DO I=1,N-M
HO=XA(I)-X
HP=XA(I+M)-X
W=C(I+1)-D(I)
DEN=HO-HP
IF(DEN.EQ.0.0) RETURN
DEN=W/DEN
D(I)=HP*DEN
C(I)=HO*DEN
END DO
IF((2*NS).LT.(N-M)) THEN
DY=C(NS+1)
ELSE
DY=D(NS)
NS=NS-1
END IF
Y=Y+DY
END DO
END
主程序!
PROGRAM MAIN
INTEGER I,J,K,M,N
REAL XA(8),YA(8),YYA(11),XA1(11),YA1(11)
REAL X,Y
DOUBLE PRECISION DY
REAL LANDRF(8),TOA(11),REF(6,7,11,8)
LANDRF(1)=0
LANDRF(2)=0.05
LANDRF(3)=0.15
LANDRF(4)=0.20
LANDRF(5)=0.25
LANDRF(6)=0.30
LANDRF(7)=0.35
LANDRF(8)=0.40
DO I=1,11
TOA(I)=(I-1)*0.1
END DO
CALL READ_LOOK(REF)
DO J=1,11
X=0.15
DO I=1,8
XA(I)=LANDRF(I)
YA(I)=REF(6,4,J,I)
END DO
N=8
CALL INTERPOLATION(XA,YA,N,X,Y,DY)
YYA(J)=Y
END DO
open(2,file='test.txt')
DO I=1,11
WRITE(2,*) YYA(I)
END DO
close(2)
X=0.1644
DO I=1,11
XA1(I)=YYA(I)
YA1(I)=(I-1)*0.1
END DO
N=11
CALL INTERPOLATION(XA1,YA1,N,X,Y,DY)
WRITE(*,*) Y,DY
END
读程序!
SUBROUTINE READ_LOOK(REF)
C PROGRAM MAIN
REAL REF(6,7,11,8)
INTEGER MOLD,IWAV,NTAU,LR
character*50 FILE
C LOOP IS AROUND THET0,NEW FILE FOR EACH THETA0
open(1, FILE='MA.TXT')
DO MOLD=1,6
DO IWAV=1,7
DO NTAU=1,11
DO LR=1,8
READ(1,*) REF(MOLD,IWAV,NTAU,LR)
END DO
END DO
END DO
END DO
CLOSE(1)
C 90 FORMAT(f10.4)
C DO LR=1,8
C WRITE(*,*) REF(1,2,1,LR)
C END DO
RETURN
END