44 real :: x,tmp,ser,xx,
gammln
45 real(KIND=dp) :: cof(6) = (/ &
46 76.18009173_dp, -86.50532033_dp, 24.01409822_dp, &
47 -1.231739516_dp, .120858003e-2_dp, -.536382e-5_dp /)
48 real(KIND=dp) :: stp = 2.50662827465_dp
49 real(KIND=dp) :: half = 0.5_dp, one = 1.0_dp, fpf = 5.5_dp
53 tmp=(x+half)*log(tmp)-tmp
66 real :: a, x, gln, gamser,
gammp, gammcf
68 if(x .lt. 0. .or. a .le. 0.)
then
73 call
gser(gamser,a,x,gln)
76 call
gcf(gammcf,a,x,gln)
85 real :: a, x, gln, gamser,
gammq, gammcf
87 if(x.lt.0..or.a.le.0.)
then
92 call
gser(gamser,a,x,gln)
95 call
gcf(gammcf,a,x,gln)
100 subroutine gser(gamser,a,x,gln)
105 real :: gamser, a, x, gln, ap, summ, del
108 integer,
parameter :: itmax=100
109 real,
parameter :: eps=3.e-7
127 if(abs(del).lt.abs(summ)*eps)go to 1
129 print*,
'gser: a too large, itmax too small'
131 1 gamser=summ*exp(-x+a*log(x)-gln)
134 subroutine gcf(gammcf,a,x,gln)
139 real :: gammcf, a, x, gln, gold, a0, a1, b0, b1, fac, an, anf, ana, g
142 integer,
parameter :: itmax=100
143 real,
parameter :: eps=3.e-7
163 if(abs((g-gold)/g).lt.eps)go to 1
167 print*,
'gcf: a too large, itmax too small'
169 1 gammcf=exp(-x+a*alog(x)-gln)*g
177 real,
external ::
gammp
204 real :: x, z, t,
erfcc
208 erfcc=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+ &
209 t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+ &
210 t*(1.48851587+t*(-.82215223+t*.17087277)))))))))
real function gammp(a, x)
real function gammq(a, x)
subroutine gcf(gammcf, a, x, gln)
subroutine gser(gamser, a, x, gln)