| 1 | ! ================================================================================================================================= |
|---|
| 2 | ! MODULE : gammad_inc |
|---|
| 3 | ! |
|---|
| 4 | ! CONTACT : orchidee-help _at_ ipsl.jussieu.fr |
|---|
| 5 | ! |
|---|
| 6 | ! LICENCE : IPSL (2006) |
|---|
| 7 | ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
|---|
| 8 | ! |
|---|
| 9 | !>\BRIEF This module computes the parameters for TOPMODEL, gamma distribution fitting. |
|---|
| 10 | !! |
|---|
| 11 | !!\n DESCRIPTION : contains. |
|---|
| 12 | !! |
|---|
| 13 | !! RECENT CHANGE(S) : None |
|---|
| 14 | !! |
|---|
| 15 | !! REFERENCE(S) : None |
|---|
| 16 | !! |
|---|
| 17 | !! SVN : |
|---|
| 18 | !! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/perso/shushi.peng/ORCHIDEE/src_sechiba/gammad_inc.f90 $ |
|---|
| 19 | !! $Date: 2015-03-10 10:16:04 +0100 (Tue, 10 Mar 2015) $ |
|---|
| 20 | !! $Revision: 2535 $ |
|---|
| 21 | !! \n |
|---|
| 22 | !_ ================================================================================================================================ |
|---|
| 23 | |
|---|
| 24 | !C###################################################################### |
|---|
| 25 | !C###################################################################### |
|---|
| 26 | !C |
|---|
| 27 | !C THIS ROUTINE COMPUTE THE INCOMPLETE AND THE COMPLEMENTARY |
|---|
| 28 | !C INCOMPLETE GAMMA FUNCTION DOUBLE PRECISION. IN BRIEF, THIS |
|---|
| 29 | !C FUNCTION HELP TO COMPUTE THE DIFFERENT FRACTION PRESENT IN THE |
|---|
| 30 | !C TOPMODEL FRAMEWORK. |
|---|
| 31 | !C MORE EXPLANATION ON THE SUBROUTINE USE ARE GIVEN IN THE NEXT |
|---|
| 32 | !C COMMENTARY |
|---|
| 33 | !C |
|---|
| 34 | !C THIS ROUTINE WAS FOUND ON http://www.netlib.org WEB SITE. |
|---|
| 35 | !C |
|---|
| 36 | !C REFERENCE - W. GAUTSCHI, :: A COMPUTATIONAL PROCEDURE FOR INCOMPLETE |
|---|
| 37 | !C GAMMA FUNCTIONS, ACM TRANS. MATH. SOFTWARE., (1979) 482-489 |
|---|
| 38 | !C |
|---|
| 39 | !C###################################################################### |
|---|
| 40 | !C###################################################################### |
|---|
| 41 | !C |
|---|
| 42 | !C LET GAMMA(A) DENOTE THE GAMMA FUNCTION AND GAM(A,X) THE |
|---|
| 43 | !C (COMPLEMENTARY) INCOMPLETE GAMMA FUNCTION, |
|---|
| 44 | !C |
|---|
| 45 | !C GAM(A,X)=INTEGRAL FROM T=X TO T=INFINITY OF EXP(-T)*T**(A-1). |
|---|
| 46 | !C |
|---|
| 47 | !C LET GAMSTAR(A,X) DENOTE TRICOMI:S FORM OF THE INCOMPLETE GAMMA |
|---|
| 48 | !C FUNCTION, WHICH FOR A.GT.0. IS DEFINED BY |
|---|
| 49 | !C |
|---|
| 50 | !C GAMSTAR(A,X)=(X**(-A)/GAMMA(A))*INTEGRAL FROM T=0 TO T=X OF |
|---|
| 51 | !C EXP(-T)*T**(A-1), |
|---|
| 52 | !C |
|---|
| 53 | !C AND FOR A.LE.0. BY ANALYTIC CONTINUATION. FOR THE PURPOSE OF |
|---|
| 54 | !C THIS SUBROUTINE, THESE FUNCTIONS ARE NORMALIZED AS FOLLOWS& |
|---|
| 55 | !C |
|---|
| 56 | !C GAM(A,X)/GAMMA(A), IF A.GT.0., |
|---|
| 57 | !C G(A,X)= |
|---|
| 58 | !C EXP(X)*X**(-A)*GAM(A,X), IF A.LE.0., |
|---|
| 59 | !C |
|---|
| 60 | !C GSTAR(A,X)=(X**A)*GAMSTAR(A,X) |
|---|
| 61 | !C =(1/GAMMA(A))*INTEGRAL FROM T=0 TO T=X OF EXP(-T)*T**(A-1) |
|---|
| 62 | !C |
|---|
| 63 | !C THE PROGRAM BELOW ATTEMPTS TO EVALUATE G(A,X) AND GSTAR(A,X), |
|---|
| 64 | !C BOTH TO AN ACCURACY OF ACC SIGNIFICANT DECIMAL DIGITS, FOR ARBI- |
|---|
| 65 | !C TRARY REAL VALUES OF A AND NONNEGATIVE VALUES OF X. THE SUB- |
|---|
| 66 | !C ROUTINE AUTOMATICALLY CHECKS FOR UNDERFLOW AND OVERFLOW CONDI- |
|---|
| 67 | !C TIONS AND RETURNS APPROPRIATE WARNINGS THROUGH THE OUTPUT PARA- |
|---|
| 68 | !C METERS IFLG, IFLGST. A RESULT THAT UNDERFLOWS IS RETURNED WITH |
|---|
| 69 | !C THE VALUE 0., ONE THAT OVERFLOWS WITH THE VALUE OF THE LARGEST |
|---|
| 70 | !C MACHINE-REPRESENTABLE NUMBER. |
|---|
| 71 | !C |
|---|
| 72 | !C NEAR LINES IN THE (A,X)-PLANE, A.LT.0., ALONG WHICH GSTAR |
|---|
| 73 | !C VANISHES, THE ACCURACY SPECIFIED WILL BE ATTAINED ONLY IN TERMS |
|---|
| 74 | !C OF ABSOLUTE ERROR, NOT RELATIVE ERROR. THERE ARE OTHER (RARE) |
|---|
| 75 | !C INSTANCES IN WHICH THE ACCURACY ATTAINED IS SOMEWHAT LESS THAN |
|---|
| 76 | !C THE ACCURACY SPECIFIED. THE DISCREPANCY, HOWEVER, SHOULD NEVER |
|---|
| 77 | !C EXCEED ONE OR TWO (DECIMAL) ORDERS OF ACCURACY# NO INDICATION |
|---|
| 78 | !C OF THIS IS GIVEN THROUGH ERROR FLAGS. |
|---|
| 79 | !C |
|---|
| 80 | !C PARAMETER LIST& |
|---|
| 81 | !C |
|---|
| 82 | !C A - - THE FIRST ARGUMENT OF G AND GSTAR. TYPE REAL. |
|---|
| 83 | !C X - - THE SECOND ARGUMENT OF G AND GSTAR. TYPE REAL. |
|---|
| 84 | !C ACC - - THE NUMBER OF CORRECT SIGNIFICANT DECIMAL DIGITS |
|---|
| 85 | !C DESIRED IN THE RESULTS. TYPE REAL. |
|---|
| 86 | !C G - - AN OUTPUT VARIABLE RETURNING THE VALUE OF G(A,X). |
|---|
| 87 | !C TYPE REAL. |
|---|
| 88 | !C GSTAR - - AN OUTPUT VARIABLE RETURNING THE VALUE OF |
|---|
| 89 | !C GSTAR(A,X). TYPE REAL. |
|---|
| 90 | !C IFLG - - AN ERROR FLAG INDICATING A NUMBER OF ERROR CONDI- |
|---|
| 91 | !C TIONS IN G UPON EXECUTION. TYPE INTEGER. |
|---|
| 92 | !C! IFLGST - - AN ERROR FLAG INDICATING A NUMBER OF ERROR CONDI- |
|---|
| 93 | !C TIONS IN GSTAR UPON EXECUTION. TYPE INTEGER. |
|---|
| 94 | !C THE VALUES OF IFLG AND IFLGST HAVE THE FOLLOWING |
|---|
| 95 | !C MEANINGS& |
|---|
| 96 | !C 0 - NO ERROR CONDITION. |
|---|
| 97 | !C 1 - ILLEGAL NEGATIVE ARGUMENT X. THE ROUTINE EXITS |
|---|
| 98 | !C WITH THE VALUES ZERO FOR G AND GSTAR. |
|---|
| 99 | !C 2 - INFINITELY LARGE RESULT AT X=0. THE ROUTINE |
|---|
| 100 | !C RETURNS THE LARGEST MACHINE-REPRESENTABLE NUMBER. |
|---|
| 101 | !C 3 - (ONLY FOR IFLGST) GSTAR IS INDETERMINATE AT |
|---|
| 102 | !C A=0. AND X=0. THE ROUTINE RETURNS THE VALUE 1., |
|---|
| 103 | !C WHICH IS THE LIMIT VALUE AS X TENDS TO ZERO FOR |
|---|
| 104 | !C FIXED A=0. |
|---|
| 105 | !C 4 - THE RESULT UNDERFLOWS. IT IS SET EQUAL TO 0. |
|---|
| 106 | !C 5 - THE RESULT OVERFLOWS. IT IS SET EQUAL TO THE |
|---|
| 107 | !C LARGEST MACHINE-REPRESENTABLE NUMBER, WITH |
|---|
| 108 | !C PROPER SIGN. |
|---|
| 109 | !C 6 - CONVERGENCE FAILS WITHIN 600 ITERATIONS, EITHER |
|---|
| 110 | !C IN TAYLOR:S SERIES OR IN LEGENDRE:S CONTINUED |
|---|
| 111 | !C FRACTION. REASON UNKNOWN. THE COMPUTATION IS |
|---|
| 112 | !C ABORTED AND THE ROUTINE RETURNS WITH ZERO |
|---|
| 113 | !C VALUES FOR G AND GSTAR. |
|---|
| 114 | !C |
|---|
| 115 | !C ALL MACHINE-DEPENDENT PARAMETERS ARE COLLECTED IN THE FIRST |
|---|
| 116 | !C DATA DECLARATION. THEY ARE AS FOLLOWS& |
|---|
| 117 | !C |
|---|
| 118 | !C PREC - - THE SINGLE PRECISION ACCURACY, TO BE SET APPROXI- |
|---|
| 119 | !C MATELY EQUAL TO BETA*ALOG(2.)/ALOG(10.), WHERE BETA |
|---|
| 120 | !C IS THE NUMBER OF BINARY DIGITS AVAILABLE IN THE MAN- |
|---|
| 121 | !C TISSA OF THE SINGLE PRECISION FLOATING-POINT WORD. |
|---|
| 122 | !C TYPE REAL. |
|---|
| 123 | !C TOPEXP - - APPROXIMATELY THE LARGEST POSITIVE NUMBER T SUCH |
|---|
| 124 | !C THAT 10.**T IS STILL REPRESENTABLE ON THE COMPUTER |
|---|
| 125 | !C IN SINGLE PRECISION FLOATING-POINT ARITHMETIC. |
|---|
| 126 | !C TYPE REAL. |
|---|
| 127 | !C BOTEXP - - APPROXIMATELY THE SMALLEST NEGATIVE NUMBER T SUCH |
|---|
| 128 | !C THAT 10.**T IS STILL REPRESENTABLE ON THE COMPUTER |
|---|
| 129 | !C IN SINGLE PRECISION FLOATING-POINT ARITHMETIC. |
|---|
| 130 | !C TYPE REAL. |
|---|
| 131 | !C |
|---|
| 132 | !C IN THE PROGRAM BELOW THESE PARAMETERS ARE SET TO CORRESPOND TO |
|---|
| 133 | !C THE MACHINE CHARACTERISTICS OF THE CDC 6500 COMPUTER. |
|---|
| 134 | !C |
|---|
| 135 | !C THE SECOND DATA DECLARATION CONTAINS THE SINGLE PRECISION |
|---|
| 136 | !C VALUE OF ALOG(10.). THE NEXT DATA DECLARATION CONTAINS THE SUCCES- |
|---|
| 137 | !C SIVE COEFFICIENTS IN THE MACLAURIN EXPANSION OF (1/GAMMA(A+1))-1. |
|---|
| 138 | !C THEY ARE GIVEN TO AS MANY DECIMAL PLACES AS IS NECESSARY TO ACHIEVE |
|---|
| 139 | !C MACHINE PRECISION (ON THE CDC 6500 COMPUTER) IN THE RANGE |
|---|
| 140 | !C ABS(A).LE..5. MORE ACCURATE VALUES OF THESE COEFFICIENTS (TO |
|---|
| 141 | !C 31 DECIMAL PLACES) CAN BE FOUND IN TABLE 5 OF J.W.WRENCH,JR., |
|---|
| 142 | !C CONCERNING TWO SERIES FOR THE GAMMA FUNCTION, MATH. COMPUT. |
|---|
| 143 | !C 22, 1968, 617-626. |
|---|
| 144 | !C |
|---|
| 145 | !C THE SUBROUTINE CALLS ON A FUNCTION SUBROUTINE, NAMED ALGA, |
|---|
| 146 | !C WHICH IS TO PROVIDE SINGLE PRECISION VALUES OF THE LOGARITHM OF |
|---|
| 147 | !C THE GAMMA FUNCTION FOR ARGUMENTS LARGER THAN OR EQUAL TO .5. |
|---|
| 148 | !C A POSSIBLE VERSION OF SUCH A FUNCTION SUBROUTINE IS APPENDED |
|---|
| 149 | !C TO THE PRESENT SUBROUTINE. IT IS TAYLORED TO THE ACCURACY RE- |
|---|
| 150 | !C QUIREMENTS OF THE CDC 6500 COMPUTER, AND USES RATIONAL APPROXI- |
|---|
| 151 | !C MATIONS DUE TO CODY AND HILLSTROM (MATH. COMPUT. 21, 1967, 198- |
|---|
| 152 | !C 203). |
|---|
| 153 | !C |
|---|
| 154 | !C REFERENCE - W. GAUTSCHI, ::A COMPUTATIONAL PROCEDURE FOR |
|---|
| 155 | !C INCOMPLETE GAMMA FUNCTIONS, ACM TRANS. MATH. SOFTWARE. |
|---|
| 156 | !C |
|---|
| 157 | !! ################ |
|---|
| 158 | MODULE gammad_inc |
|---|
| 159 | ! ################ |
|---|
| 160 | ! |
|---|
| 161 | ! |
|---|
| 162 | INTERFACE |
|---|
| 163 | ! |
|---|
| 164 | ! |
|---|
| 165 | SUBROUTINE DGAM(A, X, ACC, G, GSTAR, IFLG, IFLGST) |
|---|
| 166 | |
|---|
| 167 | DOUBLE PRECISION :: A, X, G, GSTAR, C, AL10, ALX, ALPHA, ALPREC, & |
|---|
| 168 | & TOP, BOT, AINF, EPS, EPS1, ES, SGA, AE, AA, AP1, AM1, AEP1, & |
|---|
| 169 | & AEM1, FMA, AEPS, SGAE, AAEPS, ALGP1, SGGA, ALGEP1, SGGS, ALGS, & |
|---|
| 170 | & ALG, SUMM, GA, Y, TERM, U, P, Q, R, V, T, H, SGT, A1X, RHO, XPA, & |
|---|
| 171 | & XMA, S, DLGA |
|---|
| 172 | ! |
|---|
| 173 | END |
|---|
| 174 | ! |
|---|
| 175 | END INTERFACE |
|---|
| 176 | ! |
|---|
| 177 | END MODULE gammad_inc |
|---|
| 178 | ! |
|---|
| 179 | !################################################################ |
|---|
| 180 | ! |
|---|
| 181 | SUBROUTINE DGAM(A, X, ACC, G, GSTAR, IFLG, IFLGST) |
|---|
| 182 | ! |
|---|
| 183 | !################################################################ |
|---|
| 184 | DOUBLE PRECISION :: A, X, G, GSTAR, C, AL10, ALX, ALPHA, ALPREC, & |
|---|
| 185 | & TOP, BOT, AINF, EPS, EPS1, ES, SGA, AE, AA, AP1, AM1, AEP1, & |
|---|
| 186 | & AEM1, FMA, AEPS, SGAE, AAEPS, ALGP1, SGGA, ALGEP1, SGGS, ALGS, & |
|---|
| 187 | & ALG, SUMM, GA, Y, TERM, U, P, Q, R, V, T, H, SGT, A1X, RHO, XPA, & |
|---|
| 188 | & XMA, S, DLGA |
|---|
| 189 | DIMENSION C(29) |
|---|
| 190 | DATA PREC, TOPEXP, BOTEXP /15.9545,308.,-308./ |
|---|
| 191 | DATA AL10 /2.3025850929940456840179914547D0/ |
|---|
| 192 | DATA C /.57721566490153286060651209008D0, & |
|---|
| 193 | & -.65587807152025388107701951515D0, & |
|---|
| 194 | & -4.200263503409523552900393488D-2, & |
|---|
| 195 | & .1665386113822914895017007951D0, & |
|---|
| 196 | & -4.21977345555443367482083013D-2, & |
|---|
| 197 | & -9.6219715278769735621149217D-3, & |
|---|
| 198 | & 7.2189432466630995423950103D-3, & |
|---|
| 199 | & -1.165167591859065112113971D-3, & |
|---|
| 200 | & -2.15241674114950972815730D-4, & |
|---|
| 201 | & 1.2805028238811618615320D-4, & |
|---|
| 202 | & -2.013485478078823865569D-5, & |
|---|
| 203 | & -1.2504934821426706573D-6, & |
|---|
| 204 | & 1.1330272319816958824D-6, & |
|---|
| 205 | & -2.056338416977607103D-7, & |
|---|
| 206 | & 6.1160951044814158D-9, & |
|---|
| 207 | & 5.0020076444692229D-9, & |
|---|
| 208 | & -1.181274570487020D-9, & |
|---|
| 209 | & 1.04342671169110D-10, & |
|---|
| 210 | & 7.782263439905D-12, & |
|---|
| 211 | & -3.696805618642D-12, & |
|---|
| 212 | & 5.1003702875D-13, & |
|---|
| 213 | & -2.058326054D-14, -5.34812254D-15, & |
|---|
| 214 | & 1.2267786D-15, -1.181259D-16,1.187D-18, & |
|---|
| 215 | & 1.412D-18,-2.30D-19,1.7D-20/ |
|---|
| 216 | G = 0.D0 |
|---|
| 217 | GSTAR = 0.D0 |
|---|
| 218 | IF (X.LT.0.D0) GO TO 290 |
|---|
| 219 | !C |
|---|
| 220 | !C INITIALIZATION |
|---|
| 221 | !C |
|---|
| 222 | IFLG = 0 |
|---|
| 223 | IFLGST = 0 |
|---|
| 224 | I = 0 |
|---|
| 225 | IF (X.GT.0.D0) ALX = DLOG(X) |
|---|
| 226 | ALPHA = X + .25D0 |
|---|
| 227 | IF (X.LT..25D0 .AND. X.GT.0.D0) ALPHA = DLOG(.5D0)/ALX |
|---|
| 228 | !C ALPREC = AL10*DBLE(PREC) |
|---|
| 229 | !C TOP = AL10*DBLE(TOPEXP) |
|---|
| 230 | !C BOT = AL10*DBLE(BOTEXP) |
|---|
| 231 | ALPREC =DBLE(DIGITS(1.D0))*DLOG(2.D0) |
|---|
| 232 | TOP =DLOG10(HUGE(1.D0)) |
|---|
| 233 | BOT =DLOG10(TINY(1.D0)) |
|---|
| 234 | AINF = 10.D0**TOPEXP |
|---|
| 235 | EPS = .5D0*10.D0**(-ACC) |
|---|
| 236 | EPS1 = EPS/1.D2 |
|---|
| 237 | SGA = 1.D0 |
|---|
| 238 | IF (A.LT.0.D0) SGA = -SGA |
|---|
| 239 | AE = A |
|---|
| 240 | AA = DABS(A) |
|---|
| 241 | AP1 = A + 1.D0 |
|---|
| 242 | AEP1 = AP1 |
|---|
| 243 | MA = SNGL(.5D0-A) |
|---|
| 244 | FMA = DBLE(FLOAT(MA)) |
|---|
| 245 | AEPS = A + FMA |
|---|
| 246 | SGAE = 1.D0 |
|---|
| 247 | IF (AEPS.LT.0.D0) SGAE = -SGAE |
|---|
| 248 | AAEPS = DABS(AEPS) |
|---|
| 249 | ALGP1 = 0.D0 |
|---|
| 250 | !C |
|---|
| 251 | !C EVALUATION OF THE LOGARITHM OF THE ABSOLUTE VALUE OF |
|---|
| 252 | !C GAMMA(A+1.) AND DETERMINATION OF THE SIGN OF GAMMA(A+1.) |
|---|
| 253 | !C |
|---|
| 254 | SGGA = 1.D0 |
|---|
| 255 | IF (MA.LE.0) GO TO 10 |
|---|
| 256 | IF (AEPS.EQ.0.D0) GO TO 20 |
|---|
| 257 | SGGA = SGAE |
|---|
| 258 | IF (MA.EQ.2*(MA/2)) SGGA = -SGGA |
|---|
| 259 | ALGP1 = DLGA(AEPS+1.D0) - DLOG(AAEPS) |
|---|
| 260 | IF (MA.EQ.1) GO TO 20 |
|---|
| 261 | ALGP1 = ALGP1 + DLGA(1.D0-AEPS) - DLGA(FMA-AEPS) |
|---|
| 262 | GO TO 20 |
|---|
| 263 | 10 ALGP1 = DLGA(AP1) |
|---|
| 264 | 20 ALGEP1 = ALGP1 |
|---|
| 265 | IF (X.GT.0.D0) GO TO 60 |
|---|
| 266 | !C |
|---|
| 267 | !C EVALUATION OF GSTAR(A,0.) AND G(A,0.) |
|---|
| 268 | !C |
|---|
| 269 | IF (A) 30, 40, 50 |
|---|
| 270 | 30 IFLGST = 2 |
|---|
| 271 | GSTAR = AINF |
|---|
| 272 | G = 1.D0/AA |
|---|
| 273 | RETURN |
|---|
| 274 | 40 IFLGST = 3 |
|---|
| 275 | GSTAR = 1.D0 |
|---|
| 276 | IFLG = 2 |
|---|
| 277 | G = AINF |
|---|
| 278 | RETURN |
|---|
| 279 | 50 G = 1.D0 |
|---|
| 280 | RETURN |
|---|
| 281 | 60 IF (A.GT.ALPHA) GO TO 220 |
|---|
| 282 | IF (X.GT.1.5D0) GO TO 240 |
|---|
| 283 | IF (A.LT.-.5D0) GO TO 170 |
|---|
| 284 | !C |
|---|
| 285 | !C DIRECT EVALUATION OF G(A,X) AND GSTAR(A,X) FOR X.LE.1.5 |
|---|
| 286 | !C AND -.5.LE.A.LE.ALPHA(X) |
|---|
| 287 | !C |
|---|
| 288 | GSTAR = 1.D0 |
|---|
| 289 | IF (A.GE..5D0) GO TO 110 |
|---|
| 290 | 70 SUMM = C(29) |
|---|
| 291 | DO 80 K=1,28 |
|---|
| 292 | K1 = 29 - K |
|---|
| 293 | SUMM = AE*SUMM + C(K1) |
|---|
| 294 | 80 CONTINUE |
|---|
| 295 | GA = -SUMM/(1.D0+AE*SUMM) |
|---|
| 296 | Y = AE*ALX |
|---|
| 297 | IF (DABS(Y).GE.1.D0) GO TO 100 |
|---|
| 298 | SUMM = 1.D0 |
|---|
| 299 | TERM = 1.D0 |
|---|
| 300 | K = 1 |
|---|
| 301 | 90 K = K + 1 |
|---|
| 302 | IF (K.GT.600) GO TO 330 |
|---|
| 303 | TERM = Y*TERM/DBLE(FLOAT(K)) |
|---|
| 304 | SUMM = SUMM + TERM |
|---|
| 305 | IF (DABS(TERM).GT.EPS1*SUMM) GO TO 90 |
|---|
| 306 | U = GA - SUMM*ALX |
|---|
| 307 | GO TO 120 |
|---|
| 308 | 100 U = GA - (DEXP(Y)-1.D0)/AE |
|---|
| 309 | GO TO 120 |
|---|
| 310 | 110 U = DEXP(DLGA(A)) - (X**A)/A |
|---|
| 311 | 120 P = AE*X |
|---|
| 312 | Q = AEP1 |
|---|
| 313 | R = AE + 3.D0 |
|---|
| 314 | TERM = 1.D0 |
|---|
| 315 | SUMM = 1.D0 |
|---|
| 316 | K = 1 |
|---|
| 317 | 130 K = K + 1 |
|---|
| 318 | IF (K.GT.600) GO TO 330 |
|---|
| 319 | P = P + X |
|---|
| 320 | Q = Q + R |
|---|
| 321 | R = R + 2.D0 |
|---|
| 322 | TERM = -P*TERM/Q |
|---|
| 323 | SUMM = SUMM + TERM |
|---|
| 324 | IF (DABS(TERM).GT.EPS1*SUMM) GO TO 130 |
|---|
| 325 | V = (X**AEP1)*SUMM/AEP1 |
|---|
| 326 | G = U + V |
|---|
| 327 | IF (I.EQ.1) GO TO 180 |
|---|
| 328 | IF (A) 140, 150, 160 |
|---|
| 329 | 140 T = DEXP(X)*X**(-A) |
|---|
| 330 | G = T*G |
|---|
| 331 | GSTAR = 1.D0 - A*G*DEXP(-ALGP1)/T |
|---|
| 332 | RETURN |
|---|
| 333 | 150 G = DEXP(X)*G |
|---|
| 334 | RETURN |
|---|
| 335 | 160 G = A*G*DEXP(-ALGP1) |
|---|
| 336 | GSTAR = 1.D0 - G |
|---|
| 337 | RETURN |
|---|
| 338 | !C |
|---|
| 339 | !C RECURSIVE EVALUATION OF G(A,X) FOR X.LE.1.5 AND A.LT.-.5 |
|---|
| 340 | !C |
|---|
| 341 | 170 I = 1 |
|---|
| 342 | AE = AEPS |
|---|
| 343 | AEP1 = AEPS + 1.D0 |
|---|
| 344 | IF (X.LT..25D0 .AND. AE.GT.ALPHA) GO TO 210 |
|---|
| 345 | GO TO 70 |
|---|
| 346 | 180 G = G*DEXP(X)*X**(-AE) |
|---|
| 347 | DO 190 K=1,MA |
|---|
| 348 | G = (1.D0-X*G)/(DBLE(FLOAT(K))-AE) |
|---|
| 349 | 190 CONTINUE |
|---|
| 350 | ALG = DLOG(G) |
|---|
| 351 | !C |
|---|
| 352 | !C EVALUATION OF GSTAR(A,X) IN TERMS OF G(A,X) |
|---|
| 353 | !C |
|---|
| 354 | 200 GSTAR = 1.D0 |
|---|
| 355 | IF (MA.GE.0 .AND. AEPS.EQ.0.D0) RETURN |
|---|
| 356 | SGT = SGA*SGGA |
|---|
| 357 | T = DLOG(AA) - X + A*ALX + ALG - ALGP1 |
|---|
| 358 | IF (T.LT.-ALPREC) RETURN |
|---|
| 359 | IF (T.GE.TOP) GO TO 320 |
|---|
| 360 | GSTAR = 1.D0 - SGT*DEXP(T) |
|---|
| 361 | RETURN |
|---|
| 362 | 210 I = 2 |
|---|
| 363 | ALGEP1 = DLGA(AEP1) |
|---|
| 364 | !C |
|---|
| 365 | !C EVALUATION OF GSTAR(A,X) FOR A.GT.ALPHA(X) BY TAYLOR |
|---|
| 366 | !C EXPANSION |
|---|
| 367 | !C |
|---|
| 368 | 220 G = 1.D0 |
|---|
| 369 | TERM = 1.D0 |
|---|
| 370 | SUMM = 1.D0 |
|---|
| 371 | K = 0 |
|---|
| 372 | 230 K = K + 1 |
|---|
| 373 | IF (K.GT.600) GO TO 340 |
|---|
| 374 | TERM = X*TERM/(AE+DBLE(FLOAT(K))) |
|---|
| 375 | SUMM = SUMM + TERM |
|---|
| 376 | IF (DABS(TERM).GT.EPS*SUMM) GO TO 230 |
|---|
| 377 | ALGS = AE*ALX - X + DLOG(SUMM) - ALGEP1 |
|---|
| 378 | IF (ALGS.LE.BOT) GO TO 310 |
|---|
| 379 | GSTAR = DEXP(ALGS) |
|---|
| 380 | G = 1.D0 - GSTAR |
|---|
| 381 | IF (I.NE.2) RETURN |
|---|
| 382 | G = G*DEXP(ALGEP1)/AE |
|---|
| 383 | GO TO 180 |
|---|
| 384 | !C |
|---|
| 385 | !C EVALUATION OF G(A,X) FOR X.GT.1.5 AND A.LE.ALPHA(X) BY |
|---|
| 386 | !C MEANS OF THE LEGENDRE CONTINUED FRACTION |
|---|
| 387 | !C |
|---|
| 388 | 240 GSTAR = 1.D0 |
|---|
| 389 | XPA = X + 1.D0 - A |
|---|
| 390 | XMA = X - 1.D0 - A |
|---|
| 391 | P = 0.D0 |
|---|
| 392 | Q = XPA*XMA |
|---|
| 393 | R = 4.D0*XPA |
|---|
| 394 | S = -A + 1.D0 |
|---|
| 395 | TERM = 1.D0 |
|---|
| 396 | SUMM = 1.D0 |
|---|
| 397 | RHO = 0.D0 |
|---|
| 398 | K = 1 |
|---|
| 399 | 250 K = K + 1 |
|---|
| 400 | IF (K.GT.600) GO TO 330 |
|---|
| 401 | P = P + S |
|---|
| 402 | Q = Q + R |
|---|
| 403 | R = R + 8.D0 |
|---|
| 404 | S = S + 2.D0 |
|---|
| 405 | T = P*(1.D0+RHO) |
|---|
| 406 | RHO = T/(Q-T) |
|---|
| 407 | TERM = RHO*TERM |
|---|
| 408 | SUMM = SUMM + TERM |
|---|
| 409 | IF (DABS(TERM).GT.EPS*SUMM) GO TO 250 |
|---|
| 410 | IF (A) 260, 270, 280 |
|---|
| 411 | 260 G = SUMM/XPA |
|---|
| 412 | ALG = DLOG(G) |
|---|
| 413 | GO TO 200 |
|---|
| 414 | 270 G = SUMM/XPA |
|---|
| 415 | RETURN |
|---|
| 416 | 280 ALG = A*ALX - X + DLOG(A*SUMM/XPA) - ALGP1 |
|---|
| 417 | IF (ALG.LE.BOT) GO TO 300 |
|---|
| 418 | G = DEXP(ALG) |
|---|
| 419 | GSTAR = 1.D0 - G |
|---|
| 420 | RETURN |
|---|
| 421 | 290 IFLG = 1 |
|---|
| 422 | IFLGST = 1 |
|---|
| 423 | RETURN |
|---|
| 424 | 300 IFLG = 4 |
|---|
| 425 | RETURN |
|---|
| 426 | 310 IFLGST = 4 |
|---|
| 427 | RETURN |
|---|
| 428 | 320 IFLGST = 5 |
|---|
| 429 | GSTAR = -SGT*AINF |
|---|
| 430 | RETURN |
|---|
| 431 | 330 IFLG = 6 |
|---|
| 432 | RETURN |
|---|
| 433 | 340 IFLGST = 6 |
|---|
| 434 | RETURN |
|---|
| 435 | END |
|---|
| 436 | |
|---|
| 437 | DOUBLE PRECISION FUNCTION DLGA(DX) |
|---|
| 438 | DOUBLE PRECISION DBNUM, DBDEN, DX, DC, DP, DY, DT, DS |
|---|
| 439 | DIMENSION DBNUM(8), DBDEN(8) |
|---|
| 440 | DATA DBNUM /-3.617D3,1.D0,-6.91D2,1.D0,-1.D0,1.D0,-1.D0,1.D0/, & |
|---|
| 441 | & DBDEN /1.224D5,1.56D2,3.6036D5,1.188D3,1.68D3,1.26D3,3.6D2,1.2D1/ |
|---|
| 442 | DC = .5D0*DLOG(8.D0*DATAN(1.D0)) |
|---|
| 443 | DP = 1.D0 |
|---|
| 444 | DY = DX |
|---|
| 445 | Y = SNGL(DY) |
|---|
| 446 | !C |
|---|
| 447 | !C THE CONDITIONAL CLAUSE IN THE NEXT STATEMENT EXPRESSES THE |
|---|
| 448 | !C INEQUALITY Y.GT.EXP(.121189*DPREC+.053905), WHERE DPREC IS THE |
|---|
| 449 | !C NUMBER OF DECIMAL DIGITS CARRIED IN DOUBLE PRECISION FLOATING-POINT |
|---|
| 450 | !C ARITHMETIC. |
|---|
| 451 | !C |
|---|
| 452 | 10 IF (Y.GT.35.027) GO TO 20 |
|---|
| 453 | DP = DY*DP |
|---|
| 454 | DY = DY + 1.D0 |
|---|
| 455 | Y = SNGL(DY) |
|---|
| 456 | GO TO 10 |
|---|
| 457 | 20 DT = 1.D0/(DY*DY) |
|---|
| 458 | DS = 4.3867D4/2.44188D5 |
|---|
| 459 | DO 30 I=1,8 |
|---|
| 460 | DS = DT*DS + DBNUM(I)/DBDEN(I) |
|---|
| 461 | 30 CONTINUE |
|---|
| 462 | DLGA = (DY-.5D0)*DLOG(DY) - DY + DC + DS/DY - DLOG(DP) |
|---|
| 463 | RETURN |
|---|
| 464 | END |
|---|
| 465 | |
|---|
| 466 | |
|---|
