1! ************************************************************************* 2! 3! This is a slightly modified version of Log(Gamma) function 4! program from Num. Rec. 5! Serves as backup if intrinic Gamma function 6! is not available 7! ************************************************************************ 8 9 DOUBLE PRECISION FUNCTION PAW_LN_GAMMA (XX) 10 implicit none 11 DOUBLE PRECISION XX 12 INTEGER J 13 DOUBLE PRECISION SER 14 DOUBLE PRECISION STP 15 DOUBLE PRECISION TMP 16 DOUBLE PRECISION X 17 DOUBLE PRECISION Y, COF(6) 18 SAVE STP, COF 19 DATA COF, STP/ 76.18009172947146D0, -86.50532032941677D0, 20 1 24.01409824083091D0, -1.231739572450155D0, 21 2 0.1208650973866179D-2, -.5395239384953D-5, 2.5066282746310005D0 22 3 / 23 X = XX 24 Y = X 25 TMP = X + 5.5D0 26 TMP = (X + 0.5D0)*DLOG(TMP) - TMP 27 SER = 1.000000000190015D0 28 DO J = 1, 6 29 Y = Y + 1.0D0 30 SER = SER + COF(J)/Y 31 END DO 32 33 PAW_LN_GAMMA = TMP + DLOG(STP*SER/X) 34 35 RETURN 36 END 37 38! ************************************************* 39! 40! Name : 41! 42! 43! Purpose : 44! 45! 46! Created : 47! 48! ************************************************* 49 DOUBLE PRECISION FUNCTION PAW_GAMMA(X) 50 implicit none 51 DOUBLE PRECISION X 52 53 DOUBLE PRECISION XX 54 DOUBLE PRECISION PAW_LN_GAMMA 55 56 XX = X 57 PAW_GAMMA = DEXP(PAW_LN_GAMMA(XX)) 58 59 return 60 END 61 62 63 64! ************************************************* 65! 66! Name : 67! 68! 69! Purpose : 70! 71! 72! Created : 73! 74! ************************************************* 75 DOUBLE PRECISION FUNCTION PAW_GAMMP (A, X) 76 implicit none 77 DOUBLE PRECISION A, X 78 79 DOUBLE PRECISION GAMMCF, GAMSER, GLN 80 81 IF (X .LT. A+1.0D0) THEN 82 83 CALL PAW_GSER(GAMSER, A, X, GLN) 84 PAW_GAMMP = GAMSER 85 86 ELSE 87 88 CALL PAW_GCF(GAMMCF, A, X, GLN) 89 PAW_GAMMP = 1.0D0 - GAMMCF 90 91 ENDIF 92 93 return 94 END 95 96! ************************************************* 97! 98! Name : 99! 100! 101! Purpose : 102! 103! 104! Created : 105! 106! ************************************************* 107 SUBROUTINE PAW_GCF(GAMMCF, A, X, GLN) 108 implicit none 109 INTEGER ITMAX 110 DOUBLE PRECISION A, GAMMCF, GLN, X, EPS, FPMIN 111 PARAMETER (ITMAX = 100, EPS = 3.D-16, FPMIN = 1.D-30) 112 DOUBLE PRECISION AN, B, C, D, DEL, H 113 INTEGER I 114 115 DOUBLE PRECISION PAW_LN_GAMMA 116 EXTERNAL PAW_LN_GAMMA 117 118 GLN = PAW_LN_GAMMA(A) 119 B = X + 1.0D0 - A 120 C = 1.0D0/1.D-30 121 D = 1.0D0/B 122 H = D 123 DO I = 1, 100 124 AN = -I*(I - A) 125 B = B + 2.0D0 126 D = AN*D + B 127 IF (DABS(D) .LT. 1.D-30) D = 1.D-30 128 C = B + AN/C 129 IF (DABS(C) .LT. 1.D-30) C = 1.D-30 130 D = 1.0D0/D 131 DEL = D*C 132 H = H*DEL 133 IF (DABS(DEL - 1.0D0) .LT. 3.D-16) GO TO 1 134 END DO 135 PAUSE 'a too large, ITMAX too small in gcf' 136 1 CONTINUE 137 GAMMCF = DEXP((-X) + A*DLOG(X) - GLN)*H 138 139 return 140 END 141 142! ************************************************* 143! 144! Name : 145! 146! 147! Purpose : 148! 149! 150! Created : 151! 152! ************************************************* 153 SUBROUTINE PAW_GSER(GAMSER, A, X, GLN) 154 implicit none 155 DOUBLE PRECISION A, X 156 DOUBLE PRECISION GAMSER, GLN 157 158! *** local variables *** 159 INTEGER ITMAX 160 PARAMETER (ITMAX = 100) 161 DOUBLE PRECISION EPS 162 PARAMETER (EPS = 3.0D-16) 163 INTEGER N 164 DOUBLE PRECISION AP, DEL, SUM 165 166 DOUBLE PRECISION PAW_LN_GAMMA 167 EXTERNAL PAW_LN_GAMMA 168 169 170 GLN = PAW_LN_GAMMA(A) 171 172 IF (X .LE. 0.0D0) THEN 173 IF(X.lt.0.0d0) CALL errquit("x < 0 in PAW_GSER",0,1) 174 GAMSER = 0.0D0 175 RETURN 176 ENDIF 177 178 AP = A 179 SUM = 1.0D0/A 180 DEL = SUM 181 DO N = 1, 100 182 AP = AP + 1.0D0 183 DEL = DEL*X/AP 184 SUM = SUM + DEL 185 IF (DABS(DEL) .LT. DABS(SUM)*3.0D-16) GO TO 1 186 187 END DO 188 189 CALL errquit 190 > ("a too large,ITMAX too small in PAW_GSER",0,1) 191 192 1 CONTINUE 193 GAMSER = SUM*DEXP((-X) + A*DLOG(X) - GLN) 194 195 return 196 END 197 198 199 200c $Id$ 201