program lsp
use fvn
implicit none
integer,parameter :: npoints=13,deg=3
 integer :: status,i
 real(kind=8) :: xm(npoints),ym(npoints),xstep,xc,yc
 real(kind=8) :: coeff(deg+1)
 xm = (/ -3.8,-2.7,-2.2,-1.9,-1.1,-0.7,0.5,1.7,2.,2.8,3.2,3.8,4. /)
 ym = (/ -3.1,-2.,-0.9,0.8,1.8,0.4,2.1,1.8,3.2,2.8,3.9,5.2,7.5 /)
 open(2,file='fvn_lsp_double_mesure.dat')
 open(3,file='fvn_lsp_double_poly.dat')
 do i=1,npoints
    write(2,44) xm(i),ym(i)
 end do
 close(2)
 call fvn_lspoly(npoints,xm,ym,deg,coeff,status)
 xstep=(xm(npoints)-xm(1))/1000.
 do i=1,1000
    xc=xm(1)+(i-1)*xstep
    yc=poly(xc,coeff)
    write(3,44) xc,yc
 end do
 close(3)
write(*,*) "All done, plot results with gnuplot using command :"
write(*,*) "pl 'fvn_lsp_double_mesure.dat' u 1:2 w p,'fvn_lsp_double_poly.dat' u 1:2 w l"
44      FORMAT(4(1X,1PE22.14))
contains
function poly(x,coeff)
    implicit none
    real(8) :: x
    real(8) :: coeff(deg+1)
    real(8) :: poly
    integer :: i
    poly=0.
    do i=1,deg+1
         poly=poly+coeff(i)*x**(i-1)
    end do
end function
end program