	subroutine  cubspl( tau   , c	  , n	  , ibcbeg, ibcend)
c
c	n		number of data points, ge. 2
c	(tau(i),c(1,i), i = 1,...,n)
c			abscissae and ordinates of the data points.
c			tau is assumed to be strictly increasing.
c	ibcbeg, ibcend	boundary conditions indicators, and
c	c(2,1), c(2,n)	boundary condition information. specifically,
c		ibcbeg	=  0
c			means no boundary condition at tau(1) is given.
c			in this case, the not-a-knot condition is used,
c			i.e. the jump in the third derivative across
c			tau(2) is forced to zero, thus the first and
c			the second cubic polynomial pieces are made
c			to coincide.
c		ibcbeg	=  1
c			means that the slope at tau(1) is made
c			to equal c(2,1), supplied by input.
c		ibcbeg	=  2
c			means that the second derivative at tau(1)
c			is made equal to c(2,1), supplied by input.
c		ibcend	=  0 , 1 , 2
c			has analogus meaning concerning the boundary
c			condition at tau(n), with the addional
c			information taken from c(2,n).
c
c
c	c(j,i), j = 1,...,4; i = 1,...,l := n-1
c			the polynomial coefficients of the cubic
c			interpolating spline with interior knots
c			(or joints) tau(2),...,tau(n-1).
c			precisely, in the interval
c			(tau(i),tau(i+1)), the spline f is given by
c
c			f(x)  =  c(1,i) + h*(c(2,i) + h*(c(3,i)
c				 + h*c(4,i)/3.)/2.)
c			where  h  := x - tau(i). the function
c			program  *ppvalu*  may be used to
c			evaluate  f  or its derivatives from  tau,
c			c , l = n-1 , and  k = 4.
c
c
        implicit double precision (a-h,o-z)
	integer  ibcbeg, ibcend, n     , i     , j     , l	,
     *		 m
c
	real*8	 c     ( 4     , n    ), tau   ( n    ), divdf1 ,
     *		 divdf3, dtau  , g
c
c
c	a tridiagonal linear system for the unknown slopes
c	s(i) of  f  at	tau(i), i = 1,...,n  is generated and
c	then solved by gauss elimination, with	s(i)  ending
c	up in  c(2,i), all i.
c	c(3,.) and c(4,.)  are used initially for temporary
c	storage.
c
	l  =  n - 1
c
c	compute first differences of tau sequence and store in c(3,.).
c	also, compute first divided difference of data and
c	store in c(4,.).
c
	do  10	m = 2,n
		c(3,m)	=  tau(m) - tau(m-1)
		c(4,m)	=  (c(1,m) - c(1,m-1))/c(3,m)
10	continue
c
c	constuct first equation from the boundary condition,
c	of the form
c		c(4,1)*s(1) + c(3,1)*s(2)  =  c(2,1)
c
	if  (ibcbeg .eq. 2)  go to  16
c
		if  (ibcbeg .eq. 1)  go to 15
			if  (n .gt. 2)	go to  12
c
c				no condition at left and  n = 2.
c
				c(4,1)	=  1.d0
				c(3,1)	=  1. d0
				c(2,1)	=  2.d0*c(4,2)
			go to  25
12			continue
c
c			not-a-knot condition at left and n.gt. 2
c
			c(4,1)	=  c(3,3)
			c(3,1)	=  c(3,2) + c(3,3)
			c(2,1)	=  ((c(3,2) + 2.d0*c(3,1))*c(4,2)*
     *				    c(3,3) + c(3,2)**2*c(4,3))/c(3,1)
		go to  19
15		continue
c
c		slope prescribed at left end.
c
		c(4,1)	=  1.d0
		c(3,1)	=  0.d0
		go to  18
16	continue
c
c	second derivative prescribed at left end.
c
	c(4,1)	=  2.d0
	c(3,1)	=  1.d0
	c(2,1)	=  3.d0*c(4,2) - c(3,2)/2.d0*c(2,1)
c
18	continue
	if  (n .eq. 2)	go to  25
c
c		if there are interior knots, generate the corresp. equations
c		and carry out the forward pass of gauss elimination, after
c		which the m-th equation reads
c			c(4,m)*s(m) + c(3,m)*s(m+1)  =	c(2,m).
c
19		continue
		do  20	m = 2,l
			g  =  - c(3,m+1)/c(4,m-1)
			c(2,m)	=  g*c(2,m-1) + 3.d0*(c(3,m)*c(4,m+1) +
     *						    c(3,m+1)*c(4,m))
			c(4,m)	=  g*c(3,m-1) + 2.d0*(c(3,m) + c(3,m+1))
20		continue
c
c		construct last equation from the second boundary condition,
c		of the form
c			(-g*c(4,n-1))*s(n-1) + c(4,n)*s(n)  =  c(2,n)
c
c		if the slope is prescribed at right end, one can go directly
c		to back substitution, since c array happens to be set up
c		just right for it at this point.
c
		if ( ibcend .eq. 1)  go to  30
			if  (ibcend .eq.2)  go to  24
c
c
				if  ((n .eq. 3).and.(ibcbeg.eq.0))
     *				    go	to  22
c
c					knot-a-knot and n.ge.3, and either
c					n.gt.3 or also knot-a-knot at
c					left end point.
c
					g  =  c(3,n-1) + c(3,n)
					c(2,n)	=  ((c(3,n) + 2.d0*g)*
     *						    c(4,n)*c(3,n-1) +
     *						    c(3,n)**2*
     *						    (c(1,n-1)-c(1,n-2))
     *						    /c(3,n-1))/g
					g  =  -g/c(4,n-1)
					c(4,n)	=  c(3,n-1)
				go  to	29
22				continue
c
c				either (n=3 and not-a-knot also at left)
c				or (n=2 and not-a-knot at left end point).
c
				c(2,n)	=  2.d0*c(4,n)
				c(4,n)	=  1.d0
				go  to	28
24			continue
			c(2,n)	=  3.d0*c(4,n) + c(3,n)/2.d0*c(2,n)
			c(4,n)	=  2.d0
			go  to	28
25			continue
			if  (ibcend .eq. 1)  go to  30
c
				if  (ibcend .eq. 0)  go to  26
c
					go  to	24
c
26				continue
 				if  (ibcbeg .gt. 0)  go to  22
c
c					not-a-knot at right endpoint and
c					at left endpoint and n = 2.
c
					c(2,n)	=  c(4,n)
				go  to	30
28			continue
			g  =  -1.d0/c(4,n-1)
29		continue
c
c		complete forward pass of gauss elimination.
c
		c(4,n)	=  g*c(3,n-1) + c(4,n)
		c(2,n)	=  (g*c(2,n-1) + c(2,n))/c(4,n)
c
30	continue
c
c	carry out back substitution
c
	do  40	j = l,1,-1
		c(2,j)	=  (c(2,j) - c(3,j)*c(2,j+1))/c(4,j)
40	continue
c
c	generate cubic coefficients in each interval, i.e. the deriv.s
c	at its left endpoint, from value and slope at its endpoints.
c
 	do  50	i = 2,n
		dtau  =  c(3,i)
		divdf1	=  (c(1,i) - c(1,i-1))/dtau
		divdf3	=  c(2,i-1) + c(2,i) - 2.d0*divdf1
		c(3,i-1)  =  2.d0*(divdf1 - c(2,i-1) - divdf3)/dtau
		c(4,i-1)  =  (divdf3/dtau)*(6.d0/dtau)
50	continue
c
	return
	end
 
