!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Rational Interpolation  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!  A variation of polynomial interpolation is used.  Rational  !!
!!  interpolation approximates a given function by considering  !! 
!!  the ratio of two polynomials.  Unlike generic polynomial    !!
!!  interpolation, rational interpolation is much more stable   !!
!!  when considering equidistant step sizes.  A side by side    !!
!!  comparison is made with generic polynomial interpolation.   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!  compiled with gfortran inverse.f95 rat_inter.f95  !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  ABW 2013 - 2015  !
!!!!!!!!!!!!!!!!!!!!!
program rational
  use inv_mod
  implicit none
  
  integer,parameter :: n = 10  !  Number of points to interpolate
  real(kind=16), parameter :: pi = 3.14159265358979323846264338327950q0  !  Pi to 32 digits
  real(kind=16), parameter :: a = -2.0q0*pi, b = 2.0q0*pi  !  Limits of interpolation [a,b]
  real(kind=16),dimension(n+1) :: xx,yy,bR,u,v
  real(kind=16),dimension(n+1,n+1) :: AA, AAinv
  real(kind=16) :: h, x0, x1, x, errmax1,errmax2, h0, err1,err2, xj
  integer :: i,j, xit
  
  !  Step size
  h = (b - a)/n
  
  !  Collecting x and y data
  do i = 1,n+1
     xx(i) = a + (i - 1.0q0)*h
     yy(i) = f(xx(i))
  enddo
  
  x0 = xx(1)
  x1 = xx(n+1)
  
  do i = 1,n+1
     AA(i,1) = 1.0q0
     do j = 2,n+1
        AA(i,j) = (AA(i,j-1))*xx(i)
     enddo
     bR(i) = (Q(xx(i)))*(yy(i))
  enddo
  
  !  Inverse by Gauss-Jordan Elimin. with quad precision
  call matinv(AA,AAinv,n+1)
  
  u = matmul(AAinv,yy)  !  to be used with P(x)
  v = matmul(AAinv,bR)  !  to be used with R(x)
  
  
  !!!!!!!!!!!!!!!!!!!!!!!  
  !  Data for Plotting  !
  !!!!!!!!!!!!!!!!!!!!!!!
   
  ! The function itself
  open (unit=10,file="func_data.dat",action="write",access="append")  
  do xit = -628,628,1
     x = real(xit)/100.0q0
     write(10,*) x, f(x)
  enddo
  close (10)
  
  ! The Rational interpolation of the function
  open (unit=11,file="n20_inter.dat",action="write",access="append")
  do xit = -628,628,1
     x = real(xit)/100.0q0
     write(11,*) x, R(x)
  enddo
  close (11)
  
  !!!!!!!!!!!!!!!!!!!!!!!  
  !  Done taking data   !
  !!!!!!!!!!!!!!!!!!!!!!!
  
  !  Calculating error in the interpolation
  errmax1 = 0.0q0
  errmax2 = 0.0q0
  do i = 1, n
     h0 = (xx(i+1) - xx(i))/5
     do j = 1, 5
        xj = xx(i) + real(j)*h0
        err1 = abs(f(xj) - R(xj))  !  Error in rational interpolation
        err2 = abs(f(xj) - P(xj))  !  Error in generic polynomial interpolation
        if (err1 > errmax1) then
           errmax1 = err1
        endif
        if (err2 > errmax2) then
	   errmax2 = err2
        endif
     enddo
  enddo
  
  write(*,120) 'Maximum Error in Rational Interpolation:',errmax1,' with ',n,' points.'
  write(*,*)
  write(*,120) 'Maximum Error in Polynomial Interpolation:',errmax2,' with ',n,' points.'
120 format (A,F20.16,A,I2,A)
  
contains
  
  !  Denominator of rational function
  function Q(x)
    implicit none    
    real(kind=16) :: Q, x    
    
    Q = (x - x0)**n + (x1 - x)**n
    return    
  end function Q
  
  !  Complete rational function
  function R(x)
    implicit none    
    real(kind=16) :: Rsum, R, x
    integer :: k  
    
    Rsum = 0.0q0
    do k = 0,n
       Rsum = Rsum + (x**k)*(v(k+1))
    enddo
    R = Rsum/(Q(x))    
    return   
  end function R
  
  !  Polynomial function
  function P(x)
    implicit none
    real(kind=16) :: Psum,P,x
    integer :: k
    
    Psum = 0.0q0
    do k = 0,n
       Psum = Psum + (x**k)*(u(k+1))
    enddo
    P = Psum
    return
  end function P
  
  !  Function to interpolate
  function f(x)
    implicit none    
    real(kind=16) :: f, x
    
!    f = x*x*sin(x)  !  For low n, rational and polynomial interpolations are about the same
    f = 1.0q0 + sin(x)*sin(x)  !  For low n, rational interpolation is much better
    return    
  end function f
  
end program rational
