1      Subroutine igamma(Fj,N,L)
2c $Id$
3
4      Implicit none
5      integer n
6      integer l
7      Double precision Fj(N,0:L)
8
9c Compute the incomplete gamma function.
10c
11c               /1   2j
12c     Fj(T)  =  |   x  exp( - T x**2 ) dx
13c               /0
14c
15c******************************************************************************
16      double precision x,faci
17      integer lmax
18      parameter (lmax=40)
19      double precision ftex(0:lmax)
20      integer m
21      parameter(faci=1d0/1.128379167096d0)
22c
23      if(l.eq.0) then
24        do  m = 1,N
25           call fm(fj(m,0),0,ftex)
26           fj(m,0)=ftex(0)*faci
27        enddo
28      elseif(l.eq.1) then
29        do  m = 1,N
30           call fm(fj(m,0),1,ftex)
31           fj(m,0)=ftex(0)*faci
32           fj(m+N,0)=ftex(1)*faci
33        enddo
34      elseif(l.eq.2) then
35        do  m = 1,N
36           call fm(fj(m,0),2,ftex)
37           fj(m,0)=ftex(0)*faci
38           fj(m+N,0)=ftex(1)*faci
39           fj(m+N+N,0)=ftex(2)*faci
40        enddo
41      elseif(l.eq.3) then
42        do  m = 1,N
43           call fm(fj(m,0),3,ftex)
44           fj(m,0)=ftex(0)*faci
45           fj(m+N,0)=ftex(1)*faci
46           fj(m+N+N,0)=ftex(2)*faci
47           fj(m+N+N+N,0)=ftex(3)*faci
48        enddo
49      else
50        do  m = 1,N
51           call fm(fj(m,0),l,ftex)
52           call dcopy(l+1,ftex,1,fj(m,0),N)
53        enddo
54        call dscal(N*(L+1),faci,fj(1,0),1)
55      endif
56      return
57      end
58      Subroutine igammao(Fj,N,L)
59
60      Implicit real*8 (a-h,o-z)
61      Implicit integer (i-n)
62      Parameter (EXPLIM=100.D0)
63      Parameter (TMAX=30.D0,ACCY=1.D-14)
64*
65      Parameter (PI=3.141592653589793D0)
66      Parameter (PI4=0.25D0*PI)
67
68      Dimension Fj(N,0:L)
69
70c Compute the incomplete gamma function.
71c
72c               /1   2j
73c     Fj(T)  =  |   x  exp( - T x**2 ) dx
74c               /0
75c
76c******************************************************************************
77      do 100 m = 1,N
78
79      Tm = Fj(m,0)
80
81      if( Tm.GT.TMAX )then
82
83       T2inv = 0.5D0/Tm
84       expT = exp(-min(EXPLIM,Tm))
85
86c ... upward recursion.
87
88       Fj(m,0) = sqrt(PI4/Tm)
89       do 10 j = 1,L
90        Fj(m,j) = ((2*j-1)*Fj(m,j-1) - expT)*T2inv
91  10   continue
92
93      else
94
95       T2 = 2.D0*Tm
96       expT = exp(-min(EXPLIM,Tm))
97
98c ... series expansion.
99
100       TERM  = 1 / dble(2*L+1)
101
102       sum  = TERM
103       do 20 i = 1,1000
104        TERM = TERM*T2/dble(2*(L+i)+1)
105        sum = sum + TERM
106        if( TERM.lt.ACCY ) go to 25
107  20   continue
108
109       write(*,*) 'IGAMMA:  Unable to achieve required accuracy.'
110       write(*,*) '         TERM = ',TERM,' ACCY = ',ACCY
111       stop
112
113  25   continue
114
115c ... downward recursion.
116
117       Fj(m,L) = expT*sum
118       do 30 j = L-1,0,-1
119        Fj(m,j) = (T2*Fj(m,j+1) + expT) / dble(2*j+1)
120  30   continue
121
122      end if
123
124  100 continue
125
126*acc_debug:      end
127*acc_debug:      Subroutine igamma_acc(Fj,N,L,accy,met)
128*acc_debug:
129*acc_debug:      Implicit real*8 (a-h,o-z)
130*acc_debug:      Implicit integer (i-n)
131*acc_debug:
132*acc_debug:      Parameter (EXPLIM=100.D0)
133*acc_debug:      Parameter (TMAX=30.D0)
134*acc_debug:*
135*acc_debug:      Parameter (PI=3.141592653589793D0)
136*acc_debug:      Parameter (PI4=0.25D0*PI)
137*acc_debug:
138*acc_debug:      Dimension Fj(N,0:L)
139*acc_debug:      Logical met
140*acc_debug:
141*acc_debug:c Compute the incomplete gamma function.
142*acc_debug:c
143*acc_debug:c               /1   2j
144*acc_debug:c     Fj(T)  =  |   x  exp( - T x**2 ) dx
145*acc_debug:c               /0
146*acc_debug:c
147*acc_debug:c******************************************************************************
148*acc_debug:
149*acc_debug:      met = .true.
150*acc_debug:
151*acc_debug:      do 100 m = 1,N
152*acc_debug:
153*acc_debug:      Tm = Fj(m,0)
154*acc_debug:
155*acc_debug:      if( Tm.GT.TMAX )then
156*acc_debug:
157*acc_debug:       T2inv = 0.5D0/Tm
158*acc_debug:       expT = exp(-min(EXPLIM,Tm))
159*acc_debug:
160*acc_debug:c ... upward recursion.
161*acc_debug:
162*acc_debug:       Fj(m,0) = sqrt(PI4/Tm)
163*acc_debug:       do 10 j = 1,L
164*acc_debug:        Fj(m,j) = ((2*j-1)*Fj(m,j-1) - expT)*T2inv
165*acc_debug:  10   continue
166*acc_debug:
167*acc_debug:      else
168*acc_debug:
169*acc_debug:       T2 = 2.D0*Tm
170*acc_debug:       expT = exp(-min(EXPLIM,Tm))
171*acc_debug:
172*acc_debug:c ... series expansion.
173*acc_debug:
174*acc_debug:       TERM  = 1 / dble(2*L+1)
175*acc_debug:
176*acc_debug:       sum  = TERM
177*acc_debug:       do 20 i = 1,1000
178*acc_debug:        TERM = TERM*T2/dble(2*(L+i)+1)
179*acc_debug:        sum = sum + TERM
180*acc_debug:        if( TERM.lt.ACCY ) go to 25
181*acc_debug:  20   continue
182*acc_debug:
183*acc_debug:       write(*,*) 'IGAMMA_ACC:  Unable to achieve required accuracy.'
184*acc_debug:       write(*,*) '         TERM = ',TERM,' ACCY = ',ACCY
185*acc_debug:       met = .false.
186*acc_debug:       return
187*acc_debug:
188*acc_debug:  25   continue
189*acc_debug:
190*acc_debug:c ... downward recursion.
191*acc_debug:
192*acc_debug:       Fj(m,L) = expT*sum
193*acc_debug:       do 30 j = L-1,0,-1
194*acc_debug:        Fj(m,j) = (T2*Fj(m,j+1) + expT) / dble(2*j+1)
195*acc_debug:  30   continue
196*acc_debug:
197*acc_debug:      end if
198*acc_debug:
199*acc_debug:  100 continue
200*acc_debug:
201      end
202      subroutine igamma_init
203c
204c     call need when txs is not active
205c
206      call fprep
207      return
208      end
209