# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module spline_mod
  implicit none
  real, pointer,private :: deriv(:),x_in(:),y_in(:)
  contains
  subroutine spline_close()
    deallocate(deriv)
  end subroutine
  logical function spline_init(xin,yin,yp1,ypn)
    real,pointer ::xin(:),yin(:)
    real, optional :: yp1,ypn
    integer :: i,k,n
    real :: p,q,qn,sig,un,u(size(xin))
    n=size(xin)
    allocate(deriv(n))
    x_in=>xin
    y_in=>yin
    if (present(yp1)) then
      deriv(1)=-.5
      u(1)=(3./(xin(2)-xin(1)))*((yin(2)-yin(1))/(xin(2)-xin(1))-yp1)
    else
      deriv(1)=-1
      u(1)=0
    end if
    write(0,*) "WJGDJKGD ",n
    do i=2,n-1
      sig=(xin(i)-xin(i-1))/(xin(i+1)-xin(i-1))
      p=sig*deriv(i-1)+2
      deriv(i)=(sig-1.)/p
      u(i)=(6.*((yin(i+1)-yin(i))/(xin(i+1)-xin(i))-(yin(i)-yin(i-1))/&
        & (xin(i)-xin(i-1)))/(xin(i+1)-xin(i-1))-sig*u(i-1))/p
    end do
    if (present(ypn)) then
      qn=.5
      un=(3./(xin(n)-xin(n-1)))*(ypn-(yin(n)-yin(n-1))/(xin(n)-xin(n-1&
        &)))
    else
      qn=0
      un=0.
    end if
    deriv(n)=(un-qn*u(n-1))/(qn*deriv(n-1)+1.)
    do k=n-1,1,-1
      deriv(k)=deriv(k)*deriv(k+1)+u(k)
    end do
    spline_init=.true.
  end function
  real function spline_interp(x)
    integer :: k,khi,klo,n
    real :: a,b,h,x
    logical :: found
    n=size(x_in)
    klo=1
    khi=n
    found=.false.
    do while (.not. found)
      if (khi-klo >1) then
        k=(khi+klo)/2.
        if (x_in(k)>x) then
          khi=k
        else
          klo=k
        end if
      else
        found=.true.
      end if
    end do
    h=x_in(khi)-x_in(klo)
    a=(x_in(khi)-x)/h
    b=(x-x_in(klo))/h
    spline_interp=a*y_in(klo)+b*y_in(khi)+((a**3-a)*deriv(klo)+(b**3-b&
      &)*deriv(khi))*(h**2)/6.
  end function
end module
