!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   Natural Cubic Spline Interpolation  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!  A cubic spline function is constructed that        !!
!!  interpolates a given function at given intervals.  !!
!!  A further restriction for natural cubic splines    !!
!!  is that the second derivative of the function at   !!
!!  the end points is set to zero, allowing for a      !!
!!  variable slope at the end points but also          !!
!!  stabilizing the curvature at the end points.       !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!  compiled with gfortran -fdefault-real-8 nat_spline.f95  !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  ABW 2012-2015  !
!!!!!!!!!!!!!!!!!!!
program natspline
  implicit none
  
  integer :: i
  real(kind=8) :: h
  integer,parameter :: n = 7  !  Number of interpolating cubic splines
  real(kind=8),parameter :: pi = 3.1415926535897932  !  pi to 16 digits
  real(kind=8),parameter :: a = -2.0*pi, b = 2.0*pi  !  limits of interpolation, [a,b]
  real(kind=8),dimension(0:n) :: xx,yy
  real(kind=8),dimension(0:n-1) :: hi,g,at,bt,ct
  real(kind=8),dimension(0:n) :: zans
  
  !  Write function to interpolate
  write(*,*) 'Function to interpolate:  x^2*sin(x)'
  write(*,*)
  
  !  step size
  h = (b - a)/real(n)
  
  !  taking (n+1) points of data
  write(*,*) '      X','            Y'
  do i = 0,n
     xx(i) = a + real(i)*h
     yy(i) = f(xx(i))
     write(*,"(F13.9,F13.9)") xx(i),yy(i)
  enddo
  write(*,*)
  
  !  writing parts of matrix equation
  do i = 0,n-1
     hi(i) = xx(i+1) - xx(i)
     g(i) = 6.0*(yy(i+1) - yy(i))/hi(i)
  enddo
  
  !  simplified Gaussian Elimination to solve tridiagonal matrix equation
  call gelim(n,hi,g,zans)
  
  !  z(0) = 0 = z(n) for natural spline
  !  Printing out z(0) to z(n)
  do i = 0,n
     if (i < 10) then
        write(*,"(A,I1.1,A,F20.16,F20.16)") 'z(',i,')  = ',zans(i)
     else
        write(*,"(A,I2.1,A,F20.16,F20.16)") 'z(',i,') = ',zans(i)
     endif
  enddo
  write(*,*)

  !  Coefficients of S(x)_i
  do i = 0,n-1
     at(i) = (1.0/(6.0*hi(i)))*(zans(i+1) - zans(i))
     bt(i) = zans(i)/2.0
     ct(i) = (-hi(i)/6.0)*zans(i+1) - (hi(i)/3.0)*zans(i) + (1.0/hi(i))*(yy(i+1) - yy(i))  
  enddo
  
  !  used to print interpolating spline S(x)_i
  !  plot S(x)_i with this or print data with interS() function
  write(*,*) 'Interpolating Cubic Spline:'
  write(*,*)
  do i = 0,n-1
     write(*,110) 'S(x)_',i,' = ',yy(i),' + (x -',xx(i),')*{',ct(i),' + (x -',xx(i), &
	  ')*[',bt(i),' + (x -',xx(i),')*',at(i),']} for',xx(i),' < x <',xx(i+1)
110 format (A,I2,A,F8.4,A,F8.4,A,F8.4,A,F8.4,A,F8.4,A,F8.4,A,F8.4,A,F8.4,A,F8.4)
  enddo
  
  !  Subroutine to write data from S(x)_i and f(x)
  call writedata(a,h)
  
contains
  
  !  Function to interpolate
  function f(x)
    implicit none
    
    real(kind=8) :: f,x
    
    f = x*x*sin(x)
  end function f
  
  !  Simplified Gaussian Elimination to solve [hi D hi]*zans = Bvec
  subroutine gelim(n,hi,g,zans)
    implicit none
    
    integer :: i
    integer,intent(in) :: n
    real(kind=8),dimension(1:n-1) :: Bvec,D
    real(kind=8),dimension(0:n-1) :: hi,g
    real(kind=8),dimension(0:n),intent(inout) :: zans
    
    D(1) = 2.0*(hi(0) + hi(1))  !  Diagonal terms
    Bvec(1) = g(1) - g(0)  !  right side of equation
    do i = 2, n-1
       D(i) = 2.0*(hi(i) + hi(i-1)) - (hi(i-1)**2)/(D(i-1))
       Bvec(i) = g(i) - g(i-1) - (hi(i-1)*Bvec(i-1))/(D(i-1))
    enddo
  
    zans(n) = 0.0  !  For natural cubic spline
    do i = n-1,1,-1
       zans(i) = (Bvec(i) - hi(i)*zans(i+1))/(D(i))
    enddo
    zans(0) = 0.0  !  For natural cubic spline
    
  end subroutine gelim  

  !  Function containing the interpolating cubic splines. Use to print data from spline.
  !  the i^th cubic spline is in the range [xx(i), xx(i+1)]
  function interS(x,i2)
    implicit none
    
    integer :: i2
    real(kind=8) :: interS,x
    
    interS = yy(i2) + (x - xx(i2))*(ct(i2) + (x - xx(i2))*(bt(i2) + (x - xx(i2))*at(i2)))
  end function interS 
  
  !  Write data from interpolating cubic spline S(x)_i at intervals [x(i),x(i+1)]
  subroutine writedata(a,h)
    implicit none
    
    integer :: i
    real(kind=8) :: xt,a,h
    
    open (unit=10,file="interS_data_n7.dat",action="write")
    open (unit=11,file="function_data.dat",action="write")
    do i = 0,100
       xt = a + h*(real(i)/100.0)
       write(10,*) xt,interS(xt,0)
       write(11,*) xt,f(xt)
    enddo
    do i = 101,200
       xt = a + h*(real(i)/100.0)
       write(10,*) xt,interS(xt,1)
       write(11,*) xt,f(xt)
    enddo 
    do i = 201,300
       xt = a + h*(real(i)/100.0)
       write(10,*) xt,interS(xt,2)
       write(11,*) xt,f(xt)
    enddo  
    do i = 301,400
       xt = a + h*(real(i)/100.0)
       write(10,*) xt,interS(xt,3)
       write(11,*) xt,f(xt)
    enddo  
    do i = 401,500
       xt = a + h*(real(i)/100.0)
       write(10,*) xt,interS(xt,4)
       write(11,*) xt,f(xt)
    enddo  
    do i = 501,600
       xt = a + h*(real(i)/100.0)
       write(10,*) xt,interS(xt,5)
       write(11,*) xt,f(xt)
    enddo  
    do i = 601,700
       xt = a + h*(real(i)/100.0)
       write(10,*) xt,interS(xt,6)
       write(11,*) xt,f(xt)
    enddo   
    close(10)
    close(11)
  end subroutine writedata
  
end program natspline
