spline.f90

FORTRANvisualfortran常用数值算法集及源码

SUBROUTINE spline(x,y,n,yp1,ypn,y2)
INTEGER n,NMAX
REAL yp1,ypn,x(n),y(n),y2(n)
PARAMETER (NMAX=500)
INTEGER i,k
REAL p,qn,sig,un,u(NMAX)
if (yp1>.99e30) then
  y2(1)=0.
  u(1)=0.
else
  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))
  p=sig*y2(i-1)+2.
  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
end do
if (ypn>.99e30) then
  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
  y2(k)=y2(k)*y2(k+1)+u(k)
end do
END SUBROUTINE spline