program akima
 use fvn_interpol
use Kind_Definition
 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