IMPLEMENTATION MODULE CubicSplines;
(*
This Module implements the cubic splines. See "Numerical Methods for
Scientists and Engineers" by R. W. Hamming, page 349.
*)

FROM InOut IMPORT WriteString, WriteReal, WriteLn, WriteInt, WriteCard;
FROM Storage IMPORT ALLOCATE, DEALLOCATE;

CONST
    MAXPTS = 100;  (* maximum data points that can be splined *)

TYPE
    Spline = POINTER TO SplineType;
    SplineType = RECORD
                     x, y, ydd : ARRAY [0..MAXPTS] OF REAL;
                     n         : CARDINAL; (* last index *)
                 END;
             
PROCEDURE MakeSpline(
        x : ARRAY OF REAL (* in *);
        y : ARRAY OF REAL (* in *);
        un: CARDINAL      (* in *);
    VAR S : Spline        (* out *) ) : BOOLEAN;
    (*
    This procedure takes the user x,y data and creates a cubic spline S.
    The x values must be increasing either in the direction x[1]..x[un] or
    x[un]..x[1] but must never be such that two x values overlap or cross.
    *)
    
    PROCEDURE FindOrder(
            x : ARRAY OF REAL (* in *);
            n : CARDINAL      (* in *);
        VAR from, to, by : INTEGER (* out *) ) : BOOLEAN;
        (*
        This function scans the x array to find the direction in which
        the values are increasing, i.e, 0..n or n..0. If such an ordering
        is found, then "from", "to" and "by" are set to reflect this
        order, and the function returns TRUE; FALSE otherwise.
        *)
        VAR
            i : INTEGER;
            done : BOOLEAN;
        BEGIN
            done := TRUE;
            from := 0; to := n; by := 1;
            i := 0;
            WHILE ( done ) AND ( i < to ) DO
                IF x[i+1] <= x[i] THEN
                    done := FALSE
                ELSE
                    i := i + by
                END;
            END;
            IF done THEN
                RETURN TRUE
            END;
            (*
            Try backwards.
            *)
            done := TRUE;
            i := n;
            from := n; to := 0; by := (-1);
            WHILE ( done ) AND ( i > to ) DO
                IF x[i] >= x[i-1] THEN
                    done := FALSE
                ELSE
                    i := i + by
                END;
            END;
            IF done THEN
                RETURN TRUE;
            END;
            RETURN FALSE;
        END FindOrder;

    PROCEDURE MakeAxb(
        x,y : ARRAY OF REAL (* in *);
        n   : CARDINAL      (* in *);
        VAR sub,diag,sup,b : ARRAY OF REAL (* out *) );
        (*
        Local procedure. Sets up the tridiagonal system Ax = b
        *)
        VAR
            h1, h2 : REAL;
            i      : CARDINAL;
        BEGIN
            FOR i := 1 TO n-1 DO
                h1 := x[i] - x[i-1];
                h2 := x[i+1] - x[i];
                sub[i] := h1;
                diag[i] := 2.0*( h1 + h2 );
                sup[i] := h2;
                
                b[i] := ( y[i+1]-y[i] )/h2 - ( y[i]-y[i-1] )/h1;
                b[i] := 6.0*b[i];
            END;
            sub[0] := 0.0; diag[0] := 1.0; sup[0] := 0.0;
            sub[n] := 0.0; diag[n] := 1.0; sup[n] := 0.0;
        END MakeAxb;
    
    PROCEDURE SolveTridiagonal(
        sub, diag, sup : ARRAY OF REAL (* in *);
        n              : CARDINAL      (* in *);
        VAR b          : ARRAY OF REAL (* in/out *) ) : BOOLEAN;
        (*
        Solves a tridiagonal system. returns false if trouble occurs.
        Solution is returned in "b" which contains the rhs upon entry.
        *)
        VAR
            k : CARDINAL;
            m : REAL;
        BEGIN
            FOR k := 1 TO n DO
                IF diag[k] = 0.0 THEN
                    WriteString("MakeSpline -- trouble, diagonal entry " );
                    WriteString("is zero at i " );
                    WriteCard( k, 4 );
                    WriteLn;
                    RETURN FALSE;
                END;
                m := sub[k]/diag[k-1];
                diag[k] := diag[k] - m*sup[k-1];
                b[k]    := b[k]    - m*b[k-1];
            END; (* FOR *)
            IF diag[n] = 0.0 THEN
                WriteString("MakeSpline -- trouble, diagonal entry " );
                WriteString("is zero at i " );
                WriteCard( n, 4 );
                WriteLn;
                RETURN FALSE;
            END;
        
            b[n] := b[n]/diag[n];
            FOR k := n-1 TO 0 BY -1 DO
                b[k] := ( b[k] - sup[k]*b[k+1] )/ diag[k];
            END;
            RETURN TRUE;
        END SolveTridiagonal;
    
    VAR
        i : INTEGER;
        j : CARDINAL;
        from, to, by : INTEGER;
        sub, diag, sup, b : ARRAY [0..MAXPTS] OF REAL;
    BEGIN (* MakeSpline *)
        IF un <= 2 THEN
            WriteString("MakeSpline -- No use making spline out of two " );
            WriteString("or less points" );
            WriteLn;
            RETURN FALSE;
        END;
        (*
        Make sure un-1 <= MAXPTS.
        *)
        IF un-1 > MAXPTS THEN
            WriteString("MakeSpline -- Can handle " );
            WriteCard( MAXPTS,4 );
            WriteString(" data points only ");
            WriteLn;
            RETURN FALSE;
        END;
    
        IF NOT FindOrder( x, un-1, from,to,by ) THEN
            WriteString("MakeSpline -- x values are not increasing in ");
            WriteString("either direction"); WriteLn;
            RETURN FALSE;
        END;
        NEW( S );
        S^.n := un-1;  (* last index *);
        
        i := from;
        j := 0;
        LOOP
            S^.x[j] := x[i];
            S^.y[j] := y[i];
            S^.ydd[j] := 0.0;
            IF i = to THEN
                EXIT
            END;
            INC( j );
            i := i + by;
        END;
    
        MakeAxb( S^.x, S^.y, S^.n, sub, diag, sup, b );
        IF NOT SolveTridiagonal( sub,diag,sup,S^.n, b ) THEN
            WriteString("MakeSpline -- Trouble in making spline" );
            WriteLn;
            DISPOSE( S );
            RETURN FALSE;
        END;
    
        FOR i := 1 TO INTEGER( S^.n-1 ) DO
            S^.ydd[i] := b[i];
        END;
        S^.ydd[0] := 0.0;
        S^.ydd[S^.n] := 0.0;
        RETURN TRUE;
    END MakeSpline;

PROCEDURE FindIndex(
    x : ARRAY OF REAL (* in *);
    ux: REAL          (* in *);
    n : CARDINAL      (* in *)  ): INTEGER;
    (*
    Internal proc. Returns index i such that x[i] <= ux <= x[i+1]. 
    i would range 0-n. If no bound found, a -1 is returned.
    *)
    VAR
        i : CARDINAL;
    BEGIN
        IF ux = x[0] THEN
            RETURN 0
        END;
        IF ux = x[n] THEN
            RETURN n
        END;
        FOR i := 0 TO n-1 DO
            IF ( ux >= x[i] ) AND ( ux <= x[i+1] ) THEN
                RETURN i;
            END;
        END;
        RETURN (-1);
    END FindIndex;

PROCEDURE EvalSpline(
          x : REAL   (* in *);
          S : Spline (* in *);
    VAR val : REAL   (* out *) ) : BOOLEAN;
    (*
    Evaluate spline S at x and return value in "val". If trouble, e.g,
    x is out of range, RETURN FALSE.
    *)
    VAR
        index : INTEGER;
        i     : CARDINAL;
        h, r1, r2 : REAL;
    BEGIN
        IF S = NIL THEN
            WriteString("EvalSpline -- S is not a spline" );
            WriteLn;
            RETURN FALSE;
        END;
        index := FindIndex( S^.x, x, S^.n );
        IF index < 0 THEN
            WriteString("EvalSpline -- " );
            WriteReal( x, 10 );
            WriteString(" is out of x-range ");
            WriteLn;
            RETURN FALSE;
        END;
        i := CARDINAL( index );
        IF x = S^.x[0] THEN
            val := S^.y[0];
            RETURN TRUE;
        END;
        IF x = S^.x[S^.n] THEN
            val := S^.y[S^.n];
            RETURN TRUE;
        END;
        (* interior *)
        h := S^.x[i+1] - S^.x[i];
        r1 := ( x - S^.x[i] )/h;
        r2 := ( S^.x[i+1] - x )/h;
        
        val := S^.y[i]*r2 + S^.y[i+1]*r1;
        h := h*h/6.0;
        val := val - h*S^.ydd[i]*( r2 - r2*r2*r2 );
        val := val - h*S^.ydd[i+1]*( r1 - r1*r1*r1 );
        
        RETURN TRUE;
    END EvalSpline;

PROCEDURE DerivSpline(
            x : REAL     (* in *);
            S : Spline   (* in *);
    VAR Deriv : REAL     (* out *) ): BOOLEAN;
    (*
    Evaluates the derivative of the spline "S" at x. The function value
    returned is TRUE if ok, otherwise if x is out of range or if other
    problems occur, FALSE is returned.
    *)
    VAR
        index : INTEGER;
        i     : CARDINAL;
        h, r1, r2 : REAL;
        n     : CARDINAL;
    BEGIN
        IF S = NIL THEN
            WriteString("DerivSpline -- S is not a spline" );
            WriteLn;
            RETURN FALSE;
        END;

        index := FindIndex( S^.x, x, S^.n );
        IF index < 0 THEN
            WriteString("DerivSpline -- " );
            WriteReal( x, 10 );
            WriteString(" is out of x-range ");
            WriteLn;
            RETURN FALSE;
        END;
    
        i := CARDINAL( index );
        IF i = S^.n THEN
            n := S^.n;      (* right boundary *)
            h := S^.x[n] - S^.x[n-1];
            Deriv := ( S^.y[n]-S^.y[n-1] )/h;
            Deriv := Deriv + h*( S^.ydd[n-1] + 2.0*S^.ydd[n] )/6.0;
            RETURN TRUE;
        END;
    
        h := S^.x[i+1] - S^.x[i];
        r1 := ( x - S^.x[i] )/h;
        r2 := ( S^.x[i+1] - x )/h;
        Deriv := ( S^.y[i+1] - S^.y[i] )/h;
        Deriv := Deriv + h*S^.ydd[i]*( 1.0 - 3.0*r2*r2 )/6.0;
        Deriv := Deriv - h*S^.ydd[i+1]*( 1.0 - 3.0*r1*r1 )/6.0;
        
        RETURN TRUE;
    END DerivSpline;

END CubicSplines.