program laginterp 

    implicit none
    real (kind=8), dimension(0:10) :: xpt, ypt
    real (kind=8) ::  x, y, dx
    integer :: i, errval, cntr

    interface
        function interp ( N, xpt, ypt, x ) result(y)
                real (kind=8) :: xpt(0:*), ypt(0:*), x, y
        end function interp
    end interface

    open (unit=15, file="testdata.dat", status="old")

    errval = 0
    cntr = 0
    do while ( errval /= -1 )
        read ( 15, *, iostat=errval ) x, y
        if ( errval == 0 ) then
            xpt(cntr) = x
            ypt(cntr) = y
            cntr = cntr + 1 
        else
            close(15)
        endif
    enddo
    print *, "Read ", cntr, " data points."

    dx = (maxval(xpt) - minval(xpt)) / 1000.0 
    x =  minval(xpt) 
    open (unit=15, file="interpolated.dat", status="unknown")
    do while (x <= maxval(xpt) ) 
       write (15, * ) x, interp(cntr, xpt, ypt, x)
       x = x + dx
    enddo
    close(15) 

end program laginterp


function interp ( N, xpt, ypt, x ) result(y)

        real (kind=8) :: xpt(0:*), ypt(0:*), x, y
        real (kind=8) :: buildsum, prod 

        buildsum = 0.0D0 
        do k = 0, N-1
            prod =1.0D0
            do i = 0, N-1
                if (i .ne. k) prod = prod * (x-xpt(i))/(xpt(k)-xpt(i))
            enddo
            buildsum = buildsum + ypt(k) * prod
        enddo
        y = buildsum 
        return
end function interp