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