PageRenderTime 4391ms CodeModel.GetById 29ms RepoModel.GetById 0ms app.codeStats 0ms

/DOC/Derek Nielson Thesis/Research Documents/UAFWLIS-V1_01_FiniteSizeGround/OMEGA_complexZ/Che0.f

http://ua-fwlis.googlecode.com/
FORTRAN Legacy | 317 lines | 153 code | 31 blank | 133 comment | 0 complexity | 6fa1dcb9abd3db95817bb83142b20b6e MD5 | raw file
Possible License(s): GPL-3.0
  1. SUBROUTINE CHE0(A,Z,NUMA,MNUMA,PCKEXP1,PCKEXP2,SQASP11,SQASP12,
  2. & EXPAZ,H01,H02,H11,H12,HE01CX,HE02CX,HE01EX,
  3. & HE02EX,HE01SX,HE02SX,HE01D0,HE02D0,FLAG,RELERR)
  4. * *
  5. *******************************************************************************
  6. * REFERENCES: *
  7. * *
  8. * [2] S. L. Dvorak and E. F. Kuester, ``Numerical Computation of the *
  9. * Incomplete Lipschitz-Hankel Integral Je0(a,z),'' J. Comput. Phys., *
  10. * Vol. 87, No. 2, pp. 301--327, April 1990. *
  11. * [3] M. Abramowitz and I. E. Stegun, Handbook of Mathematical *
  12. * Functions with Formulas, Graphs, and Mathematical Tables, *
  13. * Washington, D. C.: U.S. Government Printing Office, 1972. *
  14. * [4] David L. Heckmann and S.L. Dvorak,"Numerical Computation of Hankel *
  15. * Functions of Integral Order for Complex-valued Arguments *
  16. * [6] S. L. Dvorak, ``Numerical Computation of the Incomplete *
  17. * Lipschitz-Hankel Integral Ye0(a,z),'' submitted to the SIAM *
  18. * Journal on Mathematical Analysis. *
  19. * [7] Z.ZHU, S.L.DVORAK AND DAVID HECKMAN,"Numerical Computation of *
  20. * Incomplete Lipschitz-Hankel Integrals of the Hankel Type for *
  21. * Complex-Valued Arguments *
  22. * *
  23. *******************************************************************************
  24. * PURPOSE: *
  25. * *
  26. * This subroutine determines which expansion will provide the most *
  27. * efficient computation of the incomplete Lipschitz-Hankel integrals *
  28. * AND OUTPUT THE NUMERICAL VALUE FOR: *
  29. * HE0*exp(az), CHEO*exp(az) and SCHE0*exp(az) *
  30. * HE0 means standard ILHI of Hankel type. *
  31. * CHE0 means complimentary ILHI of Hankel type *
  32. * SCHE0 means CHEO(a,*,z) *
  33. * *
  34. * where A and Z are complex numbers. and Real(Z)>=0
  35. * The integrals are evaluated to SD significant digits in most cases. *
  36. * The relative error is given by RELERR = .5*10.**(-SD). *
  37. *******************************************************************************
  38. * INPUTS: *
  39. * *
  40. * A = Factor in the exponential of the incomplete Lipschitz-Hankel *
  41. * integral. (Complex array, Size is set by MNUMA) *
  42. * Z = Upper limit of integration. (Complex, Real(Z)>0) *
  43. * NUMA = Number of actual terms in the A array. *
  44. * FLAG = 0; NEED TO COMPUTE BOTH HE01 AND HE02 *
  45. * = 1; ONLY NEED TO COMPUTE HE01 *
  46. * = 2; ONLY NEED TO COMPUTE HE02 *
  47. *************************************************************************************************
  48. * OUTPUTS: *
  49. * *
  50. * PCKEXP1/PCKEXP2 = Tells which expansion was used. *
  51. * HE01EX/HE02EX = Computed value for the incomplete Lipschitz-Hankel integral. *
  52. * HE01CX/HE02CX = Computed value for the complimentary incomplete Lipschitz-Hankel integral.*
  53. * HE01SX/HE02SX = Computed value for the incomplete Lipschitz-Hankel star integral. *
  54. * HE01D0/HE02D0 = RETURN HE01(A,DELTA,0)*EXPAZ/HE02(A,DELTA,0)*EXPAZ IF COMPUTED. *
  55. * SQASP11/SQASP12= SQRT(a^2+1) CONSIDERING THE BRANCH CUT WHEN USING ASYMPTOTIC OR *
  56. * QUASI-NEUMANN.OTHERWISE, ALWAYS POSITVE REAL PART. *
  57. *************************************************************************************************
  58. * SUBROUTINES AND FUNCTIONS CALLED: *
  59. * *
  60. * CHANKL = Computes Jn, Yn AND HN1, HN2 *
  61. * CFNDEXP1/CFNDEXP2 = Determine which expansion to use. *
  62. * CFDHE01/CFDHE02 = Compute HE01/HE02 using the proper expansion. *
  63. *******************************************************************************
  64. * VARIABLE DECLARATIONS:
  65. *
  66. $Declare
  67. PARAMETER (MAXNOJ=500)
  68. INTEGER,INTENT(IN):: NUMA,MNUMA,FLAG
  69. COMPLEX*16,INTENT(IN)::Z,A(MNUMA)
  70. COMPLEX*16,INTENT(OUT):: HE01D0(MNUMA),HE02D0(MNUMA)
  71. COMPLEX*16,INTENT(INOUT)::SQASP11(MNUMA),EXPAZ(MNUMA),
  72. & SQASP12(MNUMA)
  73. COMPLEX*16,INTENT(OUT):: HE01CX(MNUMA), HE02CX(MNUMA)
  74. COMPLEX*16,INTENT(OUT):: HE01EX(MNUMA), HE02EX(MNUMA)
  75. COMPLEX*16,INTENT(OUT):: HE01SX(MNUMA), HE02SX(MNUMA)
  76. COMPLEX*16,INTENT(OUT):: H01, H02, H11, H12
  77. INTEGER,INTENT(OUT):: PCKEXP1(MNUMA),PCKEXP2(MNUMA)
  78. DOUBLE PRECISION,INTENT(IN):: RELERR
  79. *********LOCAL VARIABLE************************************************
  80. DOUBLE PRECISION ZJASY,ZSP,COMPAR,ABASP1
  81. COMPLEX*16 JBESS(0:MAXNOJ),YBESS(0:MAXNOJ)
  82. COMPLEX*16 JBESSAZ(0:MAXNOJ),YBESSAZ(0:MAXNOJ)
  83. COMPLEX*16 HN1(0:MAXNOJ), HN2(0:MAXNOJ)
  84. COMPLEX*16 BFH0,BFH1,HN1AZ(0:MAXNOJ), HN2AZ(0:MAXNOJ)
  85. COMPLEX*16 ASP1,U(0:MAXNOJ),LOGA(MNUMA),LOGNA(MNUMA),FRONT(MNUMA)
  86. COMPLEX*16 EK(0:MAXNOJ),P(0:MAXNOJ), HE01, HE02,CE1
  87. INTEGER NOJ,KMAX_1(NUMA),INTZ,NOH1,NOH2,N,FWD_FLAG,NOY
  88. INTEGER NOJAZ,NOYAZ,NOH1AZ,NOH2AZ,KMAX_2(NUMA),K,FWD_FLAGAZ
  89. DOUBLE PRECISION APROX1, APROX2,AZ
  90. DOUBLE PRECISION IMZ
  91. LOGICAL NEWZ, NEEDJ,NEEDSTRU,NEEDAZ,NEEDUEP,NEEDJ_FLAG
  92. INTEGER FWD_FLAG1,FWD_FLAG2
  93. PARAMETER (E=2.71828)
  94. ****************************************************************************
  95. * NEEDJ_FLAG=0 ALL THE SERIES IS SECOND CONVERGENT, SO NO NEED TO COMPUTE Jn(Z)
  96. * =1 AT LEAST ONE OF THE SERIES NEED COMPUTE Jn(Z) BY BACKWARD
  97. * =2 ALL THE SERIES CAN USE FORWARD TO GET HN1 AND HN2, NO Jn(Z) NEEDED
  98. *******************************************************************************
  99. * COMMON BLOCKS
  100. *
  101. * FACTOR(K)=SQPI/(2.**(K+1)*GAMMA(K+1.5))
  102. DOUBLE PRECISION FACTOR(0:500)
  103. COMMON/FACT/FACTOR
  104. DOUBLE PRECISION PI, TWOPI
  105. COMPLEX*16 IM
  106. COMMON/CONST/ PI, TWOPI,IM
  107. INTEGER ZLO, ALO(2)
  108. COMMON /ZALO/ ZLO, ALO
  109. INTEGER SD
  110. COMMON /SIGDIG/ SD
  111. *******************************************************************************
  112. EXTERNAL CHANKL,CFDHEXP1,CFDHE01,CFDHEXP2,CFDHE02,CHAE01,CHAE02
  113. ****************************************************************************
  114. * Calculate an array of Hankel functions. This subroutine calculates
  115. * as many Hankel functions that are required for forward recurrence
  116. * to be stable for higher orders. The higher order Hankel functions
  117. * are calculated as needed.[4]
  118. NEEDJ=.TRUE.
  119. DO 10 N = 1, NUMA
  120. *SQASP1 is always choose the positive real part.
  121. SQASP12(N)=SQASP11(N)
  122. 10 CONTINUE
  123. ZSP=ABS(Z)
  124. INTZ=NINT(ZSP+2)
  125. ZJASY = REAL(SD+5)
  126. **************** Determine which expansion to use.
  127. IF(FLAG==2)THEN
  128. DO N=1,NUMA
  129. CALL CFDHEXP2(A(N),Z,SQASP12(N),RELERR,KMAX_2(N),PCKEXP2(N))
  130. END DO
  131. * Determine if the CHANKL.F need to be called, i.e. if need to compute Jn(Z) by backward
  132. * THOUGH PCKEXP=12 NEED NOT COMPUTE HN1/2(Z), HOWEVER, WE NEED RETURN HN1/2(Z) FIRST AND SECOND ORDER
  133. * SET DEFAULT TO USE ASYMPTOTIC EXPRESION TO COMPUTE HN1/2(0) AND HN1/2(1)
  134. NEEDJ_FLAG=.FALSE.
  135. DO N=1,NUMA
  136. IF (PCKEXP2(N)>=14.OR. ZSP<20.0) THEN
  137. * USE BACKWARD TO COMPUTE JBESS, HN1/2(0) AND HN1/2(1)
  138. NEEDJ_FLAG=.TRUE.
  139. END IF
  140. END DO
  141. IF(NEEDJ_FLAG==.FALSE.) THEN
  142. * USE ASYMPTOTIC EXPRESION TO COMPUTE HN1/2(0) AND HN1/2(1)
  143. CALL CHAE02(Z,RELERR*1.D-6,HN2(0),HN2(1))
  144. IF (REAL(Z)==0. .AND. IMAG(Z)<0.) THEN
  145. HN2(0)=IM*IMAG(HN2(0))
  146. HN2(1)=DBLE(HN2(1))
  147. END IF
  148. NOH2=1
  149. FWD_FLAG2=2
  150. ELSE
  151. * USE BEWARD TO COMPUTE Jn AND HN1/HN2
  152. CALL CHANKL(Z,RELERR,JBESS,NEEDJ,MAXNOJ,NOJ,HN1,HN2,
  153. & NOH1,NOH2,YBESS,NOY,FWD_FLAG)
  154. FWD_FLAG2=FWD_FLAG
  155. END IF
  156. ***********Compute the ILHIs using the proper expansion.
  157. NEEDSTRU=.TRUE.
  158. DO N = 1, NUMA
  159. NEEDAZ=.TRUE.
  160. NEEDUEP=.TRUE.
  161. IF(PCKEXP2(N)==14) THEN
  162. DO K = 0, MAXNOJ
  163. U(K) = 0.D0
  164. END DO
  165. END IF
  166. CALL CFDHE02(Z,A(N),RELERR,SQASP12(N),PCKEXP2(N),KMAX_2(N),
  167. & EXPAZ(N),JBESS,YBESS,MAXNOJ,NOJ,NOY,HN1,HN2,NOH1,NOH2,
  168. & NEEDJ,HE02CX(N),HE02EX(N),HE02SX(N),HE02D0(N),U,
  169. & EK,P,NEEDUEP,ZJASY,FWD_FLAG2,BFH0,BFH1,NEEDSTRU,JBESSAZ,
  170. & YBESSAZ,NOJAZ,NOYAZ,HN1AZ,HN2AZ,NOH1AZ,NOH2AZ,NEEDAZ,
  171. & FWD_FLAGAZ,CE1)
  172. END DO
  173. ELSE IF(FLAG==1)THEN
  174. DO N=1,NUMA
  175. CALL CFDHEXP1(A(N),Z,SQASP11(N),RELERR,KMAX_1(N),PCKEXP1(N))
  176. END DO
  177. * Determine if the CHANKL.F need to be called, i.e. if need to compute Jn(Z) by backward
  178. * THOUGH PCKEXP=12 NEED NOT COMPUTE HN1/2(Z), HOWEVER, WE NEED RETURN HN1/2(Z) FIRST AND SECOND ORDER
  179. * SET DEFAULT TO USE ASYMPTOTIC EXPRESION TO COMPUTE HN1/2(0) AND HN1/2(1)
  180. NEEDJ_FLAG=.FALSE.
  181. DO N=1,NUMA
  182. IF (PCKEXP1(N)>=14.OR. ZSP<20.0) THEN
  183. * USE BACKWARD TO COMPUTE JBESS, HN1/2(0) AND HN1/2(1)
  184. NEEDJ_FLAG=.TRUE.
  185. END IF
  186. END DO
  187. IF(NEEDJ_FLAG==.FALSE.) THEN
  188. * USE ASYMPTOTIC EXPRESION TO COMPUTE HN1/2(0) AND HN1/2(1)
  189. CALL CHAE01(Z,RELERR*1.D-6,HN1(0),HN1(1))
  190. IF(REAL(Z)==0. .AND. IMAG(Z)>0.) THEN
  191. HN1(0)=IM*IMAG(HN1(0))
  192. HN1(1)=DBLE(HN1(1))
  193. END IF
  194. NOH1=1
  195. FWD_FLAG1=1
  196. ELSE
  197. * USE BEWARD TO COMPUTE Jn AND HN1/HN2
  198. CALL CHANKL(Z,RELERR,JBESS,NEEDJ,MAXNOJ,NOJ,HN1,HN2,
  199. & NOH1,NOH2,YBESS,NOY,FWD_FLAG)
  200. FWD_FLAG1=FWD_FLAG
  201. END IF
  202. ***********Compute the ILHIs using the proper expansion.
  203. NEEDSTRU=.TRUE.
  204. DO N = 1, NUMA
  205. NEEDAZ=.TRUE.
  206. NEEDUEP=.TRUE.
  207. IF(PCKEXP1(N)==14) THEN
  208. DO K = 0, MAXNOJ
  209. U(K) = 0.D0
  210. END DO
  211. END IF
  212. CALL CFDHE01(Z,A(N),RELERR,SQASP11(N),PCKEXP1(N),KMAX_1(N),
  213. & EXPAZ(N),JBESS,YBESS,MAXNOJ,NOJ,NOY,HN1,HN2,NOH1,NOH2,
  214. & NEEDJ,HE01CX(N),HE01EX(N),HE01SX(N),HE01D0(N),U,
  215. & EK,P,NEEDUEP,ZJASY,FWD_FLAG1,BFH0,BFH1,NEEDSTRU,JBESSAZ,
  216. & YBESSAZ,NOJAZ,NOYAZ,HN1AZ,HN2AZ,NOH1AZ,NOH2AZ,NEEDAZ,
  217. & FWD_FLAGAZ,CE1)
  218. END DO
  219. ELSE
  220. DO N=1,NUMA
  221. CALL CFDHEXP1(A(N),Z,SQASP11(N),RELERR,KMAX_1(N),PCKEXP1(N))
  222. CALL CFDHEXP2(A(N),Z,SQASP12(N),RELERR,KMAX_2(N),PCKEXP2(N))
  223. END DO
  224. * Determine if the CHANKL.F need to be called, i.e. if need to compute Jn(Z) by backward
  225. * THOUGH PCKEXP=12 NEED NOT COMPUTE HN1/2(Z), HOWEVER, WE NEED RETURN HN1/2(Z) FIRST AND SECOND ORDER
  226. * SET DEFAULT TO USE ASYMPTOTIC EXPRESION TO COMPUTE HN1/2(0) AND HN1/2(1)
  227. NEEDJ_FLAG=.FALSE.
  228. DO N=1,NUMA
  229. IF (PCKEXP1(N)>=14 .OR. PCKEXP2(N)>=14
  230. & .OR. ZSP<20.0) THEN
  231. * USE BACKWARD TO COMPUTE JBESS, HN1/2(0) AND HN1/2(1)
  232. NEEDJ_FLAG=.TRUE.
  233. END IF
  234. END DO
  235. IF(NEEDJ_FLAG==.FALSE.) THEN
  236. * USE ASYMPTOTIC EXPRESION TO COMPUTE HN1/2(0) AND HN1/2(1)
  237. CALL CHAE01(Z,RELERR*1.D-6,HN1(0),HN1(1))
  238. CALL CHAE02(Z,RELERR*1.D-6,HN2(0),HN2(1))
  239. IF(REAL(Z)==0. .AND. IMAG(Z)>0.) THEN
  240. HN1(0)=IM*IMAG(HN1(0))
  241. HN1(1)=DBLE(HN1(1))
  242. ELSE IF (REAL(Z)==0. .AND. IMAG(Z)<0.) THEN
  243. HN2(0)=IM*IMAG(HN2(0))
  244. HN2(1)=DBLE(HN2(1))
  245. END IF
  246. NOH1=1
  247. NOH2=1
  248. FWD_FLAG1=1
  249. FWD_FLAG2=2
  250. ELSE
  251. * USE BEWARD TO COMPUTE Jn AND HN1/HN2
  252. CALL CHANKL(Z,RELERR,JBESS,NEEDJ,MAXNOJ,NOJ,HN1,HN2,
  253. & NOH1,NOH2,YBESS,NOY,FWD_FLAG)
  254. FWD_FLAG1=FWD_FLAG
  255. FWD_FLAG2=FWD_FLAG
  256. END IF
  257. ***********Compute the ILHIs using the proper expansion.
  258. NEEDSTRU=.TRUE.
  259. DO N = 1, NUMA
  260. NEEDAZ=.TRUE.
  261. NEEDUEP=.TRUE.
  262. IF(PCKEXP2(N)==14 .OR. PCKEXP1(N)==14) THEN
  263. DO K = 0, MAXNOJ
  264. U(K) = 0.D0
  265. END DO
  266. END IF
  267. CALL CFDHE01(Z,A(N),RELERR,SQASP11(N),PCKEXP1(N),KMAX_1(N),
  268. & EXPAZ(N),JBESS,YBESS,MAXNOJ,NOJ,NOY,HN1,HN2,NOH1,NOH2,
  269. & NEEDJ,HE01CX(N),HE01EX(N),HE01SX(N),HE01D0(N),U,
  270. & EK,P,NEEDUEP,ZJASY,FWD_FLAG1,BFH0,BFH1,NEEDSTRU,JBESSAZ,
  271. & YBESSAZ,NOJAZ,NOYAZ,HN1AZ,HN2AZ,NOH1AZ,NOH2AZ,NEEDAZ,
  272. & FWD_FLAGAZ,CE1)
  273. CALL CFDHE02(Z,A(N),RELERR,SQASP12(N),PCKEXP2(N),KMAX_2(N),
  274. & EXPAZ(N),JBESS,YBESS,MAXNOJ,NOJ,NOY,HN1,HN2,NOH1,NOH2,
  275. & NEEDJ,HE02CX(N),HE02EX(N),HE02SX(N),HE02D0(N),U,
  276. & EK,P,NEEDUEP,ZJASY,FWD_FLAG2,BFH0,BFH1,NEEDSTRU,JBESSAZ,
  277. & YBESSAZ,NOJAZ,NOYAZ,HN1AZ,HN2AZ,NOH1AZ,NOH2AZ,NEEDAZ,
  278. & FWD_FLAGAZ,CE1)
  279. END DO
  280. END IF
  281. H01=HN1(0)
  282. H11=HN1(1)
  283. H02=HN2(0)
  284. H12=HN2(1)
  285. RETURN
  286. END