|
a/SpezFunkt2.mod.m2pp |
|
b/SpezFunkt2.mod.m2pp |
|
... |
|
... |
| 45 |
(* 14.02.18, MRi: Erstellen der ersten Version von CdGamma *)
|
45 |
(* 14.02.18, MRi: Erstellen der ersten Version von CdGamma *)
|
| 46 |
(* 23.05.18, MRi: Erweiterung von CdGamma und Aufnahme in SpezFunkt2 *)
|
46 |
(* 23.05.18, MRi: Erweiterung von CdGamma und Aufnahme in SpezFunkt2 *)
|
| 47 |
(* 19.12.18, MRi: Ersetzen von GammaIn durch Neuuebersetzng von AS239 *)
|
47 |
(* 19.12.18, MRi: Ersetzen von GammaIn durch Neuuebersetzng von AS239 *)
|
| 48 |
(* da die alte Routine (basierend auf AS32) fuer grosse X *)
|
48 |
(* da die alte Routine (basierend auf AS32) fuer grosse X *)
|
| 49 |
(* teilweise falsche Ergebnisse erzielt hatte. *)
|
49 |
(* teilweise falsche Ergebnisse erzielt hatte. *)
|
|
|
50 |
(* 19.12.18, MRi: Erstellen der ersten Version von AlNorm *)
|
| 50 |
(*------------------------------------------------------------------------*)
|
51 |
(*------------------------------------------------------------------------*)
|
| 51 |
(* Testroutinen vorhanden, IBeta gegen unabhaengigen Fortran-Code ge- *)
|
52 |
(* Testroutinen vorhanden, IBeta gegen unabhaengigen Fortran-Code ge- *)
|
| 52 |
(* testet, IGamma & JGamma nur gegen Pascal-Aequivalent *)
|
53 |
(* testet, IGamma & JGamma nur gegen Pascal-Aequivalent *)
|
| 53 |
(* Testcode (parially against Fortran versions) in TstSpzFun2 *)
|
54 |
(* Testcode (parially against Fortran versions) in TstSpzFun2 *)
|
| 54 |
(* dGamit wird in BoysInt2Lib intensiv genutzt *)
|
55 |
(* dGamit wird in BoysInt2Lib intensiv genutzt *)
|
|
... |
|
... |
| 60 |
(*------------------------------------------------------------------------*)
|
61 |
(*------------------------------------------------------------------------*)
|
| 61 |
(* Implementation : Michael Riedl *)
|
62 |
(* Implementation : Michael Riedl *)
|
| 62 |
(* Licence : GNU Lesser General Public License (LGPL) *)
|
63 |
(* Licence : GNU Lesser General Public License (LGPL) *)
|
| 63 |
(*------------------------------------------------------------------------*)
|
64 |
(*------------------------------------------------------------------------*)
|
| 64 |
|
65 |
|
| 65 |
(* $Id: SpezFunkt2.mod.m2pp,v 1.4 2018/02/17 09:17:19 mriedl Exp $ *)
|
66 |
(* $Id: SpezFunkt2.mod.m2pp,v 1.6 2019/02/01 22:24:03 mriedl Exp mriedl $ *)
|
| 66 |
|
67 |
|
| 67 |
FROM Errors IMPORT FatalError,ErrOut,Fehler,Fehlerflag;
|
68 |
FROM Errors IMPORT FatalError,ErrOut,Fehler,Fehlerflag;
|
| 68 |
IMPORT LowLong;
|
69 |
IMPORT LowLong;
|
| 69 |
FROM TestReal IMPORT Real8NaNquite,IsNearInteger;
|
70 |
FROM TestReal IMPORT Real8NaNquite,IsNearInteger;
|
| 70 |
IMPORT TestReal;
|
71 |
IMPORT TestReal;
|
|
... |
|
... |
| 104 |
|
105 |
|
| 105 |
(* machine dependent constants *)
|
106 |
(* machine dependent constants *)
|
| 106 |
|
107 |
|
| 107 |
CONST ltone = 7.0; utzero = 18.66; con = 1.28;
|
108 |
CONST ltone = 7.0; utzero = 18.66; con = 1.28;
|
| 108 |
|
109 |
|
| 109 |
p = 0.398942280444; q = 0.39990348504;
|
110 |
p = 0.398942280444; q = 0.399903485040;
|
| 110 |
r = 0.398942280385;
|
111 |
r = 0.398942280385;
|
| 111 |
a1 = 5.75885480458; a2 = 2.62433121679;
|
112 |
a1 = 5.758854804580; a2 = 2.624331216790;
|
| 112 |
a3 = 5.92885724438;
|
113 |
a3 = 5.928857244380;
|
| 113 |
b1 = -29.8213557807; b2 = 48.6959930692;
|
114 |
b1 = -29.821355780700; b2 = 48.695993069200;
|
|
|
115 |
|
| 114 |
c1 = -3.8052E-08; c2 = 3.98064794E-04;
|
116 |
c1 = -3.80520000E-08; c2 = 3.98064794E-04;
|
| 115 |
c3 = -0.151679116635; c4 = 4.8385912808;
|
117 |
c3 = -0.151679116635; c4 = 4.838591280800;
|
| 116 |
c5 = 0.742380924027; c6 = 3.99019417011;
|
118 |
c5 = 0.742380924027; c6 = 3.990194170110;
|
|
|
119 |
|
| 117 |
d1 = 1.00000615302; d2 = 1.98615381364;
|
120 |
d1 = 1.000006153020; d2 = 1.986153813640;
|
| 118 |
d3 = 5.29330324926; d4 = -15.1508972451;
|
121 |
d3 = 5.293303249260; d4 = -15.150897245100;
|
| 119 |
d5 = 30.789933034;
|
122 |
d5 = 30.789933034000;
|
| 120 |
|
123 |
|
| 121 |
VAR z,y,alnorm : LONGREAL;
|
124 |
VAR z,y,alnorm : LONGREAL;
|
| 122 |
up : BOOLEAN;
|
125 |
up : BOOLEAN;
|
| 123 |
BEGIN
|
126 |
BEGIN
|
| 124 |
up := Upper;
|
127 |
up := Upper;
|
|
... |
|
... |
| 149 |
PROCEDURE GammaIn( X : LONGREAL;
|
152 |
PROCEDURE GammaIn( X : LONGREAL;
|
| 150 |
P : LONGREAL;
|
153 |
P : LONGREAL;
|
| 151 |
VAR Ifault : INTEGER) : LONGREAL;
|
154 |
VAR Ifault : INTEGER) : LONGREAL;
|
| 152 |
|
155 |
|
| 153 |
(*----------------------------------------------------------------*)
|
156 |
(*----------------------------------------------------------------*)
|
| 154 |
(* Auxiliary functions required: ALOGAM := logarithm of the gamma *)
|
157 |
(* This is a M2 verions of AS239. Auxiliary functions required *)
|
| 155 |
(* function, and AlNorm := algorithm AS66 *)
|
158 |
(* are the logarithm of the gamma function and AlNorm (AS66). *)
|
| 156 |
(*----------------------------------------------------------------*)
|
159 |
(*----------------------------------------------------------------*)
|
| 157 |
|
160 |
|
| 158 |
CONST Tol = 1.0E-14;
|
161 |
CONST Tol = 1.0E-14;
|
| 159 |
Xbig = 1.0E+08; Oflow = MAX(REAL);
|
162 |
Xbig = 1.0E+08; Oflow = LFLOAT(MAX(REAL));
|
| 160 |
Elimit = -88.0; Plimit = 1000.0;
|
163 |
Elimit = -88.0; Plimit = 1000.0;
|
| 161 |
|
164 |
|
| 162 |
VAR pn1,pn2,pn3,pn4,pn5,pn6 : LONGREAL;
|
165 |
VAR pn1,pn2,pn3,pn4,pn5,pn6 : LONGREAL;
|
| 163 |
a,b,c,an,rn,arg,gamma : LONGREAL;
|
166 |
a,b,c,an,rn,arg,gamma : LONGREAL;
|
| 164 |
BEGIN
|
167 |
BEGIN
|
|
... |
|
... |
| 171 |
IF (X = 0.0) THEN
|
174 |
IF (X = 0.0) THEN
|
| 172 |
RETURN 0.0;
|
175 |
RETURN 0.0;
|
| 173 |
END; (* IF *)
|
176 |
END; (* IF *)
|
| 174 |
(* Use a normal approximation if P > Plimit *)
|
177 |
(* Use a normal approximation if P > Plimit *)
|
| 175 |
IF (P > Plimit) THEN
|
178 |
IF (P > Plimit) THEN
|
| 176 |
pn1 := 3.0*sqrt(P)*((X/P)**(1.0/3.0)+1.0/(9.0*P)-1.0);
|
179 |
pn1 := 3.0*sqrt(P)*(power((X/P),(1.0/3.0)) + 1.0/(9.0*P) - 1.0);
|
| 177 |
RETURN AlNorm(pn1,FALSE);
|
180 |
RETURN AlNorm(pn1,FALSE);
|
| 178 |
END; (* IF *)
|
181 |
END; (* IF *)
|
| 179 |
(* If X is extremely large compared to P then set gamma = 1 *)
|
182 |
(* If X is extremely large compared to P then set gamma = 1 *)
|
| 180 |
IF (X > Xbig) THEN
|
183 |
IF (X > Xbig) THEN
|
| 181 |
RETURN 1.0;
|
184 |
RETURN 1.0;
|