program akima

   use fvn_interpol

   implicit none
   integer :: nbpoints,nppoints,i

   real(kind=dp_kind),dimension(:),allocatable :: x_d,y_d,breakpoints_d
   real(kind=dp_kind),dimension(:,:),allocatable :: coeff_fvn_d
   real(kind=dp_kind) :: xstep_d,xp_d,ty_d,fvn_y_d

   open(2,file='fvn_akima_double.dat')
   open(3,file='fvn_akima_breakpoints_double.dat')
   nbpoints=30
   allocate(x_d(nbpoints))
   allocate(y_d(nbpoints))
   allocate(breakpoints_d(nbpoints))
   allocate(coeff_fvn_d(4,nbpoints))

   xstep_d=20./real(nbpoints,dp_kind)

   do i=1,nbpoints

      x_d(i)=-10.+real(i,dp_kind)*xstep_d

      y_d(i)=dsin(x_d(i))
      write(3,44) x_d(i),y_d(i)
   end do
   close(3)
   call fvn_akima(nbpoints,x_d,y_d,breakpoints_d,coeff_fvn_d)
   nppoints=1000

   xstep_d=22./real(nppoints,dp_kind)

   do i=1,nppoints

      xp_d=-11.+real(i,dp_kind)*xstep_d

      ty_d=dsin(xp_d)
      fvn_y_d=fvn_spline_eval(xp_d,nbpoints-1,breakpoints_d,coeff_fvn_d)
      write(2,44) xp_d,ty_d,fvn_y_d
   end do
   close(2)
  deallocate(coeff_fvn_d,breakpoints_d,y_d,x_d)

  write(6,*) "All done, plot results with gnuplot using command :"
  write(6,*) "pl 'fvn_akima_double.dat' u 1:2 w l,'fvn_akima_breakpoints_double.dat' w p"

  
  44       FORMAT(4(1X,1PE22.14))
  end program