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