program akima

   use fvn_interpol

   implicit none
   integer :: nbpoints,nppoints,i
   real(8),dimension(:),allocatable :: x_d,y_d,breakpoints_d
   real(8),dimension(:,:),allocatable :: coeff_fvn_d
   real(8) :: 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./dfloat(nbpoints)
   do i=1,nbpoints
      x_d(i)=-10.+dfloat(i)*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./dfloat(nppoints)
   do i=1,nppoints
      xp_d=-11.+dfloat(i)*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