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