--- 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.