program gls
use fvn_misc
use excursion
implicit none

    interface
        subroutine zef(m, n, x, fvec, iflag)
            integer(4), intent(in) :: m,n
            real(8), dimension(:), intent(in) :: x
            real(8), dimension(:), intent(inout) :: fvec
            integer(4), intent(inout) :: iflag
        end subroutine
    end interface


 integer,parameter :: npoints=13,deg=3
 integer :: status,i
 real(kind=dp_kind) :: xstep,xc,yc
 real(kind=dp_kind) :: coeff(deg+1)
 real(kind=dp_kind) :: fvec(npoints)
 integer(4) :: iwa(deg+1)
 integer(4) :: info
 real(8) :: tol
 
 allocate(xm(npoints),ym(npoints))
 
 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_lm.dat')
 do i=1,npoints
    write(2,44) xm(i),ym(i)
 end do
 close(2)
 
 coeff=1.
 call fvn_lm(zef,npoints,deg+1,coeff,info)
 write(*,*) "info : ",info
 
 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_lm.dat' u 1:2 w l"
44      FORMAT(4(1X,1PE22.14))
contains
function poly(x,coeff)
    use fvn_common
    implicit none
    real(kind=dp_kind) :: x
    real(kind=dp_kind) :: coeff(deg+1)
    real(kind=dp_kind) :: poly
    integer :: i
    poly=0.
    do i=1,deg+1
         poly=poly+coeff(i)*x**(i-1)
    end do
end function
end program

subroutine zef(m, n, x, fvec, iflag)
    use excursion
    implicit none
    integer(4), intent(in) :: m,n
    real(8), dimension(:), intent(in) :: x
    real(8), dimension(:), intent(inout) :: fvec
    integer(4), intent(inout) :: iflag
    
    integer(4) :: i
        
    do i=1,m
        fvec(i)=ym(i)-(x(1)+x(2)*xm(i)+x(3)*xm(i)**2+x(4)*xm(i)**3)
    enddo
end subroutine