/*--------------------------------------------------------------------------
 * Spline3w.ox - code for cubic spline weights
 *
 *       code (C) Siem Jan Koopman 1996
 *		 20-09-97 (JAD): made more efficient as CubicSpline.
 *
 * This code computes a cubic spline weight matrix.
 *
 * Syntax:
 *
 * Spline3w(const vKn, const pmZ);
 *		vKn[k]    in:  vector with k knots
 *		pmZ		  in:  address of variable
 *				  out: T * k matrix with weights,
 *                     where T = vKn[k-1] - vKn[0] + 1
 * 
 *	Returns 1 if successful, 0 if failed.
 *--------------------------------------------------------------------------*/

#include <oxstd.h>

Spline3w(const vKn, const pmZ)
{
    decl k, i, j, t, ct, t1, t2;
    decl P, Q, S;
    decl r, s;
    decl a, b, u, zi, zi1;

    k = sizerc(vKn) - 1;
	if (k == 0)
		return 0;
		
    t1 = vKn[0];  t2 = vKn[k];
    ct = t2 - t1 + 1;

    P = 2.0 * unit(k+1);
    Q = zeros(k+1, k+1);

	for (i = 1, zi = vKn[1] - t1; i < k; i++)
    {
		zi1 = zi;  zi = vKn[i+1] - vKn[i]; 
		u = zi1 + zi;
		b = zi1 / u; 
        P[i][i-1] = b;
        P[i][i+1] = 1 - b;
        Q[i][i-1] = 6.0 / (zi1 * u);
        Q[i][i]  = -6.0 / (zi1 * zi);
        Q[i][i+1] = 6.0 / (zi * u);
    }
    S = invert(P) * Q;

    pmZ[0] = zeros(ct, k+1);

	for (i = j = 0, t = t1; t <= t2; ++t, ++j)
    {
        for (; t > vKn[i+1]; ) i++;

        if (t < vKn[i] || i > k)
        {
            print("\n Error: Spline is wrongly defined...\n");
            return 0;
        }
        a = vKn[i+1] - t; b = t - vKn[i];
		zi = a + b;
        r = s = zeros(1, k + 1); 
        r[i] = a / zi;
		r[i+1] = b / zi;
        s[i]   = a * ((a * a) - (zi * zi)) / (6.0 * zi);
        s[i+1] = b * ((b * b) - (zi * zi)) / (6.0 * zi);
        
        pmZ[0][j][] = r + (s * S);
    }
return 1;
}

