!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  PROGRAM Intpol
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  (c) 2005, Christian B. Lang, Mathematische Methoden in der Physik
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Numerische Interpolation mit Hilfe der Formel
!  von Lagrange. Der Grad des Interpolationspolynoms
!  n ist frei waehlbar (n=np-1 bei der np-Punkt Formel).
!  Die Stuetzstellen und zu interpolierenden Funktionswerte
!  werden als Wertepaare eingegeben. Hier sind Testdaten:
!  3                  np
!  6                  ndat
!  0. 0.              daten
!  0.1 0.01           daten
!  0.2 0.04           daten
!  0.3 0.091          daten
!  0.4 0.165          daten
!  0.5 0.30           daten
!  -0.1               x low
!  0.6                x high
!  -0.1               y low
!  0.4                y high
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

PROGRAM Intpol

   implicit none
   
   integer, parameter::  datamax=100
   real, dimension(0:datamax):: x,y
   real:: xa,xb,ya,yb,xx,yy,xdel
   integer:: np,ndat,i
   character (len=256)::gnucmd,gnudat1,gnudat2
   
   write(*,*) 'np-Punkt Interpolation: np(=Grad +1 ) = '
   read(*,*) np
   
!  Daten werden eingelesen (ndat=Anzahl der Datenpunkte)
   ndat=-1
   do while ((ndat<0).or.(ndat>datamax+1))
     write(*,*) 'Wieviele Datenpunkte : (maximal ',datamax+1,')'
     read(*,*) ndat
   end do
   ndat = ndat - 1
   write(*,*) 'Daten zeilenweise eingeben: x y '
   do i = 0,ndat
     read(*,*) x(i), y(i)
   end do
   
   write(*,*) 'Darstellungsbereich:'
   
   
   write(*,*) 'Untere Schranke fuer die x-Werte= '
   read(*,*) xa
   write(*,*) 'Obere  Schranke fuer die x-Werte= '
   read(*,*) xb
   write(*,*) 'Untere Schranke fuer die y-Werte= '
   read(*,*) ya
   write(*,*) 'Obere  Schranke fuer die y-Werte= '
   read(*,*) yb

   gnudat1='intpol_plot1.dat'
   gnudat2='intpol_plot2.dat'
   gnucmd='intpol_plot.cmd'
   
   print *,'Ich erzeuge drei Files:'
   print *,trim(gnudat1)
   print *,trim(gnudat2)
   print *,trim(gnucmd)
   print *,'Falls diese Files schon existieren, erfolgt eine Fehlermeldung.'
   open(unit=9 ,file=trim(gnudat1),status='new')
   open(unit=10,file=trim(gnudat2),status='new')
   open(unit=11,file=trim(gnucmd),status='new')
   write(11,'(a,i2,a)') 'set title "Interpolation der Ordnung ',np,'"'
   write(11,'(a)') 'set xlabel "x"'
   write(11,'(a)') 'set ylabel "y"'
   write(11,'(a,f5.1,a,f5.1,a)') 'set yrange [',ya,":",yb,']'
   write(11,'(a)') 'plot "'//trim(gnudat2)//'" using 1:2 with lines t "Interpolationsfunktion", "'&
   &                      //trim(gnudat1)//'" using 1:2 with points pt 4 ps 2 t "Daten"'
   write(11,'(a)') 'pause -1'
   write(11,'(a)') 'quit'
   close(11)
   
   xdel = (xb - xa) / 100
!  Zeichne die Daten:
   do i = 0,ndat
     write(9,*)x(i), y(i), xdel
   end do
   
!  Nun berechne die stueckweise aus Interpolationspolynomen zusammengesetzte Interpolationskurve:
   xx = xa - xdel
   do i = 0,100
      xx = xx + xdel
      call Interpol(x, y, ndat, np, xx, yy)
     write(10,*)xx,yy
   end do

   close(9)
   close(10)

!  Wenn das Unterprogramm "system" nicht verfuegbar ist, so kommentieren
!  Sie die folgende Zeile aus. Dann muessen Sie die Grafik selbst mittels
! "gnuplot intpol_plot.cmd" aufreufen.

   call system(trim('gnuplot '//trim(gnucmd)))

   print *,'Berechnung und Darstellung sind abgeschlossen.'
   
END PROGRAM Intpol

SUBROUTINE Interpol (x, y,ndat,np,xx,yy)

   implicit none
   
   integer, parameter::  datamax=100
   real, dimension(0:datamax), intent(in):: x,y
   integer, intent(in)::  ndat,np
   real, intent(in)::xx
   real, intent(out) ::yy
   integer:: d1,d2,dd,i,j,na,nb
   real:: fj

   d1 = int(0.01 + np / 2)
   d2 = np - 1 - d1
   dd = d1 + d2
   
!  Finde das gewuenschte Intervall:
   
   j=0
   do while ((xx >= x(j)).and.( j < ndat))
     j=j+1
   end do
   
!  Bestimme die Indizes der Interpolationspunkte
   
   na = 0
   
!  Weiter links geht es nicht! 

   if(j - d1 > na) na = j - d1
   nb = na + dd
!  Zu weit rechts?
   if(nb > ndat) nb = ndat
   na = nb - dd
   yy = 0
   
!  Summe ueber die Terme f_j(xx)*yy(j) 

   do j = na,nb
     fj = y(j)
     do i = na,nb
!       Berechne Gewichtsfaktor f_j(x)}
        if(i.ne.j) fj = fj * (xx - x(i)) / (x(j) - x(i))
     end do
     yy = yy + fj
   end do
   
END SUBROUTINE Interpol

