C
C PAGE 277-278: NUMERICAL MATHEMATICS AND COMPUTING, CHENEY/KINCAID, 1985
C
C FILE: SPL3.FOR
C
C NATURAL CUBIC SPLINE FUNCTION FOR SIN(X) AT EQUIDISTANT POINTS (SPL3,ZSPL3)
C
      PARAMETER  (N = 10)   
      DIMENSION T(N),Y(N),U(N),V(N),Z(N),W(N),B(N)
      H = 1.6875/REAL(N-1)
      DO 2 I=1,N  
        T(I) = REAL(I-1)*H  
        Y(I) = SIN(T(I))    
   2  CONTINUE    
      CALL ZSPL3(N,T,Y,W,B,U,V,Z)     
      DO 3 I = 1,4*N-3      
        X = REAL(I-1)*H*0.25
        D = SIN(X) - SPL3(N,T,Y,Z,X)  
        PRINT 4,I,X,D       
   3  CONTINUE    
   4  FORMAT(2X,I5,F22.14,E10.3)      
      STOP
      END 
  
      SUBROUTINE ZSPL3(N,T,Y,H,B,U,V,Z) 
      DIMENSION  T(N),Y(N),H(N),B(N),U(N),V(N),Z(N)       
      DO 2 I = 1,N-1
        H(I) = T(I+1) - T(I)
        B(I) = (Y(I+1) -Y(I))/H(I)    
   2  CONTINUE    
      U(2) = 2.0*(H(1) + H(2))
      V(2) = 6.0*(B(2) - B(1))
      DO 3 I = 3,N-1
        U(I) = 2.0*(H(I) + H(I-1)) - H(I-1)**2/U(I-1)     
        V(I) = 6.0*(B(I) - B(I-1)) - H(I-1)*V(I-1)/U(I-1) 
   3  CONTINUE    
      Z(N) = 0.0  
      DO 4 I = N-1,2,-1     
        Z(I) = (V(I) - H(I)*Z(I+1))/U(I)
   4  CONTINUE    
      Z(1) = 0.0
      RETURN
      END 
  
      FUNCTION SPL3(N,T,Y,Z,X)
      DIMENSION  T(N),Y(N),Z(N)       
      DO 2 I = N-1,1,-1     
        DIFF = X - T(I)     
        IF(DIFF .GE. 0.0)  GO TO 3    
   2  CONTINUE    
      I = 1 
   3  H = T(I+1) - T(I)     
      B = (Y(I+1) - Y(I))/H - H*(Z(I+1) + 2.0*Z(I))/6.0 
      P = 0.5*Z(I) + DIFF*(Z(I+1) - Z(I))/(6.0*H) 
      P = B + DIFF*P
      SPL3 = Y(I) + DIFF*P  
      RETURN      
      END