!*********************************************************************
! driver program for dspsp20.f90
! date       : october 25, 2003
! programmer : makoto nakajima
!
! driver program for dspsp (verison 2.0).
! given a function "ridicurous", which is differentiable over [0,1],
! shows the true values of the function and the values obtained by
! using schumaker shape-preserving quadratic spline, using 10 points.
! there are two versions of interpolation. one uses also derivative
! of the function and the other not. the results from both methods
! are shown.
! NOTICE THE DIFFERENCE BETWEEN WHEN THE DERIVATIVE INFO IS USED AND
! WHEN IT IS NOT USED
!(i) dspspline is the one without derivatives.
!(ii) dspspline2 is the one where you provide with derivatives
!        of two end points. I made it before.
!(iii) dspspline3 is the one for which you provide derivatives at
!        all the grid points.
!
!
!There are some other subroutines as well.  dspspleasyd
!is basically the one without derivatives (dspspline), but the 
!subroutine also returns the deirvative at the evaluation point.



!********************************************************************
program driver_dspsp20

  !******* activate module *******
  use dspsp
  implicit none

  !******* parameters *******
  integer,parameter:: prec=8
  !double precision
  integer,parameter:: ngrid=10
  !number of grid points used for interpolation
  integer,parameter:: neval=100
  !number of points for function evaluation
  real(prec),parameter:: xub=1.0,xlb=0.0
  !upperbound and lowerbound of domain
  real(prec),parameter:: xstep=2.0
  !controlling the diustance of grid points
  real(prec),parameter:: xd=1.0e-10
  !small step to calculate numerical derivative

  !******* variables *******
  real(prec),dimension(ngrid):: sx,sy,sd
  real(prec),dimension(neval):: x,y,y1,y2
  real(prec):: yd
  integer:: p0,ierr

  !******* set evaluation points *******
  call set_grids(xlb,xub,1.0_prec,neval,x)

  !******* set grid points *******
  call set_grids(xlb,xub,xstep,ngrid,sx)

  !******* calculate the value at evaluation points *******
  do p0=1,neval
     y(p0)=ridiculous(x(p0))
  end do

  !******* calculate values and drivatives at grid points *******
  do p0=1,ngrid-1
     sy(p0)=ridiculous(sx(p0))
     yd=ridiculous(sx(p0)+xd)
     sd(p0)=(yd-sy(p0))/xd
  end do
  p0=ngrid
  sy(p0)=ridiculous(sx(p0))
  yd=ridiculous(sx(p0)-xd)
  sd(p0)=(sy(p0)-yd)/xd

  !******* calculate interpolated values *******
  do p0=1,neval
     call dspspleasy(ngrid,sx,sy,x(p0),y1(p0),ierr)
     if (ierr>0) then
        print 100
        stop
     end if
  end do
  do p0=1,neval
     call dspspleasy3(ngrid,sx,sy,x(p0),sd,y2(p0),ierr)
     if (ierr>0) then
        print 102
        stop
     end if
  end do

  !******* write out results *******
  print 110
  print 112
  do p0=1,neval
     print 114,x(p0),y(p0),y1(p0),y2(p0)
  end do

  !******* end of the program *******

100 format('something wrong with dspspleasy...')
102 format('something wrong with dspspleasy3...')
110 format('******* results *******')
112 format('      x',1x,'      true',1x,' sp_spline',1x,' sp_sp_der')
114 format(f7.4,1x,f10.4,1x,f10.4,1x,f10.4)

contains

  !********************************************************************
  ! function: ridiculous
  !********************************************************************
  function ridiculous(x)
    implicit none
    integer,parameter:: prec=8
    real(prec):: x,ridiculous
    ridiculous=(x+1.0)**(1.0_prec-2.0_prec)/(1.0_prec-2.0_prec)
    return
  end function ridiculous

  !********************************************************************
  ! subroutine: set_grids
  !********************************************************************
  subroutine set_grids(xtop,xend,gpara,ngrid,gridval)
    implicit none
    integer,parameter:: prec=8
    real(prec),intent(in):: xtop,xend,gpara
    integer,intent(in):: ngrid
    real(prec),intent(out),dimension(ngrid):: gridval
    integer:: grid0

    gridval(1)=xtop
    gridval(ngrid)=xend
    do grid0=2, ngrid-1
       gridval(grid0)=gridval(1)+(gridval(ngrid)-gridval(1))*&
            (real(grid0-1,prec)/real(ngrid-1,prec))**gpara
    end do
    return
  end subroutine set_grids

end program driver_dspsp20