subroutine spline(x,y,n,yp1,ypn,y2) * * Given arrays X & Y (Y=f(X)) of length n containing a tabulated * function, and given the values yp1 and ypn of the first * derivatives at the endpoints, this routine returns an array Y2 * of length n containing the second derivatives ofthe interpolating * function at the tabulated points of X. * If yp1 and/or ypn are 10^{30} or larger, the routine sets the * BC's for a natural spline, with zero second derivatives at that * boundary. * Routine is from Numerical Recipies * implicit double precision(a-h,o-z) parameter (nmax=10000) dimension x(n),y(n),y2(n),u(nmax) if(yp1.gt.0.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 11 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)) 1 /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p 11 continue if(ypn.gt.0.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 12 k=n-1,1,-1 y2(k)=y2(k)*y2(k+1)+u(k) 12 continue return end ************************************************************************