.
下面是一个样条拟合程序:
C=======================================================================
C
C THIRD SPLINE INTERPLATION
C
C=======================================================================
SUBROUTINE SPL(N,K,T,Y,Y1,Y2,SUM,U1,V1)
DIMENSION U1(N),V1(N)
REAL M(100),H(100)
DIMENSION A(100),B(100),C(100),D(100),U(100),V(100)
IF(N.EQ.1) THEN
Y=V1(1)
GOTO 90
ENDIF
IF(N.EQ.2) THEN
Y=(V1(2)-V1(1))/(U1(2)-U1(1))*(T-U1(1))+V1(1)
GOTO 90
ENDIF
IF(N.EQ.3) THEN
IF(T.GE.U1(2)) THEN
Y=(V1(3)-V1(2))/(U1(3)-U1(2))*(T-U1(2))+V1(2)
ELSE
Y=(V1(2)-V1(1))/(U1(2)-U1(1))*(T-U1(1))+V1(1)
ENDIF
GOTO 90
ENDIF
IF(U1(1).GT.U1(N)) THEN
DO 12 I=1,N
U(I)=U1(N-I+1)
12 V(I)=V1(N-I+1)
ELSE
DO 8 I=1,N
U(I)=U1(I)
8 V(I)=V1(I)
ENDIF
MA=0
DO 73 I=1,N-1
IF(U(I).GE.U(I+1)) THEN
MA=MA+1
N=N-1
ENDIF
U(I)=U(I+MA)
73 V(I)=V(I+MA)
DO 5 I=2,N
H(I)=U(I)-U(I-1)
5 B(I)=2.
B(1)=2.
C(1)=0.
D(1)=0.
A(N)=0.
D(N)=0.
DO 10 I=2,N-1
A(I)=H(I)/(H(I)+H(I+1))
10 C(I)=1.0-A(I)
DO 15 I=2,N-1
15 D(I)=6.*((V(I+1)-V(I))/H(I+1)-(V(I)-V(I-1))/
* H(I))/(H(I)+H(I+1))
D(1)=D(1)/B(1)
W=B(1)
DO 20 I=2,N
I1=I-1
B(I1)=C(I-1)/W
W=B(I)-A(I)*B(I-1)
20 D(I)=(D(I)-A(I)*D(I-1))/W
DO 25 I=1,N-1
NI=N-I
D(NI)=D(NI)-B(NI)*D(NI+1)
25 M(NI)=D(NI)
IF((T-U(1)).LT.0..OR.(T-U(N)).GT.0.) GOTO 35
GOTO 40
35 CONTINUE
IF(T.LT.U(1)) THEN
AS=(V(3)-V(1))/(U(3)-U(1))
BS=(V(2)-V(1))/(U(2)-U(1))
AS=(AS-BS)/(U(3)-U(2))
Y=V(1)+(T-U(1))*BS+(T-U(1))*(T-U(2))*AS
ELSE
AS=(V(N)-V(N-2))/(U(N)-U(N-2))
BS=(V(N-1)-V(N-2))/(U(N-1)-U(N-2))
AS=(AS-BS)/(U(N)-U(N-1))
Y=V(N-2)+(T-U(N-2))*BS+(T-U(N-2))*(T-U(N-1))*AS
ENDIF
RETURN
40 I=2
47 IF(T.GT.U(I)) GOTO 50
GOTO 55
50 I=I+1
GOTO 47
55 W=(T-U(I-1))/H(I)
M(N)=0.
GOTO (75,70,60,65,65) ,K
60 SUM=0.
DO 62 I=2,N
62 SUM=SUM+(V(I-1)+V(I))*H(I)/2.-(M(I-1)+M(I))*H(I)**3/24.
IF(K.EQ.3.) GOTO 80
65 Y2=M(I-1)*(1.-W)+M(I)*W
IF(K.EQ.5) GOTO 80
70 Y1=(V(I)-V(I-1))/H(I)+H(I)*(M(I-1)
* *(-W*W/2.+W-1./3.)+M(I)*(W*W-1./3.)/2.)
IF(K.EQ.2) GOTO 80
75 Y=V(I-1)*(1.-W)+V(I)*W+W*(W-1.)*H(I)*H(I)*
* (M(I-1)*(2.-W)+M(I)*(W+1.))/6.
80 CONTINUE
90 RETURN
END
C======================================================================= |