三次样条插值

!==================================================================================
!                                                                                 =     
!            本程序中均为三次样条插值计算                                         =
!----------------------------------------------------------------------------------
!==================================================================================        
! 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