Berechnung eines interpolierenden kubischen Splines

program Splines(input, output);

    {Berechnung eines interpolierten kubischen natuerlichen Splines}
    
    uses
        stdfile;
        
    const
        MaxArray = 255;
    
    type
        MyArray = array[0..MaxArray] of double;
        
    var
        t, y:    MyArray;
        x:       double;
        n, i:    integer;
        h0, h1:  double;
        a, b:    double;
        u, v, w: MyArray;
        nd, hd:  MyArray;
        c:       MyArray;
        q:       double;
        
    procedure init;
        var
            fName: string;
            f:     text;
        begin {init}
            writeln('Berechnet ein interpoliertes kubisches natuerliches Spline');
            writeln;
            writeln('Welche Datei enthaelt die Punktepaare? (In splines.dat liegt ein Beispiel.)');
            readln(fName);
            if StdOpen(f, fName) then begin
                n := -1;
                repeat
                    n := n + 1;
                    if n > MaxArray then begin
                         writeln('Die Datei "', fName, '" hat zu viele ',
                             'Stuertzstellen (mehr als ', MaxArray, ').');
                         halt
                     end {if}
                 until not StdRead2(f, t[n], y[n]);
                 n := n - 1
            end else begin {if}
                writeln('Kann Datei "', fName, '" nicht oeffnen.');
                halt
            end; {if}
            close(f);
            writeln('Insgesamt ', n + 1, ' Stuetzstellen gelesen.');
            writeln;
            write('An welcher Stelle soll das Spline berechnet werden? ');
            readln(x);
            writeln
        end; {init}

    procedure Randbedingungen(var c: MyArray; u, v, w, y, t: MyArray);
        var
            h1, hn:             double;
            d1, dn:             double;
            b1, b2:             double;
            alpha, beta:        double;
            a11, a12, a21, a22: double;
            det:                double;
            i:                  integer;
        begin {Randbedingungen}
            h1 := t[1] - t[0];
            hn := t[n] - t[n - 1];
            d1 := (y[1] - y[0]) / h1;
            dn := (y[n] - y[n - 1]) / hn;
            b1 := 3 * d1 - u[1];
            b2 := 3 * dn - u[n - 1];
            a11 := 2 * v[1];
            a12 := w[1];
            a21 := v[n - 1];
            a22 := 2 + w[n - 1];
            det := a11 * a22 - a12 * a21;
            alpha := (b1 * a22 - a12 * b2) / det;
            beta := (a11 * b2 - b1 * a21) / det;
            c[0] := alpha;
            for i := 1 to n - 1 do begin
                c[i] := u[i] + alpha * v[i] + beta * w[i]
            end; {for}
            c[n] := beta
        end; {Randbedingungen}

    function Spline(x: double; t, y, c: MyArray): double;
        var
            k:              integer;
            h, delta:       double;
            a0, a1, a2, a3: double;
        begin {Spline}
            k := n + 1;
            repeat
                k := k - 1
            until x > t[k];
            k := k + 1;
            h := t[k] - t[k - 1];
            delta := (y[k] - y[k - 1]) / h;
            a0 := y[k];
            a1 := c[k];
            a2 := (c[k] - delta) / h;
            a3 := (c[k] - 2 * delta + c[k - 1]) / sqr(h);
            Spline := ((a3 * (x - t[k - 1]) + a2) * (x - t[k]) + a1) * (x- t[k]) + a0
        end; {Spline}
        
    procedure result;
        begin {result}
            writeln('Das Ergenis hat den Wert ', Spline(x, t, y, c):8:5)
        end; {result}

    begin {Splines}
        init;
        h0 := t[1] - t[0];
        a := (y[1] - y[0]) / sqr(h0);
        v[1] := -1 / h0;
        for i := 1 to n - 1 do begin
            h1 := t[i + 1] - t[i];
            b := (y[i + 1] - y[i]) / sqr(h1);
            u[i] := 3 * (a + b);
            hd[i] := 2 * (1 / h0 + 1 / h1);
            nd[i] := 1 / h1;
            a := b;
            h0 := h1
        end; {for}
        w[n - 1] := -1 / h1;
        {Vorwaertselimination}
        for i := 1 to n - 2 do begin
            q := -nd[i] / hd[i];
            hd[i + 1] := hd[i + 1] + nd[i] * q;
            u[i + 1] := u[i + 1] + q * u[i];
            v[i + 1] := q * v[i]
        end; {for}
        {Rueckwaertselimination}
        for i := n - 1 to 2 do begin
            q := hd[i];
            u[i] := u[i] / q;
            v[i] := v[i] / q;
            w[i] := w[i] / q;
            q := nd[i - 1];
            u[i - 1] := u[i - 1] - q * u[i];
            v[i - 1] := v[i - 1] - q * v[i];
            w[i - 1] := -q * w[i]
        end; {for}
        q := hd[1];
        u[1] := u[1] / q;
        v[1] := v[1] / q;
        w[1] := w[1] / q;
        Randbedingungen(c, u, v, w, y, t);
        result     
    end. {Splines}