!==================================================================================
! =
! 本程序中均为三次样条插值计算 =
!----------------------------------------------------------------------------------
!==================================================================================
! X-----------自变量的X坐标 *
! Y-----------因变量 *
! N-----------变量个数 *
! S2----------(插值系数)各节点对应的y的二阶导数值,即三弯矩方程的解 *
! DY----------(插值系数)(yi+1-yi)/(xi+1-xi) *
! S-----------(插值系数)Y的二阶导数在区间(xi,xi+1)上的斜率 *
! SX----------需要计算插值的各点的X坐标 *
! M-----------x==sx需要计算插值的数组容量大小 *
! F-----------x==sx处插值结果 *
! F1----------x==sx处插值结果的一阶导数值 *
! F2----------x==sx处插值结果的二阶导数值 *
!==================================================================================
SUBROUTINE SPLINE(X,Y,N,SX,M,F,F1,F2) !三次样条插值程序
IMPLICIT NONE
INTEGER:: i, j, k, M, N, N1, N2
Real(KIND=8):: X(N),Y(N),SX(M),F(M),F1(M),F2(M)
! Real(KIND=8):: S2(N), H(N1), DY(N1), S(N1), E(N2) !直接都给数组定义N的空间大小,方便调用
Real(KIND=8):: S2(N), H(N), DY(N), S(N), E(N)
Real(KIND=8):: Z, H1, H2, H3, H4
N1=N-1
N2=N-2
DO i = 1, N1
H(i) = X(i+1) - X(i)
IF(H(i)==0.0) THEN
WRITE(1000,*) '错误:SPLINE中H(I)==0.0,行数为:',I
WRITE(*,*) '错误:SPLINE中H(I)==0.0,行数为:',I
STOP
ENDIF
DY(i) = ( Y(i+1) - Y(i) ) / H(i)
ENDDO
S2(1) = 0.d0; S2(N) = 0.d0
DO i = 2, N1
S2(i) = 6.d0 * ( DY(i) - DY(i-1) )
ENDDO
Z = 0.5d0 / ( H(1) + H(2) )
S(1) = -H(2) * Z
E(1) = S2(2) * Z
DO i = 2, N2
k = i - 1
j = i + 1
Z = 1.d0 / ( 2.d0*( H(i)+H(j) ) + H(i)*S(k) )
S(i) = -H(j) * Z
E(i) = ( S2(j)-H(i)*E(k) ) * Z
ENDDO
S2(N1) = E(N2)
DO i = N2, 2, -1
k = i - 1
S2(i) = S(k)*S2(i+1) + E(k)
ENDDO
DO i = 1, N1
S(i) = ( S2(i+1) - S2(i) ) / H(i)
ENDDO
!具体求结果
i = 2
k = 1
DO j = 1, M
DO
IF(SX(j)<X(1).OR.SX(j)>X(N)) THEN
F(j)=0.0
F1(j)=0.0
F2(j)=0.0
ELSEIF ( SX(j) > X(i) ) THEN
k = i
i = i + 1
ELSE
EXIT
ENDIF
ENDDO
H1 = SX(j) - X(k)
H2 = SX(j) - X(i)
H3 = H1 * H2
H4 = S2(k) + H1*S(k)
Z = ( S2(i) + S2(k) + H4 ) / 6.d0
F(j) = Y(k) + H1*DY(k) + H3*Z
F1(j) = DY(k) + Z*( H1+H2 ) + H3 * S(k) / 6.d0
F2(j)=H4
ENDDO
END SUBROUTINE SPLINE