!!!!!!!!!!!!!  1D cubic spline interpolation !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine spline(x,y,n,yp1,ypn,y2)
implicit none
    INTEGER, PARAMETER :: iknd=SELECTED_INT_KIND(8)  ! nobody does this
    INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(14)  ! this is double precision
    INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6)  ! this is single precision   
!This subroutine is used to calculate the coefficients
!x, y, n are inputs, y2 is an output
!yp1 and ypn are inputs - the first derivative at each end; set to 1e30 for natural spline
integer :: n,NMAX
real(DP) :: x(n),y(n),y2(n)
real(DP) :: yp1,ypn,temp
parameter (NMAX = 500) !most number of points in x and y
integer :: i,k
real(DP) :: p,qn,sig,un,u(NMAX)

if(yp1 > .99e30)then
    y2(1) = 0           !natural spline boundaries
    u(1) = 0
else                    !specify first derivative
    y2(1) = -0.5
    u(1) = (3./(x(2) - x(1))) * ((y(2) - y(1))/(x(2) - x(1)) - yp1)
endif
do i = 2,n-1
    sig = (x(i) - x(i-1))/(x(i+1) - x(i-1)) !decomposition loop of tridiagonal algorithm. 
    p = sig * y2(i-1) + 2                   !y2 and u are used for temporary storage.
    y2(i) = (sig - 1.)/p
    u(i) = (6. * ((y(i+1) - y(i))/(x(i+1) - x(i)) - (y(i) - y(i-1)) &
            /(x(i) - x(i-1)))/(x(i+1) - x(i-1)) - sig * u(i-1))/p
enddo 
if (ypn > .99e30) then  !set upper boundary to natural spline or specify first deriv
    qn = 0.
    un = 0.
else
    qn = 0.5
    un = (3./(x(n) - x(n-1))) * (ypn - (y(n) - y(n-1))/(x(n) - x(n-1)))
endif
y2(n) = (un - qn * u(n-1))/(qn * y2(n-1) + 1.)
do k = n-1,1,-1 !backsubstitution of the triadiagonal algorithm
    y2(k) = y2(k) * y2(k + 1) + u(k)
enddo
return
end subroutine spline


subroutine splint(xa,ya,y2a,n,x,y)
implicit none
    INTEGER, PARAMETER :: iknd=SELECTED_INT_KIND(8)  ! nobody does this
    INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(14)  ! this is double precision
    INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6)  ! this is single precision   
!this is where the interpolation actually occurs
!give xa, ya, y2a from spline
!x is the value where you want to get y; y is the only output
integer :: n
real(DP) :: x,y,xa(n),y2a(n),ya(n)
integer :: k,khi,klo
real(DP) :: a,b,h
klo = 1
khi = n

!if(xa(klo) > x)then
!     k = klo 
!     go to 3
!endif
!if(xa(khi) < x)then
!     k = khi 
!     go to 3
!endif
!find the right place in the table using bisection
1 if ((khi - klo) > 1) then
    k = (khi + klo)/2
    if(xa(k) > x) then
        khi = k
    else
        klo = k
    endif
  goto 1
endif !khi and klo bracket x
h = xa(khi) - xa(klo)
if (h == 0.) pause 'bad xa input in splint' !xa's must be distinct
!3 if((k == klo).or.(k == khi))then
!    a = 0
!    b = 0
!else
    a = (xa(khi) - x)/h !cubic spline is now evaluated using linear interp
    b = (x - xa(klo))/h
!endif
y = a * ya(klo) + b * ya(khi) + &
    ((a**3 - a) * y2a(klo) + (b**3 - b) * y2a(khi)) * (h**2)/6.
return
end subroutine splint

!    integer :: i
!    integer, parameter :: n = 6    
!    real :: Tdp(n),Hcap(n),Hpow(n)
!    real :: yp1,ypn,y2(n),y
!       
!    
!    open (unit = 15,status = "old", &
!    file = "C:\Documents and Settings\apertzbo\My Documents\Research\standalone\Main_testing\Main_testing\interp stuff\heat.txt")
!    do i = 1,n
!        READ(15,*) Tdp(i),Hcap(i),Hpow(i)     
!    enddo
!    close(unit = 15)
!    
!!    do i = 1,n
!!        write(*,10) Tdp(i), Hcap(i), Hpow(i)
!!        10 format(f10.4,2x,f10.4,2x,f10.4)
!!    enddo
!    
!    call spline(Tdp,Hcap,n,1e30,1e30,y2)
!    
!    do i = 1,n
!        write(*,10) y2(i)
!        10 format(f10.4,2x,f10.4,2x,f10.4)
!    enddo
!
!    call splint(Tdp,Hcap,y2,n,17.6,y)
!    write(*,20) y
!    20 format(f10.4)


!!!!!!!!!!!!!  2D cubic spline interpolation !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine splie2(x1a,x2a,ya,m,n,y2a)
implicit none
    INTEGER, PARAMETER :: iknd=SELECTED_INT_KIND(8)  ! nobody does this
    INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(14)  ! this is double precision
    INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6)  ! this is single precision   
!use this subroutine to calculate the second derivatives
integer :: m,n,NN
real(DP) :: x1a(m), x2a(n), y2a(m,n), ya(m,n)
parameter (NN = 100)
integer :: j,k
real(DP) :: y2tmp(NN), ytmp(NN),temp

temp = 1.e30
do j = 1,m
    do k = 1,n
        ytmp(k) = ya(j,k)
    enddo
    call spline(x2a,ytmp,n,temp,temp,y2tmp)
    do k = 1,n
        y2a(j,k) = y2tmp(k)
    enddo
enddo

return
end subroutine splie2

subroutine splin2(x1a,x2a,ya,y2a,m,n,x1,x2,y)
implicit none
    INTEGER, PARAMETER :: iknd=SELECTED_INT_KIND(8)  ! nobody does this
    INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(14)  ! this is double precision
    INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6)  ! this is single precision   
!used to actually do the 2d interpolation
integer :: m,n,NN
real(DP) :: x1,x2,y,x1a(m),x2a(n),y2a(m,n),ya(m,n)
parameter (NN = 100)
integer :: j,k
real(DP) :: y2tmp(NN),ytmp(NN),yytmp(NN),temp

temp = 1.e30
do j = 1,m
    do k = 1,n
        ytmp(k) = ya(j,k)
        y2tmp(k) = y2a(j,k)
    enddo
    call splint(x2a,ytmp,y2tmp,n,x2,yytmp(j))
enddo
call spline(x1a,yytmp,m,temp,temp,y2tmp)
call splint(x1a,yytmp,y2tmp,m,x1,y)

return
end subroutine splin2

!!!!!!!!!!!!!!!! 2D linear interpolation !!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine linear(x1a,x2a,ya,m,n,x1,x2,y)
implicit none
    INTEGER, PARAMETER :: iknd=SELECTED_INT_KIND(8)  ! nobody does this
    INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(14)  ! this is double precision
    INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6)  ! this is single precision   
!this is where the interpolation actually occurs
!give xa, ya, y2a from spline
!x is the value where you want to get y; y is the only output
integer :: m,n,NN
real(DP) :: x1,x2,y,x1a(m),x2a(n),ya(m,n)
parameter (NN = 100)
integer :: k,khi,klo,j,jhi,jlo
real(DP) :: a,b,h
klo = 1
khi = m
if(x1a(klo) > x1)then
     k = klo 
     go to 3
endif
if(x1a(khi) < x1)then
     k = khi 
     go to 3
endif
!find the right place in the table using bisection
1 if ((khi - klo) > 1) then
    k = (khi + klo)/2
    if(x1a(k) > x1) then
        khi = k
    else
        klo = k
    endif
  goto 1
endif !khi and klo bracket x
3 jlo = 1
jhi = n
if(x2a(jlo) > x2)then
    j = jlo 
    goto 4
endif
if(x2a(jhi) < x2)then
    j = jhi 
    goto 4
endif

2 if ((jhi - jlo) > 1) then
    j = (jhi + jlo)/2
    if(x2a(j) > x2) then
        jhi = j
    else
        jlo = j
    endif
  goto 2
endif !khi and klo bracket x

4 if((k == klo).or.(k == khi))then
    a = 0
else
    k = k-1 
    a = (x1 - x1a(k))/(x1a(k+1) - x1a(k))
endif
if((j == jlo).or.(j == jhi))then
    b = 0
else
    j = j-1 
    b = (x2 - x2a(j))/(x2a(j+1) - x2a(j))
endif

y = (1-a)*(1-b)*ya(k,j) + a*(1-b)*ya(k+1,j) + a*b*ya(k+1,j+1) + (1-a)*b*ya(k,j+1)
return
end subroutine linear

!!!!!!!!!!!!!  1D polynomial interpolation !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine polint(xa,ya,n,x,y,dy)
implicit none
    INTEGER, PARAMETER :: iknd=SELECTED_INT_KIND(8)  ! nobody does this
    INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(14)  ! this is double precision
    INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6)  ! this is single precision   
integer :: n,NMAX
real(DP) :: dy,x,y,xa(n),ya(n)
parameter (NMAX = 10)
integer :: i,m,ns
real(DP) :: den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)
ns = 1
dif = abs(x - xa(1))

do i = 1,n
    dift = abs(x - xa(i))
    if (dift < dif) then
        ns = i
        dif = dift
    endif
    c(i) = ya(i)
    d(i) = ya(i)
enddo
y = ya(ns)
ns = ns - 1
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 == 0.) pause 'failure in point'
        den = w/den
        d(i) = hp*den
        c(i) = ho*den
    enddo
    if (2*ns < n-m) then
        dy = c(ns + 1)
    else
        dy = d(ns)
        ns = ns - 1
    endif
    y = y + dy
enddo

return
end subroutine polint

!!!!!!!!!!!!!!!!!!!! 2D polynomial interpolation !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine polin2(x1a,x2a,ya,m,n,x1,x2,y,dy)
implicit none
    INTEGER, PARAMETER :: iknd=SELECTED_INT_KIND(8)  ! nobody does this
    INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(14)  ! this is double precision
    INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6)  ! this is single precision   
integer :: m,n,NMAX,MMAX
real(DP) :: dy,x1,x2,y,x1a(m),x2a(n),ya(m,n)
parameter (NMAX = 20,MMAX = 20)
integer :: j,k
real(DP) :: ymtmp(MMAX), yntmp(NMAX)
do j = 1,m
    do k = 1,n
        yntmp(k) = ya(j,k)
    enddo
    call polint(x2a,yntmp,n,x2,ymtmp(j),dy)
enddo
call polint(x1a,ymtmp,m,x1,y,dy)

return
end subroutine polin2



!
!    integer, parameter :: m = 6, n = 5
!    real :: CTdp(m), CTwb(n), y2a(n*m), Ccap(m*n), Csens(m*n), Cpow(m*n)
!    real :: Ccap2(m,n), Csens2(m,n), Cpow2(m,n)
!    integer, parameter :: NN = 100
!    integer :: i,j,k,l
!    real :: y,x1,x2,y2,y3,dy
!
!    open (unit = 15,status = "old", &
!    file = "C:\Documents and Settings\apertzbo\My Documents\Research\standalone\Main_testing\Main_testing\interp stuff\cool.txt")
!    do j = 1,m*n+2
!        if (j == 1) then
!            READ(15,*) CTdp
!        elseif (j == 2) then
!            READ(15,*) CTwb
!        else
!            READ(15,*) Ccap(j-2),Csens(j-2),Cpow(j-2)  
!        endif
!    enddo  
!         
!    close(unit = 15)   
!    
!    l = 0
!    do k = 1,n
!        i = 1
!        do j = 1,m
!            Ccap2(j,k) = Ccap(i+l)
!            Csens2(j,k) = Csens(i+l)
!            Cpow2(j,k) = Cpow(i+l)
!            i = i + n
!        enddo
!        l = l+1
!    enddo
! 
!!    call WRITE_MATRIX(Ccap2)
!
!    call splie2(CTdp,CTwb,Ccap2,m,n,y2a)
!    call splin2(CTdp,CTwb,Ccap2,y2a,m,n,30.,16.944,y)   !input as x2,x1 not x1,x2 for when you want to interpolate
!    call linear(CTdp,CTwb,Ccap2,m,n,30.,16.944,y2)
!    call polin2(CTdp,CTwb,Ccap2,m,n,30.,16.944,y3,dy)
!    write(*,10) y,y2,y3
!    10 format (f10.4,2x,f10.4,2x,f10.4)