1#if defined(HPUX) 2 double precision function derf(x) 3 implicit none 4 double precision x,util_erf,derf 5 external util_erf 6 derf=util_erf(x) 7 return 8 end 9 10 double precision function derfc(x) 11 implicit none 12 double precision x,util_erfc,derfc 13 external util_erfc 14 derfc=util_erfc(x) 15 return 16 end 17c 18 double precision function erfc(x) 19 implicit none 20 double precision x,util_erfc 21 external util_erfc 22 erfc=util_erfc(x) 23 return 24 end 25c 26#endif 27 28 real*8 function util_erf(x) 29 implicit none 30#if NO_ERF 31 real*8 x 32* **** local variables **** 33 real*8 f 34* ***** external functions **** 35 real*8 gammp 36 external gammp 37 IF(x.lt.0.0d0)then 38 f = -gammp(0.5d0,x**2) 39 ELSE 40 f = +gammp(0.5d0,x**2) 41 ENDIF 42 util_erf = f 43#else 44 double precision x 45 double precision derf 46 util_erf=derf(x) 47#endif 48 return 49 end 50 51 real*8 function util_erfc(x) 52 implicit none 53#ifdef NO_ERF 54 real*8 x 55* **** local variables **** 56 real*8 f 57* ***** external functions **** 58 real*8 gammp 59 external gammp 60 IF(x.lt.0.0d0)then 61 f = -gammp(0.5d0,x**2) 62 ELSE 63 f = +gammp(0.5d0,x**2) 64 ENDIF 65 util_erfc = 1.0d0-f 66#else 67 double precision x 68 double precision derfc 69 util_erfc=derfc(x) 70#endif 71 return 72 end 73 74c************************************************************************* 75c ! 76c ! This is a slightly modified version of Log(Gamma) function 77c ! program from Num. Rec. 78c !************************************************************************ 79 80 real*8 function ln_gamma(xx) 81 implicit none 82 real*8 xx 83 84* *** local variables **** 85 integer j 86 real*8 ser,stp,tmp,x,y,cof(6) 87 SAVE cof,stp 88 DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, 89 >24.01409824083091d0,-1.231739572450155d0,0.1208650973866179d-2, 90 >-.5395239384953d-5,2.5066282746310005d0/ 91 x=xx 92 y=x 93 tmp=x+5.5d0 94 tmp=(x+0.5d0)*dlog(tmp)-tmp 95 ser=1.000000000190015d0 96 DO 11 j=1,6 97 y=y+1.0d0 98 ser=ser+cof(j)/y 9911 CONTINUE 100 101 ln_gamma=tmp+dlog(stp*ser/x) 102 103 return 104 END 105 106 real*8 function util_gamma(x) 107 implicit none 108 real*8 x 109 110 real*8 xx 111 real*8 ln_gamma 112 external ln_gamma 113 114 XX = X 115 util_gamma = dexp(ln_gamma(xx)) 116 117 return 118 END 119 120 real*8 function util_gammp(a,x) 121 implicit none 122 real*8 a,x 123 124#include "errquit.fh" 125 126* **** external functions **** 127 real*8 gammcf,gamser,gln 128 129 IF(x.LT.0.0d0 .OR. a.LE.0.0d0) THEN 130 call errquit('bad arguments in util_gammp',0, INPUT_ERR) 131 END IF 132 133 IF (x .lt. (a+1.0d0)) THEN 134 call gser(gamser,a,x,gln) 135 util_gammp=gamser 136 ELSE 137 CALL gcf(gammcf,a,x,gln) 138 util_gammp=1.0d0-gammcf 139 ENDIF 140 return 141 end 142 143 144 145 real*8 function gammp(a,x) 146 implicit none 147 real*8 a,x 148 149#include "errquit.fh" 150 151* **** external functions **** 152 real*8 gammcf,gamser,gln 153 154 IF(x.LT.0.0d0 .OR. a.LE.0.0d0) THEN 155 call errquit('bad arguments in gammp',0, INPUT_ERR) 156 END IF 157 158 IF (x .lt. (a+1.0d0)) THEN 159 call gser(gamser,a,x,gln) 160 gammp=gamser 161 ELSE 162 CALL gcf(gammcf,a,x,gln) 163 gammp=1.0d0-gammcf 164 ENDIF 165 return 166 end 167 168 SUBROUTINE gcf(gammcf,a,x,gln) 169 implicit none 170#include "errquit.fh" 171 integer ITMAX 172 real*8 a,gammcf,gln,x,EPS,FPMIN 173 PARAMETER (ITMAX=1000,EPS=3.0d-12,FPMIN=1.0d-30) 174 real*8 an,b,c,d,del,h 175 double precision undovl 176 parameter(undovl=-20d0*2.3025d0) 177 integer i 178 real*8 ln_gamma 179 external ln_gamma 180 181 gln=ln_gamma(a) 182 b=x + 1.0d0 - a 183 c=1.0d0/FPMIN 184 d=1.0d0/b 185 h=d 186 DO 11 i=1,ITMAX 187 an=-i*(i-a) 188 b=b+2.0d0 189 d=an*d+b 190 IF(DABS(d).lt.FPMIN) d=FPMIN 191 c=b+an/c 192 IF(DABS(c).lt.FPMIN) c=FPMIN 193 d=1.0d0/d 194 del=d*c 195 h=h*del 196 IF(DABS(del-1.0d0).lt.EPS) GOTO 1 19711 CONTINUE 198 call errquit('a too large, ITMAX too small in gcf',0, INPUT_ERR) 1991 if((-x+a*LOG(x)-gln).gt.undovl) then 200 gammcf=DEXP(-x+a*LOG(x)-gln)*h 201 else 202 gammcf=0d0 203 endif 204 return 205 END 206 207 208 209 subroutine gser(gamser,a,x,gln) 210 implicit none 211#include "errquit.fh" 212 real*8 a,x 213 real*8 gamser,gln 214 215* **** parameters **** 216 integer ITMAX 217 parameter (ITMAX=1000) 218 real*8 EPS 219 parameter (EPS=3.0d-12) 220 221* **** local variables **** 222 integer n 223 real*8 ap,del,sum 224 225* **** external functions *** 226 real*8 ln_gamma 227 external ln_gamma 228 229 230 gln=ln_gamma(a) 231 232 IF(x.le.0.0d0) THEN 233 IF(x.lt.0.0d0) call errquit('x < 0 in gser',0, UNKNOWN_ERR) 234 gamser=0.0d0 235 RETURN 236 ENDIF 237 238 ap=a 239 sum=1.0d0/a 240 del=sum 241 DO n=1,ITMAX 242 ap=ap+1.0d0 243 del=del*x/ap 244 sum=sum+del 245 IF(DABS(del).lt.DABS(sum)*EPS) GOTO 1 246 END DO 247 248 CALL errquit('a too large, ITMAX too small in gser',0, INPUT_ERR) 249 2501 gamser=sum*DEXP(-x+a*LOG(x)-gln) 251 252 return 253 END 254c $Id$ 255