      program Lagrange 
c     
c     stosujac metode Langrange'a program                       
c     wykonuje interpolacje funkcji, ktorej wezly czyta z dysku
c      [korzystamy ze zrodel podprogramu LTRP z biblioteki NSWC]
c      autor: Janusz Mrozek,     8.11.2007  
c

      parameter (maxwezl=20,lzag=10,maxint=(maxwezl-1)*lzag+1)

c
c       maxwezl = maks. liczba wezlow interpolacji
c       lzag    = "zageszczenie" tablicowania na uzytek wykresu
c       maxint  = maks.liczba wartosci wielomianu interpolacyjnego 
c                 przeznaczonych do sporzadzenia wykresu
c
      real*4 Xwezl(maxwezl),Ywezl(maxwezl),Xintp(maxint),Yintp(maxint),
     %       Wspwlm(maxwezl),
     %       xpocz,xkonc,xkrok
      character*80 nazwa_wezl, nazwa_pwykr
      integer i,j,Lwezl,Lprzedz,Lwint,Lwykr, Ldod
c
c     Lwezl   = liczba wezlow interpolacji 
c                czytanych ponizej z urz. nr 10
c     Lprzedz = liczba przedzialow pomiedzy wezlami
c                malo odkrywcze: Lprzedz = Lwezl-1
c     Lwint   = liczba wezlow uzywanych do interpolacji Lagrange'a
c                  od wartosci 
c               min = 2 (interpolacja odcinkami prostej lamanej) do
c                  do wartosci
c               max = Lwezl (interpolacja jednym wielomianem
c               przechodzacym przez wszystkie wezly (stopien = Lwezl-1)
c
c     Lwykr   = liczba wartosci interpolowanych gromadzona na uzytek
c               wykresu zapisywana na dysku poprzez urzadzenie nr 11
c
      write(*,*) ' jak nazywa sie plik z wezlami interpolacji?'
      read (*,'(a)') nazwa_wezl
c
c     otwieramy plik z danymi 
c
      open(unit=10,file=nazwa_wezl,status='old')
c
c      teraz czytamy wezly z dysku
c
       Lwezl = 0
10     read(10,*,end=20) x,y    ! czytamy az do napotkania znaku 
                                  ! konca pliku

       Lwezl = Lwezl+1            ! udalo sie wczytac kolejny wezel            

       Xwezl(Lwezl) = x           ! zapamietujemy wczytana pare liczb
       Ywezl(Lwezl) = y           ! do wlasciwych miejsc w tablicach

       go to 10                   ! wracamy do rozkazu czytania 
c
c      skonczyly sie dane na dysku, jestesmy w linii z etykieta 20
c      i przechodzimy do rozmowy na temat szczegolow interpolacji

20     write(*,*) ' z dysku wczytano',Lwezl,' wezlow'

       write(*,*) ' przez ile wezlow prowadzic wielomian interpol.?'
       write(*,*) ' min = 2, maks =',Lwezl
c
       read (*,*) Lwint

c
c      Na uzytek wykresu powinnismy sporzadzic dokladna tabele wartosci
c      wielomianu interpolacyjnego, dzielac kazdy przedzial miedzy dwoma
c      kolejnymi wezlami na Ldod dodatkowych podprzedzialow
c      gdzie    2 <= Ldod <= lzag 
c

30     write(*,*) ' na ile podprzedzialow podzielic odcinki pomiedzy'
       write(*,*) ' kolejnymi wezlami?'
       write(*,*) ' zakres wartosci to 2 do ',lzag

       read (*,*) Ldod

       if(Ldod .lt. 2 .or. Ldod .gt. lzag) then  ! bywamy podejrzliwi...
                                                                     ! i slusznie, jak sie okazalo
                                                                     ! klient wyraznie spi lub 
                                                                     ! nie co czyni
        write(*,*) 'bledna wartosc! trzeba wpisac jeszcze raz!'
        go to 30
       end if
c
c      wiemy juz wszystko co trzeba by przygotowac tabelke wartosci 
c      Xintp zmiennej niezaleznej x, dla ktorych wyznaczymy wartosci 
c      wielomianu interpolacyjnego - na podstawie tej tabeli mozemy
c      zbudowac wykres wielomianu. Wykres ten uzyjemy do porownania 
c      oryginalnej funkcji 
c
       Lprzedz = Lwezl-1   ! zaczynamy od obliczenia ile jest przedzialow
                           ! pomiedzy wezlami

c
c      teraz zaczynamy wycieczke po tych przedzialach
c                           
       Lwykr = 1
       do i=1,Lprzedz

        xpocz = xwezl(i)     ! ustalamy gdzie sie przedzial zaczyna
        xkonc = xwezl(i+1)   ! i gdzie sie konczy

        xkrok = (xkonc - xpocz)/Ldod ! obliczamy krok wewnatrz i-tego przedzialu

        Xintp(Lwykr) = xpocz  ! wpisujemy wartosc x dla poczatku przedzialu

        do j= 1,Ldod-1
          xpocz = xpocz+xkrok  ! oszukujemy "przesuwajac" poczatek
          Lwykr = Lwykr+1       ! przedzialu o xkrok - podczas pracy 
                                            !  w poprzednim przedziale wpisano juz 
                                            ! wyniki dla punktu x, ktory konczy poprzedni
                                            ! przedzial i stanowi zarazem pocztek biezacego.

         Xintp(Lwykr) = xpocz  ! zapisujemy kolejna wartosc x do wykresu  
        end do

          Lwykr = Lwykr+1      ! przesuwamy wskaznik w tablicy Xintp o 1        
       end do                  ! na uzytek nastepnego przedzialu

          Xintp(Lwykr) = xkonc ! nasz genialny algorytm gubi ostatnia
                               ! wartosc, stad trzeba ja wpisac "recami"
                               ! poza petla       
c
c      wszystko juz wiemy, mamy komplet danych
c      teraz mozemy zlecic robote procedurze LTRP z biblioteki NSWC
c                                         (a Kali proznowac i jesc!)
c

        write(*,*) 'Lwykr=',Lwykr,'Lprzedz=', Lprzedz,' Ldod=',Ldod

      call LTRP(Lwint,Xwezl,Ywezl,Lwezl,Xintp,Yintp,Lwykr,Wspwlm,IERR)
c
c      zawsze sprawdzamy czy procedura biblioteczna nie zglosila sygnalu
c      o klopotach
c       
      if(IERR.ne.0) then 
c       jezeli tu jestesmy to, niestety, zrealizowal sie czarny scenariusz:
c      nasz podopieczny cos  zdolal spaprac...

        write(*,*) 'LTRP zglosila blad nr',ierr  ! piszemy ze jest zle
        stop ' poczytaj dokumentacje i pomysl!'  ! i mordujemy program,
                                                                        ! pamietajac o misji
                                                                        ! edukacyjnej uczelni
      end if

c
c     sytuacja jest dobra, nie bylo sygnalu bledu z LTRP
c                    mamy wyniki!!!
c      
c        pytanie: jak by je tu nazwac?
c
      write(*,*) ' jak nazwac plik z wynikami?'
      read (*,'(a)') nazwa_pwykr
c
c     otwieramy plik nie przejmujac sie nadmiernie tym czy jego nazwa
c     juz byla uzywana (status="unknown')
c     w koncu sie uczymy i mamy prawo bladzic...
c     w najgorszym razie wymarzemy sobie wyniki w pliku o tej
c     samej nazwie
c
      open(unit=11,file=nazwa_pwykr,status='unknown')
c
c     teraz mozemy wziac sie do zapisania naszych wynikow na dysku:
c
      do i=1,Lwykr                    ! dla kazdego punktu na wykresie

       write(11,*) Xintp(i),Yintp(i)   ! zapisujemy na dysku wartosci x 
                                                        !                i obliczone dla nich   y 
                                                        ! z wielomianu interpolacyjnego
      end do
c
c       pozostaje jeszcze tylko zamkniecie pliku...
c

      close(unit=11,status='keep')
c
c       i podanie klientowi szokujacej informacji, ze to juz koniec     
c
      stop 'skonczone!'
c
      end

      SUBROUTINE LTRP (M,X,Y,N,XI,YI,NI,T,IERR)
C     ******************************************************************
C     LAGRANGE INTERPOLATION
C     ******************************************************************
      DIMENSION X(N),Y(N),XI(NI),YI(NI),T(M)
C
C     CHECK INPUT
C
      IF (M.LT.2) GO TO 130
      IF (M.GT.N) GO TO 131
      IF (NI.LT.1) GO TO 132
      IERR = 0
C
C     INITIALIZATION
C
      MM1 = M-1
      K  = 1
      XX = XI(1)
      ILOLD = 0
C
C     FIND THE SUBINTERVAL WHICH CONTAINS XX. I = 1 IF XX.LT.X(2)
C     AND I = N IF XX.GE.X(N). OTHERWISE X(I).LE.XX.LT.X(I+1).
C
      IF (XX-X(1)) 10,11,20
   10 I = 1
      IL = 1
      IR = M
      GO TO 80
   11 I = 1
      YI(K) = Y(1)
      GO TO 110
C
   20 IF (X(N)-XX) 21,22,23
   21 I = N
      IL = N-M+1
      IR = N
      GO TO 80
   22 I = N
      YI(K) = Y(N)
      GO TO 110
   23 IL = 1
      IR = N
C
C     BISECTION SEARCH
C
   30 I = (IL+IR)/2
      IF (I.EQ.IL) GO TO 50
      IF (XX-X(I)) 31,32,33
   31 IR = I
      GO TO 30
   32 YI(K) = Y(I)
      GO TO 110
   33 IL = I
      GO TO 30
C
C     LINEAR FORWARD SEARCH
C
   40 IF (XX-X(I+1)) 50,41,42
   41 I = I+1
      YI(K) = Y(I)
      GO TO 110
   42 I = I+1
      GO TO 40
C
C     POINT XX LIES IN THE OPEN SUBINTERVAL (X(I),X(I+1)).
C     FIND THE M CLOSEST POINTS X(IL),...,X(IR) TO XX.
C
   50 IF (I.GT.M) GO TO 51
      IL = 1
      NUM = I
      GO TO 60
   51 IL = I-M+1
      NUM = M
C
   60 IPM = I+M
      IF (IPM.LE.N) GO TO 61
      IR = N
      NUM = NUM+N-I
      GO TO 70
   61 IR = IPM
      NUM = NUM+M
C
   70 NUM = NUM-M
      IF (NUM.EQ.0) GO TO 80
      DL = XX-X(IL)
      DR = X(IR)-XX
      DO 72 L=1,NUM
      IF (DL.LE.DR) GO TO 71
      IL = IL+1
      DL = XX-X(IL)
      GO TO 72
   71 IR = IR-1
      DR = X(IR)-XX
   72 CONTINUE
C
C     COMPUTE THE COEFFICIENTS T(1),...,T(M) OF THE BACKWARD
C     NEWTON FORM OF THE INTERPOLATING POLYNOMIAL.
C
   80 IF (IL.EQ.ILOLD) GO TO 100
      ILOLD = IL
      ILM1 = IL-1
      DO 81 J=1,M
      L = ILM1+J
   81 T(J) = Y(L)
C
      DO 91 ISTEP=1,MM1
      JMAX = M-ISTEP
      DO 90 J=1,JMAX
      II = ILM1+J
      L = II+ISTEP
   90 T(J) = (T(J)-T(J+1))/(X(II)-X(L))
   91 CONTINUE
C
C     EVALUATION OF THE INTERPOLATING POLYNOMIAL
C
  100 YI(K) = T(1)
      DO 101 J=2,M
      L = ILM1+J
  101 YI(K) = T(J)+YI(K)*(XX-X(L))
C
C     NEXT POINT
C
  110 IF (K.GE.NI) RETURN
      K = K+1
      XX = XI(K)
      IF (XX-X(1)) 10,11,111
  111 IF (X(N)-XX) 21,22,120
C
  120 IF (XX-XI(K-1)) 121,122,40
  121 IL = 1
      IR = MIN0(I+1,N)
      GO TO 30
  122 YI(K) = YI(K-1)
      GO TO 110
C
C     ERROR RETURN
C
  130 IERR = 1
      RETURN
  131 IERR = 2
      RETURN
  132 IERR = 3
      RETURN
      END


