--- a
+++ b/Integral.mod.m2pp
@@ -0,0 +1,538 @@
+IMPLEMENTATION MODULE Integral;
+
+  (*========================================================================*)
+  (* HINWEIS : Bitte nur die Datei Integral.mod.m2pp editieren              *)
+  (*========================================================================*)
+  (* Es sind 2 Versionen enthalten die mit                                  *)
+  (*                                                                        *)
+  (*   m2pp -D __{Parameter}__ < Integral.mod.m2pp > Integral.mod           *)
+  (*                                                                        *)
+  (* mit Parameter = {XDS|GM2} erzeugt werden koennen.                      *)
+  (*                                                                        *)
+  (*   XDS    : Der Export von Objekten aus lokalen Modulen wird mit        *)
+  (*            "IMPORT" angezeigt - das scheint nicht ganz "regelkonform"  *)
+  (*   GM2    : Normales "EXPORT" wird genutzt (ELSE-Fall).                 *)
+  (*                                                                        *)
+  (* ansonsten gibt es keine Aenderungen am Quellcode                       *)
+  (*                                                                        *)
+  (* Fuer die Routine InitGausInt kann Debug-Code mit der Option            *)
+  (* -D __DEBUG__ eingefuegt werden                                         *)
+  (*------------------------------------------------------------------------*)
+  (* Module providing procedures for numerical integration                  *)
+  (*------------------------------------------------------------------------*)
+  (* Letzt Bearbeitung:                                                     *)
+  (*                                                                        *)
+  (* 17.07.93, MRi: Durchsicht                                              *)
+  (* 01.03.18, MRi: Kosmetische Korrekturen, TriQR durch TriQL ersetzt      *)
+  (* 02.03.18, MRi: Korrektur in NumInt, Simson & SimCum eingefuegt         *)
+  (*------------------------------------------------------------------------*)
+  (* Testprogram in TstIntegral                                             *)
+  (*------------------------------------------------------------------------*)
+  (* Implementation : Michael Riedl                                         *)
+  (* Licence        : GNU Lesser General Public License (LGPL)              *)
+  (*------------------------------------------------------------------------*)
+
+  (* $Id: Integral.mod.m2pp,v 1.5 2018/03/02 10:48:20 mriedl Exp mriedl $ *)
+
+FROM SYSTEM    IMPORT TSIZE;
+FROM Storage   IMPORT ALLOCATE,DEALLOCATE;
+FROM LongMath  IMPORT sqrt;
+FROM LMathLib  IMPORT MachEps,Funktion1;
+FROM Deklera   IMPORT VEKTOR,CARDVEKTOR,MaxDim;
+FROM Errors    IMPORT Fehler,Fehlerflag,ErrOut,FatalError;
+FROM MatLib    IMPORT QuickSort;
+FROM EigenLib1 IMPORT ListenZeiger,TriQL,BerechneEV;
+<* IF (__DEBUG__) THEN *>
+IMPORT TIO;
+<* END *>
+
+MODULE GAUSSINT;
+
+IMPORT VEKTOR,CARDVEKTOR,MaxDim;
+IMPORT QuickSort;
+IMPORT ListenZeiger,TriQL,BerechneEV;
+IMPORT sqrt;
+IMPORT MachEps,Funktion1;
+<* IF (__DEBUG__) THEN *>
+IMPORT TIO;
+<* END *>
+
+<* IF (__XDS__) THEN *>
+IMPORT InitGaussInt,GaussInt; (* EXPORT *)
+<* ELSE *>
+EXPORT InitGaussInt,GaussInt;
+<* END *>
+
+(*------------------ Lokale Objekte ------------------------------*)
+
+CONST     genau   = 10.0*MachEps;
+
+VAR       A,X     : VEKTOR;
+          Anzahl  : CARDINAL; (* Anzahl der Stuetzstellen *)
+
+(*----------------------------------------------------------------*)
+
+PROCEDURE InitGaussInt(NSt : CARDINAL);
+
+          VAR  HD,ND,EVi : VEKTOR;
+               I         : LONGREAL;
+               i         : CARDINAL;
+               ListPtr   : ListenZeiger;
+               II        : CARDVEKTOR;
+BEGIN
+      IF (NSt <= 1) THEN X[1]:=0.0; A[1]:=2.0; Anzahl:=1; RETURN; END;
+      IF (NSt > MaxDim-2) THEN NSt:=MaxDim-2; END;
+      Anzahl:=NSt;
+      FOR i:=1 TO NSt-1 DO
+        I:=VAL(LONGREAL,i);
+        ND[i]:=I / sqrt(4.0*I*I - 1.0);
+        HD[i]:=0.0;
+      END;
+      ND[NSt]:=0.0;
+      HD[NSt]:=0.0;
+      (* Die X[i] sind die Eigenwerte der Trilinearmatrix *)
+      TriQL(HD,ND,X,NSt,ListPtr,genau);
+      QuickSort(X,II,NSt,1);
+      FOR i:=1 TO NSt DO
+        BerechneEV(HD,ND,X[i],EVi,II[i],genau,30,ListPtr,NSt);
+        A[i]:=2.0*EVi[1]*EVi[1];
+<* IF (__DEBUG__) THEN *>
+        TIO.WrLngReal(X[i],24,-16);
+        TIO.WrLngReal(A[i],24,-16);
+        TIO.WrLn;
+<* END *>
+      END;
+END InitGaussInt;
+
+PROCEDURE GaussInt(Untergr,Obergr : LONGREAL;
+                   N              : CARDINAL;
+                   Funk           : Funktion1) : LONGREAL;
+
+          VAR      Sum : LONGREAL;
+                   p,h : LONGREAL;
+                   i,k : CARDINAL;
+BEGIN
+      h:=ABS(Obergr-Untergr) / VAL(LONGREAL,2*N);
+      Sum:=0.0; p:=Untergr + h;
+      FOR i:=1 TO N DO
+        FOR k:=1 TO Anzahl DO Sum:=Sum + A[k]*h*Funk(p + X[k]*h); END;
+        p:=p + 2.0*h;
+      END;
+      RETURN Sum;
+END GaussInt;
+
+BEGIN (* Zur Sicherheit vorbelegt *)
+      Anzahl:=2;
+      X[1]:=0.577350269189626E+00; X[2]:=-X[1];
+      A[1]:=1.000000000000000E+00; A[2]:=A[1];
+END GAUSSINT;
+
+PROCEDURE GaussIntKonst(Untergr,Obergr : LONGREAL;
+                        N              : CARDINAL;
+                        NSt            : CARDINAL; (* Ordnung der Integration *)
+                        Funk           : Funktion1) : LONGREAL;
+
+          VAR Sum : LONGREAL;
+              p,h : LONGREAL;
+              A,X : ARRAY [1..24] OF LONGREAL;
+              i,k : CARDINAL;
+BEGIN
+      IF NOT (NSt IN {4,6,8,10,12,16,24}) THEN
+        ErrOut("Anzahl St��tzpunkte nicht implementiert, nehme 4");
+        NSt:=4;
+      END;
+      CASE NSt OF
+        2: (* GAUSS-02 DAT *)
+        
+        X[1]:=0.577350269189626E+00; X[2]:=-X[1];
+        A[1]:=1.000000000000000E+00; A[2]:=A[1];
+
+      | 4: (* GAUSS-04 DAT *)
+  
+        X[1]:=0.861136311594053E+00; X[4]:=-X[1];
+        X[2]:=0.339981043584856E+00; X[3]:=-X[2];
+  
+        A[1]:=0.347854845137454E+00; A[4]:=A[1];
+        A[2]:=0.652145154862546E+00; A[3]:=A[2];
+  
+      | 6: (* GAUSS-06 DAT *)
+ 
+        A[1]:=0.171324492379170E+00; A[6]:=A[1];
+        A[2]:=0.360761573048139E+00; A[5]:=A[2];
+        A[3]:=0.467913934572691E+00; A[4]:=A[3];
+  
+        X[1]:=0.932469514203152E+00; X[6]:=-X[1];
+        X[2]:=0.661209386466265E+00; X[5]:=-X[2];
+        X[3]:=0.238619186083197E+00; X[4]:=-X[3];
+  
+      | 8: (* GAUSS-08 DAT *)
+  
+        X[1]:=0.960289856497536E+00; X[8]:=-X[1];
+        X[2]:=0.796666477413627E+00; X[7]:=-X[2];
+        X[3]:=0.525532409916329E+00; X[6]:=-X[3];
+        X[4]:=0.183434642495650E+00; X[5]:=-X[4];
+  
+        A[1]:=0.101228536290376E+00; A[8]:=A[1];
+        A[2]:=0.222381034453374E+00; A[7]:=A[2];
+        A[3]:=0.313706645877887E+00; A[6]:=A[3];
+        A[4]:=0.362683783378362E+00; A[5]:=A[4];
+        
+      |10: (* GAUSS-10 DAT *)
+  
+        X[1]:=0.973906528517172E+00; X[10]:=-X[1];
+        X[2]:=0.865063366688985E+00; X[ 9]:=-X[2];
+        X[3]:=0.679409568299024E+00; X[ 8]:=-X[3];
+        X[4]:=0.433395394129247E+00; X[ 7]:=-X[4];
+        X[5]:=0.148874338981631E+00; X[ 6]:=-X[5];
+  
+        A[1]:=0.666713443086880E-01; A[10]:=A[1];
+        A[2]:=0.149451349150581E+00; A[ 9]:=A[2];
+        A[3]:=0.219086362515982E+00; A[ 8]:=A[3];
+        A[4]:=0.269266719309996E+00; A[ 7]:=A[4];
+        A[5]:=0.295524224714753E+00; A[ 6]:=A[5];
+      |12: (* GAUSS-12 DAT *)
+  
+        X[1]:=0.981560634246719E+00; X[12]:=-X[1];
+        X[2]:=0.904117256370475E+00; X[11]:=-X[2];
+        X[3]:=0.769902674194305E+00; X[10]:=-X[3];
+        X[4]:=0.587317954286617E+00; X[ 9]:=-X[4];
+        X[5]:=0.367831498998180E+00; X[ 8]:=-X[5];
+        X[6]:=0.125233408511469E+00; X[ 7]:=-X[6];
+  
+        A[1]:=0.471753363865120E-01; A[12]:=A[1];
+        A[2]:=0.106939325995318E+00; A[11]:=A[2];
+        A[3]:=0.160078328543346E+00; A[10]:=A[3];
+        A[4]:=0.203167426723066E+00; A[ 9]:=A[4];
+        A[5]:=0.233492536538355E+00; A[ 8]:=A[5];
+        A[6]:=0.249147045813403E+00; A[ 7]:=A[6];
+  
+      | 16: (* GAUSS-16 DAT *)
+  
+        X[1]:=0.989400934991649939E+00; X[16]:=-X[1];
+        X[2]:=0.944575023073232573E+00; X[15]:=-X[2];
+        X[3]:=0.865631202387831741E+00; X[14]:=-X[3];
+        X[4]:=0.755404408355003026E+00; X[13]:=-X[4];
+        X[5]:=0.617876244402643743E+00; X[12]:=-X[5];
+        X[6]:=0.458016777657227384E+00; X[11]:=-X[6];
+        X[7]:=0.281603550779258915E+00; X[10]:=-X[7];
+        X[8]:=0.950125098376374405E-01; X[ 9]:=-X[8];
+  
+        A[1]:=0.271524594117540947E-01; A[16]:=A[1];
+        A[2]:=0.622535239386478928E-01; A[15]:=A[2];
+        A[3]:=0.951585116824927857E-01; A[14]:=A[3];
+        A[4]:=0.124628971255533877E+00; A[13]:=A[4];
+        A[5]:=0.149595988816576736E+00; A[12]:=A[5];
+        A[6]:=0.169156519395002536E+00; A[11]:=A[6];
+        A[7]:=0.182603415044923584E+00; A[10]:=A[7];
+        A[8]:=0.189450610455068502E+00; A[ 9]:=A[8];
+  
+      |24: (* GAUSS-24 DAT *)
+
+        A[ 1]:=0.123412297999871993E-01; A[24]:=A[1];
+        A[ 2]:=0.285313886289336634E-01; A[23]:=A[2];
+        A[ 3]:=0.442774388174198060E-01; A[22]:=A[3];
+        A[ 4]:=0.592985849154367807E-01; A[21]:=A[4];
+        A[ 5]:=0.733464814110802998E-01; A[20]:=A[5];
+        A[ 6]:=0.861901615319532743E-01; A[19]:=A[6];
+        A[ 7]:=0.976186521041138844E-01; A[18]:=A[7];
+        A[ 8]:=0.107444270115965634E+00; A[17]:=A[8];
+        A[ 9]:=0.115505668053725599E+00; A[16]:=A[9];
+        A[10]:=0.121670472927803391E+00; A[15]:=A[10];
+        A[11]:=0.125837456346828303E+00; A[14]:=A[11];
+        A[12]:=0.127938195346752159E+00; A[13]:=A[12];
+ 
+        X[ 1]:=0.995187219997021366E+00; X[24]:=-X[1];
+        X[ 2]:=0.974728555971309502E+00; X[23]:=-X[2];
+        X[ 3]:=0.938274552002732756E+00; X[22]:=-X[3];
+        X[ 4]:=0.886415527004401030E+00; X[21]:=-X[4];
+        X[ 5]:=0.820001985973902919E+00; X[20]:=-X[5];
+        X[ 6]:=0.740124191578554358E+00; X[19]:=-X[6];
+        X[ 7]:=0.648093651936975573E+00; X[18]:=-X[7];
+        X[ 8]:=0.545421471388839535E+00; X[17]:=-X[8];
+        X[ 9]:=0.433793507626045141E+00; X[16]:=-X[9];
+        X[10]:=0.315042679696163369E+00; X[15]:=-X[10];
+        X[11]:=0.191118867473616311E+00; X[14]:=-X[11];
+        X[12]:=0.640568928626056300E-01; X[13]:=-X[12];
+      ELSE
+      END;    
+      h:=ABS(Obergr-Untergr) / VAL(LONGREAL,2*N);
+      Sum:=0.0; p:=Untergr + h;
+      FOR i:=1 TO N DO
+        FOR k:=1 TO NSt DO Sum:=Sum + A[k]*h*Funk(p + X[k]*h); END;
+        p:=p + 2.0*h;
+      END;
+      RETURN Sum;
+END GaussIntKonst;
+
+PROCEDURE SimpInt(Untergr,Obergr : LONGREAL;
+                  genau          : LONGREAL;
+                  MaxIter        : CARDINAL;
+              VAR Iter           : CARDINAL;
+                  Funk           : Funktion1) : LONGREAL;
+
+          VAR  i,n                     : CARDINAL;
+               Sum,ss,S1,S2,S4, H,H2,x : LONGREAL;
+BEGIN
+      n:=2; Fehler:=FALSE;
+      H:=(Obergr-Untergr)*0.5;
+      S1:=H*(Funk(Untergr) + Funk(Obergr));
+      S2:=0.0;
+      S4:=4.0*H*Funk(Untergr+H);
+      Sum:=S1*S2*S4;
+      Iter:=0;
+      REPEAT
+        INC(Iter);
+        ss:=Sum;
+        n:=2*n;
+        S2:=0.5*S2 + 0.25*S4;
+        S1:=0.5*S1; H:=0.5*H;
+        i:=1;       S4:=0.0;
+        H2:=2.0*H;  x:=Untergr+H;
+        REPEAT
+          S4:=S4 + Funk(x);
+          x:=x+H2;
+          INC(i,2);
+        UNTIL (i > n);
+        S4:=4.0*H*S4;
+        Sum:=S1 + S2 + S4;
+      UNTIL (ABS(Sum-ss) < genau*ABS(Sum)) OR (Iter >= MaxIter);
+      IF (Iter >= MaxIter) THEN
+        Fehler:=TRUE;
+        Fehlerflag:='MaxIter "uberschritten (SimpInt)';
+        ErrOut(Fehlerflag);
+      END;
+      RETURN Sum / 3.0;
+END SimpInt;
+
+PROCEDURE Romberg(UnterGr  : LONGREAL;
+                  OberGr   : LONGREAL;
+                  genau    : LONGREAL;
+                  MaxIter  : CARDINAL;
+              VAR Iter     : CARDINAL;
+                  Funk     : Funktion1) : LONGREAL;
+
+          CONST  eps          = 1.0E-12;
+                 Max          = 100;
+
+          TYPE   Tableau      = ARRAY [1..Max] OF LONGREAL;
+
+          VAR    Xneu,Xalt    : Tableau;  (* Integrationstableaus *)
+                 i,j,ProdIter : CARDINAL;
+                 dH,Sum,Prod  : LONGREAL; (* dH : Abstand d. Int.-punkte *)
+BEGIN
+      IF (genau <= eps) THEN genau:=eps;  END;
+      IF (MaxIter = 0) THEN MaxIter:=100; END;
+      IF (MaxIter > Max) THEN MaxIter:=Max; END;
+      dH:=OberGr - UnterGr;
+      Xalt[1]:=0.5*dH*(Funk(UnterGr) + Funk(OberGr));
+      Iter:=1; ProdIter:=1;
+      REPEAT
+        INC(Iter);
+        Sum:=0.0;
+        FOR i:=1 TO ProdIter DO
+          Sum:=Sum + Funk(UnterGr + (VAL(LONGREAL,i) - 0.5)*dH);
+        END;
+        Xneu[1]:=0.5*(Xalt[1] + dH*Sum);
+        ProdIter:=2*ProdIter;
+        Prod:=1.0;
+        FOR j:=2 TO Iter DO
+          Prod:=4.0*Prod;
+          Xneu[j]:=(Prod*Xneu[j-1] - Xalt[j-1]) / (Prod - 1.0);
+        END;
+        dH:=0.5*dH;
+        Xalt:=Xneu;
+      UNTIL (ABS(Xneu[Iter - 1] - Xneu[Iter]) <= ABS(genau * Xneu[Iter]))
+         OR (Iter >= MaxIter);
+
+      IF Iter >= MaxIter THEN
+        Fehler:=TRUE; Fehlerflag:='MaxIter "uberschritten (Romberg).';
+        ErrOut(Fehlerflag);
+      END;
+      RETURN Xneu[Iter];
+END Romberg;
+
+PROCEDURE Quadrature(Untergr,Obergr : LONGREAL;
+                     eps            : LONGREAL;    (* Geforderte Genauigkeit *)
+                 VAR err            : LONGREAL;    (* Fehlerabsch"atzung     *)
+                     Funkt          : Funktion1) : LONGREAL;
+
+          CONST  eta                              = 1.0E-15;
+
+          VAR    H,c,d1,den,ddt,e,gr,hm,
+                 nt,sm,t,t1,t2,t2a,ta,tab,tb,v,w  : LONGREAL;
+                 int                              : LONGREAL;
+                 dt                               : ARRAY [0..6] OF LONGREAL;
+                 d                                : ARRAY [1..6] OF LONGREAL;
+                 i,m,mr,n,nn                      : CARDINAL;
+                 bo,bu,odd                        : BOOLEAN;
+BEGIN
+      Fehler:=FALSE;
+      IF (eps < eta) THEN eps:=eta; END;
+      FOR i:=0 TO 6 DO dt[i]:=0.0; END;
+      FOR i:=1 TO 6 DO d[i]:=0.0; END;
+      n:=2; nn:=3;
+      H:=Obergr-Untergr;
+      t1:=0.0; gr:=0.0; sm:=0.0;
+      t2:=0.5*(Funkt(Untergr) + Funkt(Obergr)); t2a:=t2; tb:=ABS(t2);
+      c:=t2*H; dt[0]:=c;
+      odd:=TRUE; bu:=FALSE;
+      m:=0;
+      LOOP
+        INC(m);
+        bo:= (m >= 7);
+        hm:=H/VAL(LONGREAL,n);
+        IF odd THEN
+          FOR i:=1 TO n BY 2 DO
+            w:=Funkt(Untergr+VAL(LONGREAL,i)*hm); (**)
+            t2:=t2+w;
+            tb:=tb+ABS(w);
+          END;
+          nt:=t2;
+          tab:=tb*ABS(hm);
+          d[1]:= 16.0/9.0;
+          d[3]:= 64.0/9.0;
+          d[5]:=256.0/9.0;
+        ELSE
+          FOR i:=1 TO n BY 6 DO
+            w:=VAL(LONGREAL,i)*hm;
+            t1:=Funkt(Untergr+w) + Funkt(Obergr-w) + t1;
+          END;
+          nt:=t1+t2a;
+          t2a:=t2;
+          d[1]:=9.0/4.0;
+          d[3]:=9.0;
+          d[5]:=36.0;
+        END; (* IF odd *)
+        ddt:=dt[0];
+        nt:=nt*hm; t:=nt; dt[0]:=nt;
+        IF bo THEN
+          mr:=6;
+          d[6]:=64.0;
+          w:=144.0;
+        ELSE
+          mr:=m;
+          w:=VAL(LONGREAL,n*n); d[m]:=w;
+        END; (* IF bo *)
+        FOR i:=1 TO mr DO
+          d1:=d[i]*ddt;
+          den:=d1-nt;
+          e:=nt-ddt;
+          IF (den <> 0.0) THEN
+            e:=e/den;
+            v:=nt*e;
+            nt:=d1*e;
+            t:=t+v;
+          ELSE
+            v:=0.0;
+            nt:=0.0;
+          END; (* IF *)
+          ddt:=dt[i];
+          dt[i]:=v;
+        END; (* FOR i *)
+        ta:=c;
+        int:=t; c:=t;
+        IF NOT bo THEN t:=t-v; END;
+        v:=t-ta;
+        t:=t+v;
+        err:=ABS(v);
+        IF (ta < t) THEN
+          d1:=ta;
+          ta:=t;
+          t:=d1;
+        END;
+        bo:=bo AND ((ta < gr) OR (t > sm));
+        IF bu OR bo OR (err < eps*tab*w) THEN EXIT END;
+        gr:=ta;
+        sm:=t;
+        odd:=NOT odd;
+        i:=n;
+        n:=nn;
+        nn:=i*i;
+        bu:=bo;
+        d[2]:=4.0;
+        d[4]:=16.0;
+        IF (m = 15) THEN
+          Fehler:=TRUE;
+          Fehlerflag:='keine Konvergenz (Quadrature)';
+          ErrOut(Fehlerflag);
+          EXIT;
+        END;
+      END; (* LOOP m *)
+      v:=eta*tab;
+      IF (err < v) THEN err:=v; END;
+      IF (ABS(err) > ABS(eps)) AND NOT Fehler THEN
+        Fehler:=TRUE;
+        Fehlerflag:='|err| > |eps| (Quadrature)';
+        ErrOut(Fehlerflag);
+      END;
+      RETURN int;
+END Quadrature;
+
+PROCEDURE NumInt(VAR X,Y : ARRAY OF LONGREAL;
+                     N   : CARDINAL) : LONGREAL;
+
+          VAR i      : CARDINAL;
+              dX,Sum : LONGREAL;
+BEGIN
+      Sum:=0.0;
+      FOR i:=1 TO N-1 DO
+        dX := 0.5*(X[i] -X[i-1]);
+        Sum:=Sum + dX*(Y[i] + Y[i-1]);
+      END;
+      RETURN Sum;
+END NumInt;
+
+PROCEDURE Simson(VAR Y : ARRAY OF LONGREAL;
+                     H : LONGREAL;
+                     N : CARDINAL) : LONGREAL;
+
+          VAR i,j : CARDINAL;
+              sum : LONGREAL;
+BEGIN
+      sum:=0.0;
+      j:=0;
+      FOR i:=2 TO N-1 BY 2 DO
+        sum:=sum + Y[j] + 4.0*Y[i-1] + Y[i];
+        j:=i;
+      END;
+      IF NOT ODD(N) THEN
+        sum:=sum - 0.25*Y[N-3] + 2.0*Y[j] + 1.25*Y[N-1];
+      END;
+      RETURN sum*H / 3.0;
+END Simson;
+ 
+PROCEDURE SimCum(VAR Y    : ARRAY OF LONGREAL;
+                     H    : LONGREAL;
+                     N    : CARDINAL) : LONGREAL;
+
+          VAR i,j,k : CARDINAL;
+              fac,S : LONGREAL;
+              term  : ARRAY [1..5] OF LONGREAL;
+              Area  : POINTER TO ARRAY [0..MAX(INTEGER)-1] OF LONGREAL;
+BEGIN
+      ALLOCATE(Area,N*TSIZE(LONGREAL));
+      Area^[0] := 0.0;
+      fac      := H / 3.0;
+      term[1]  :=   1.25*Y[0];
+      term[2]  := - 0.25*Y[0];
+      k:=0;
+      FOR i:=2 TO N-1 BY 2 DO
+        j := i - 1;
+        term[3]  :=   2.00*Y[j];
+        term[4]  :=   1.25*Y[i];
+        term[5]  := - 0.25*Y[i];
+        Area^[j] := Area^[k] + fac*(term[1] + term[3] + term[5]);
+        Area^[i] := Area^[j] + fac*(term[2] + term[3] + term[4]);
+        term[1]  := term[4];
+        term[2]  := term[5];
+        k := i;
+      END; (* FOR *)
+      IF NOT ODD(N) THEN
+        Area^[N-1]:=Area^[k] + fac*(-0.25*Y[j] + 2.0*Y[k] + 1.25*Y[N-1]);
+      END;
+      S := Area^[N-1];
+      DEALLOCATE(Area,N*TSIZE(LONGREAL));
+      RETURN S;
+END SimCum;
+
+END Integral.