Integration einer diskreten Funktion

program Integration(input, output);

    {Integration einer diskreten Funktion}
    
    uses
        stdfile;
    
    const
        MaxArray = 1024;
    
    var
        x, fx:      array [0..MaxArray] of double;
        n:          integer;
        res:        double;
        t:          double;
        h, hfuture: double;
        k, j:       integer;
        
    procedure init;
        var
            fName: string;
            f:     text;
        begin {init}
            writeln('Integration einer diskreten Funktion');
            writeln;
            writeln('In welcher Datei liegen die Punktepaare? (In integr.dat liegt ein Beispiel.)');
            readln(fName);
            if StdOpen(f, fName) then begin
                n := -1;
                repeat
                    n := n + 1;
                    if n > MaxArray then begin
                        writeln('Zu viele Paare (mehr als ', MaxArray,
                            ') in der Datei.');
                        halt
                    end {if}
                until not StdRead2(f, x[n], fx[n]);
                n := n - 1;
                writeln('Insgesamt ', n + 1, ' Werte gelesen.')
            end else begin {if}
                writeln('Die Datei "', fName, '" kann nicht geoeffnet werden.');
                halt
            end; {if}
            close (f)
        end; {init}

    procedure result;
        begin {result}
            writeln('Die Loesung des Integrals lautet ', res:8:4)
        end; {result}

    begin {Integration}
        init;
        h := x[1] - x[0];
        k := 1;
        res := 0;
        for j := 1 to n - 1 do begin
            hfuture := x[j + 1] - x[j];
            if h = hfuture then begin
                if k = 3 then begin
                    res := res + 2 * h * (fx[j - 1] + 4 * fx[j - 2] + fx[j - 3]) / 6;
                    k := k - 1
                end else begin {if}
                    k := k + 1
                end {if}
            end else begin {if}
                if k = 1 then begin
                    res := res + h * (fx[j] + fx[j - 1]) / 2
                end else begin {if}
                    if k = 2 then begin
                        res := res + 2 * h * (fx[j] + 4 * fx[j - 1] + fx[j - 2]) / 6
                    end else begin {if}
                        t := fx[j] + 3 * (fx[j - 1] + fx[j - 2]) + fx[j - 3];
                        res := res + 3 * h * t / 8
                    end; {if}
                    k := 1
                end {if}
            end; {if}
            h := hfuture
        end; {for}
        result     
    end. {Integration}