C+
C NAME:
C	nrSpline
C PURPOSE:
C	Calculate second derivatives in X locations of interpolating spline
C	for Y fnc-values (the 2nd derivative array serves as input to
C	the interpolation subroutine nrSplInt)
C CATEGORY:
C	Math: interpolation and extrapolation
C CALLING SEQUENCE:
	subroutine nrSpline(N,X,Y,YP1,YPN,Y2)
C INPUTS:
C	N	integer		# points and fnc-values
C	X(N)	real		x-values; MUST be sorted
C	Y(N)	real		function values
C	YP1	real		derivative at point 1
C				if YP1>=1.e30 then 2nd derivative is zero in point 1
C	YPN	real		derivative at point N
C				if YP1>=1.e30 then 2nd derivative is zero in point N
C OUTPUTS:
C	Y2(N)	real		second derivatives at locations X
C CALLS:
C	Say
C RESTRICTIONS:
C	The X array must be sorted.
C	NO CHECK IS MADE TO ENSURE THAT THIS IS THE CASE
C PROCEDURE:
C	See Numerical Recipes, par. 3.4, p. 88
C MODIFICATION HISTORY:
C	JUN-1993, Paul Hick (UCSD)
C-
	    integer	N
	    real	X (N)
	    real	Y (N)
	    real	Y2(N)

	parameter      (NMAX = 1000)
	real		U(NMAX)

	if (N .gt. NMAX) call Say('nrSpline','E','Stop','Too small parameter NMAX')

	if (YP1 .gt. 0.99e30) then
	    Y2(1) = 0.0
	    U(1)  = 0.0
	else
	    Y2(1) = -0.5
	    U(1)  = (3.0/(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)
	end if

	do I=2,N-1
	    SIG = (X(I)-X(I-1))/(X(I+1)-X(I-1))
	    P = SIG*Y2(I-1)+2.0
	    Y2(I) = (SIG-1.0)/P
	    U(I) = (6.0*((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 .gt. 0.99e30) then
	    QN = 0.0
	    UN = 0.0
	else
	    QN = 0.5
	    UN = (3.0/(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))
	end if

	Y2(N) = (UN-QN*U(N-1))/(QN*Y2(N-1)+1.0)
	do K=N-1,1,-1
	    Y2(K) = Y2(K)*Y2(K+1)+U(K)
	end do

	return
	end