C======================================================================
C---  Homework #2 : Drawing window and graphics using X programming
C
C---  Programmed by Sung-Hoon Ko
C
C----------------------------------------------------------------------
C---  Main program of Polynomial interpolation
C---  
C---  r:input value
C---  a0,a1,b,c : constant variable
C---  f0,f1 : exact value
C---  index : index for data points 
C---  x1(),x2(),x3() : array which store data entries
C---  y10(),y20(),y30() : array which store value when a=0, b=c=0.1
C---  y11(),y21(),y31() : array which store value when a=1, b=c=0.1
C---  pol_val() : array which store values caculated by interpolation
C---  dTIME : library function to time
C---  tarray(), time: variable for time checking

      PROGRAM Hw2_Polynomial_Interpolation

      INCLUDE 'hw2.h'
      REAL r,index,graph_cntl,time_cntl
      REAL f0,f1
      REAL pol_val(6)
      INTEGER i,j,length

C---  Read inpur from terminal. 
  88  WRITE(6,7) '( To terminate enter 999 )'              
      WRITE(6,*) '--->>> Input r (0.1<r<4.0)'
      READ *, r
     
C---  if r==999 then teminate.      
      if(r.EQ.999) then 
        goto 999
      endif

      WRITE(6,*) '-->> Type "1" to draw a graph 
     *                  otherwise type any number'
      READ *, graph_cntl

      WRITE(6,*) '-->> Type "1" to time the routine 
     *                  otherwise type any number'
      READ *, time_cntl

C---  Check the input range which is 0.1< <4.0.
      if((r.LE.0.1) .OR. (r.GE.4.0)) then
        WRITE(6,7) '-- r should be 0.1< <4.0 try again'
  7     FORMAT(//,A,/)
        goto 88
       endif
         
C---  Calculate exact value
      CALL EXACT_VALUE(r,a0,b,c,f0)     
      CALL EXACT_VALUE(r,a1,b,c,f1)
      
C --- Set index value and compute value for 101 entries    
      index=(4.0-0.1)/100.0      
      DO i=0,100
        x1(i+1)=0.1+i*index
        CALL EXACT_VALUE(x1(i+1),a0,b,c,y10(i+1))
        CALL EXACT_VALUE(x1(i+1),a1,b,c,y11(i+1))
      ENDDO

C --- Set index value and compute value for 1001 entries
      index=(4.0-0.1)/1000.0
      DO i=0,1000
        x2(i+1)=0.1+i*index
        CALL EXACT_VALUE(x2(i+1),a0,b,c,y20(i+1))
        CALL EXACT_VALUE(x2(i+1),a1,b,c,y21(i+1))
      ENDDO

C --- Set index value and compute value for 10001 entries
      index=(4.0-0.1)/10000.0
      DO i=0,10000
        x3(i+1)=0.1+i*index
        CALL EXACT_VALUE(x3(i+1),a0,b,c,y30(i+1))
        CALL EXACT_VALUE(x3(i+1),a1,b,c,y31(i+1))
      ENDDO

C---  Calculate by 1st and 2nd order interpolation
C---  length : number of data
C---  length=2 : 1st order interpolation
C---  length=3 : 2nd order interpolation

      DO length=2,3

C---  101 entries      
        CALL LOCATE(x1,entry1,r,j)                
        CALL POLINT(x1(j),y10(j),length,r,pol_val(1))
        CALL POLINT(x1(j),y11(j),length,r,pol_val(2))

C---  1001 entries
        CALL LOCATE(x2,entry2,r,j) 
        CALL POLINT(x2(j),y20(j),length,r,pol_val(3))
        CALL POLINT(x2(j),y21(j),length,r,pol_val(4))  
        
C---  10001 entries
        CALL LOCATE(x3,entry3,r,j)
        CALL POLINT(x3(j),y30(j),length,r,pol_val(5))
        CALL POLINT(x3(j),y31(j),length,r,pol_val(6))          
      
C---  Print out the results        
        WRITE(6,22) 'Input: r =',r 
 
        WRITE(6,33) '<Case 1.  a=0,  b=c=0.1>'
        WRITE(6,55) 'Exact Value          = ',f0
        if(length .EQ. 2) then
           WRITE(6,55) ' (1st order interpolation)'
        else
           WRITE(6,55) ' (2nd order interpolation)'
        endif
        WRITE(6,55) 'Interpolation(101)   = ',pol_val(1)
        WRITE(6,55) 'Interpolation(1001)  = ',pol_val(3)
        WRITE(6,55) 'Interpolation(10001) = ',pol_val(5)
        WRITE(6,33) 'Exact - Interpolation'
        WRITE(6,55) 'Exact-Inter(101)   = ',f0-pol_val(1)
        WRITE(6,55) 'Exact-Inter(1001)  = ',f0-pol_val(3)
        WRITE(6,55) 'Exact-Inter(10001) = ',f0-pol_val(5)

        WRITE(6,33) '<Case 2.  a=1,  b=c=0.1>'
        WRITE(6,55) 'Exact Value          = ',f1
        if(length .EQ. 2) then
           WRITE(6,55) ' (1st order interpolation)'
        else
           WRITE(6,55) ' (2nd order interpolation)'
        endif
        WRITE(6,55) 'Interpolation(101)   = ',pol_val(2)
        WRITE(6,55) 'Interpolation(1001)  = ',pol_val(4)
        WRITE(6,55) 'Interpolation(10001) = ',pol_val(6)
        WRITE(6,33) 'Exact - Interpolation'
        WRITE(6,55) 'Exact-Inter(101)   = ',f1-pol_val(2)
        WRITE(6,55) 'Exact-Inter(1001)  = ',f1-pol_val(4)
        WRITE(6,55) 'Exact-Inter(10001) = ',f1-pol_val(6)

  22    FORMAT(// A,F10.5 )
  33    FORMAT(// 1X, A /)
  55    FORMAT(3X,A,E17.10)

      ENDDO


c      do i=0,entry1
c        write(6,*) 'X',i, '=', x1(i) , y11(i)
c      enddo
 

C---  Draw the graph
       if(graph_cntl .EQ. 1.) then     ! Draw a graph?
          CALL  dp_initialize
          CALL  dp_draw_graph(x3,y31,entry3)
       endif
        

C---  Time the routines
       if(time_cntl .EQ. 1.) then       
         WRITE(6,33) '[ The running time ]'  
         CALL Exact_Time(r,a0,b,c,f0)
         CALL Inter_Time(x1,y10,length,r)
       endif
 
      if(r.NE.999) goto 88
      
  999 STOP
      END

     
      
C----------------------------------------------------------
C---  SUBROUTINE  EXACT_VALUE()
C---  calculate exact value for f/r

      SUBROUTINE  EXACT_VALUE(rr,aa,bb,cc,yy)
      REAL rr,aa,bb,cc,yy
      
      yy=aa/(rr**3)-(6.0*bb)/(rr**8)+(12.0*cc)/(rr**14)

      RETURN
      END



C----------------------------------------------------------
C---  SUBROUTINE LOCATE()
C---  Given an array xx(1:n), and given a value jj
C---  such that x is between xx(j) and xx(j+1).
C---  xx(1:n) must be monotonic, either increasing or 
C---  decreasing. j=0 or j=n is returned to indicate
C---  that x is out of range.

      SUBROUTINE LOCATE(xx,n,x,jj)
      INTEGER jj,n
      REAL x,xx(n)
      INTEGER jl,jm,ju
C---  Initialize lower and upper limits
      jl=0
      ju=n+1
C---  If we are not yet done
   10 if(ju-jl .gt. 1) then
C---  Compute a midpoint
        jm=(ju+jl)/2
        if((xx(n) .gt. xx(1)) .eqv. (x .gt. xx(jm))) then
C---  Replace the lower limit
           jl=jm
        else
C---  Replace the upper limit
           ju=jm
        endif
      goto 10
      endif
C---  Set the output
      jj=jl
      RETURN
      END



C-------------------------------------------------------------
C---  SUBROUTINE  POLINT()
C---  Given  arrays xa and ya, each of length n, and given
C---  a value x, this routine returns a value y.
C---  If P(x) is the polynomial of degree N-1 such that
C---  P(XAi)=YAi, i=1,...,n, then the returned value y=P(x).

      SUBROUTINE  POLINT(xa,ya,n,x,y)
      INTEGER n,NMAX
      REAL dy,x,y,xa(n),ya(n)
C---  Largest anticipated value of n.
      PARAMETER (NMAX=10)
      INTEGER i,m,ns
      REAL den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)     
      ns=1
      dif=abs(x-xa(1))
C---  Here we find the index ns of the closest table entry.
      do i=1,n
        dift=abs(x-xa(i))
        if(dift .lt. dif) then
          ns=i
          dif=dift
        endif
C---  Initialize the tableau of c's and d's
        c(i)=ya(i)
        d(i)=ya(i)
      enddo 
C---  This is the initial approximation to y.
      y=ya(ns)
      ns=ns-1
C---  For each column of the tableau, we loop over
C---  the current c's and d's and update them.
      do m=1,n-1
        do i=1,n-m
          ho=xa(i)-x
          hp=xa(i+m)-x
          w=c(i+1)-d(i)
          den=ho-hp
          if(den .eq. 0.) pause 'failure in polint'
C---      This error can occur only if two input xa's are
C---      (to within roundoff) identical.
          den=w/den
C---      Here the c's and d's are updated.
          d(i)=hp*den
          c(i)=ho*den
        enddo 
        if(2*ns .lt. n-m) then
          dy=c(ns+1)
        else
          dy=d(ns)
          ns=ns-1
        endif
        y=y+dy
      enddo 
      RETURN
      END
     


C-----------------------------------------------------------
C---  SUBROUTINE  Exact_Time()
C---  Time the Exact_Value routine.

      SUBROUTINE  Exact_Time(rr,aa,bb,cc,yy)
      INTEGER NTest, i
      REAL dtime, tarray(2), time
      REAL rr,aa,bb,cc,yy

      NTest=1000000

      time=dtime(tarray)
      DO i=1,NTest
        CALL Exact_Value(rr,aa,bb,cc,yy)
      ENDDO
      time=dtime(tarray)
      
      WRITE(6,11) ' Exact_Value routine  :', time/NTest
  11  FORMAT(A, f12.8, ' seconds')

      RETURN
      END



C-----------------------------------------------------------
C---  SUBROUTINE  Inter_Time()
C---  Time the 2nd order Polynomial_Interpolation routine
C---  with 101 table entries.

      SUBROUTINE  Inter_Time(x1,y10,length,r)
      INTEGER NTest, i, m,j, length
      REAL dtime, tarray(2), time
      REAL r,x1(101),y10(101), pol_val(1)

      NTest=1000000
      
      m=101
   
      time=dtime(tarray)
      DO i=1,NTest
        CALL LOCATE(x1,m,r,j)                
        CALL POLINT(x1(j),y10(j),length,r,pol_val(1))
      ENDDO
      time=dtime(tarray)

      WRITE(6,11) ' Interpolation routine:', time/NTest
  11  FORMAT(A, f12.8, ' seconds',/)

      RETURN
      END

   



 


      

	


      
