!****************************************************
!*  Program to demonstrate Lagrange interpolation   *
!*     of Function SIN(X) in double precision       *
!* ------------------------------------------------ *
!* Reference: BASIC Scientific Subroutines, Vol. II *
!* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].*
!*                                                  *
!*             F90 Version by J.-P. Moreau, Paris.  *
!*                      (www.jpmoreau.fr)           *
!* ------------------------------------------------ *
!* SAMPLE RUN:                                      *
!*                                                  *
!*  Lagrange interpolation of SIN(X):               *
!*                                                  *
!*    Input X: 0.5                                  *
!*    Order of interpolation: 4                     *
!*                                                  *
!*    SIN(X) = 0.47942552 (exact value: 0.47942554) *
!*                                                  *   
!****************************************************
PROGRAM Lagrange

! Label: 100

real*8  X(0:14), Y(0:14)
real*8  xx,yy
integer iv,n

  print *,' '
  print *,'Lagrange interpolation of SIN(X):'
  iv=14
  ! Input sine table
  !-----------------------------------------------------------------
  ! Sine table values from  Handbook of mathematical functions
  ! by M. Abramowitz and I.A. Stegun, NBS, june 1964
  !-----------------------------------------------------------------
  X(1)=0.000d0; Y(1)=0.00000000d0; X(2)=0.125d0; Y(2)=0.12467473d0
  X(3)=0.217d0; Y(3)=0.21530095d0; X(4)=0.299d0; Y(4)=0.29456472d0
  X(5)=0.376d0; Y(5)=0.36720285d0; X(6)=0.450d0; Y(6)=0.43496553d0
  X(7)=0.520d0; Y(7)=0.49688014d0; X(8)=0.589d0; Y(8)=0.55552980d0
  X(9)=0.656d0; Y(9)=0.60995199d0;X(10)=0.721d0; Y(10)=0.66013615d0
  X(11)=0.7853981634d0; Y(11)=0.7071067812d0; X(12)=0.849d0; 
  Y(12)=0.75062005d0; X(13)=0.911d0; Y(13)=0.79011709d0
  X(14)=0.972d0; Y(14)=0.82601466d0
  !-----------------------------------------------------------------

  ! Input interpolation point
100 continue
  print *,' '
  write(*,"('   Input X: ')",advance='no')
  read *,xx
  write(*,"('   Order of interpolation: ')",advance='no')
  read *, n 

  call Interpol_Lagrange(iv,n,xx,yy,X,Y)

  if (n.ne.0) then
    write(*,50) yy, dsin(xx)
    goto 100
    print *,' '
  end if
  print *,' '
50 format(/'   SIN(X) = ',f10.8,'  (exact value: ',f10.8,')')
  stop

end


!********************************************************
!*          Lagrange interpolation subroutine           *
!* ---------------------------------------------------- *
!* n is the level of the interpolation ( Ex. n=2 is     *
!* quadratic ). v is the total number of table values.  *
!* X(i), Y(i) are the coordinate table values, Y(i)     *
!* being the dependant variable. The X(i) may be arbi-  *
!* trarily spaced.  x is the interpolation point which  *
!* is assumed to be in the interval  with at least one  *
!* table value to the left, and n to the right. If this *
!* is violated, n will be set to zero. It is assumed    *
!* that the table values are in ascending X(i) order.   *
!********************************************************
Subroutine Interpol_Lagrange(iv,n,xx,yy,X,Y)  
  !Labels: 100,200,300
  integer i,iv,j,k,n
  real*8  X(0:14), Y(0:14)
  real*8  xx,yy
  real*8  XL(0:9)
  ! Check to see if interpolation point is correct
  if (xx < X(1)) goto 100 
  if (xx <= X(iv-n)) goto 200
  ! An error has been encountered
100 n=0 ; return
  ! Find the relevant table interval
200 i=0
300 i = i + 1
  if (xx > X(i)) goto 300
  i = i - 1
  ! Begin interpolation
  do j = 0, n
    XL(j)=1.d0
  end do
  yy=0.d0
  do k = 0, n
    do j = 0, n
      if (j.eq.k) goto 400
      XL(k)=XL(k)*(xx-X(j+i))/(X(i+k)-X(j+i))
400 end do
    yy=yy+XL(k)*Y(i+k)
  end do
  return
end


! End of file lagrange.f90