zunk1.f 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426
  1. SUBROUTINE ZUNK1(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM,
  2. * ALIM)
  3. C***BEGIN PROLOGUE ZUNK1
  4. C***REFER TO ZBESK
  5. C
  6. C ZUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
  7. C RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
  8. C UNIFORM ASYMPTOTIC EXPANSION.
  9. C MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
  10. C NZ=-1 MEANS AN OVERFLOW WILL OCCUR
  11. C
  12. C***ROUTINES CALLED ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,ZABS
  13. C***END PROLOGUE ZUNK1
  14. C COMPLEX CFN,CK,CONE,CRSC,CS,CSCL,CSGN,CSPN,CSR,CSS,CWRK,CY,CZERO,
  15. C *C1,C2,PHI,PHID,RZ,SUM,SUMD,S1,S2,Y,Z,ZETA1,ZETA1D,ZETA2,ZETA2D,ZR
  16. DOUBLE PRECISION ALIM, ANG, APHI, ASC, ASCLE, BRY, CKI, CKR,
  17. * CONER, CRSC, CSCL, CSGNI, CSPNI, CSPNR, CSR, CSRR, CSSR,
  18. * CWRKI, CWRKR, CYI, CYR, C1I, C1R, C2I, C2M, C2R, ELIM, FMR, FN,
  19. * FNF, FNU, PHIDI, PHIDR, PHII, PHIR, PI, RAST, RAZR, RS1, RZI,
  20. * RZR, SGN, STI, STR, SUMDI, SUMDR, SUMI, SUMR, S1I, S1R, S2I,
  21. * S2R, TOL, YI, YR, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R,
  22. * ZET1DI, ZET1DR, ZET2DI, ZET2DR, ZI, ZR, ZRI, ZRR, D1MACH, ZABS
  23. INTEGER I, IB, IFLAG, IFN, IL, INIT, INU, IUF, K, KDFLG, KFLAG,
  24. * KK, KODE, MR, N, NW, NZ, INITD, IC, IPARD, J
  25. DIMENSION BRY(3), INIT(2), YR(N), YI(N), SUMR(2), SUMI(2),
  26. * ZETA1R(2), ZETA1I(2), ZETA2R(2), ZETA2I(2), CYR(2), CYI(2),
  27. * CWRKR(16,3), CWRKI(16,3), CSSR(3), CSRR(3), PHIR(2), PHII(2)
  28. DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 /
  29. DATA PI / 3.14159265358979324D0 /
  30. C
  31. KDFLG = 1
  32. NZ = 0
  33. C-----------------------------------------------------------------------
  34. C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
  35. C THE UNDERFLOW LIMIT
  36. C-----------------------------------------------------------------------
  37. CSCL = 1.0D0/TOL
  38. CRSC = TOL
  39. CSSR(1) = CSCL
  40. CSSR(2) = CONER
  41. CSSR(3) = CRSC
  42. CSRR(1) = CRSC
  43. CSRR(2) = CONER
  44. CSRR(3) = CSCL
  45. BRY(1) = 1.0D+3*D1MACH(1)/TOL
  46. BRY(2) = 1.0D0/BRY(1)
  47. BRY(3) = D1MACH(2)
  48. ZRR = ZR
  49. ZRI = ZI
  50. IF (ZR.GE.0.0D0) GO TO 10
  51. ZRR = -ZR
  52. ZRI = -ZI
  53. 10 CONTINUE
  54. J = 2
  55. DO 70 I=1,N
  56. C-----------------------------------------------------------------------
  57. C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
  58. C-----------------------------------------------------------------------
  59. J = 3 - J
  60. FN = FNU + DBLE(FLOAT(I-1))
  61. INIT(J) = 0
  62. CALL ZUNIK(ZRR, ZRI, FN, 2, 0, TOL, INIT(J), PHIR(J), PHII(J),
  63. * ZETA1R(J), ZETA1I(J), ZETA2R(J), ZETA2I(J), SUMR(J), SUMI(J),
  64. * CWRKR(1,J), CWRKI(1,J))
  65. IF (KODE.EQ.1) GO TO 20
  66. STR = ZRR + ZETA2R(J)
  67. STI = ZRI + ZETA2I(J)
  68. RAST = FN/ZABS(COMPLEX(STR,STI))
  69. STR = STR*RAST*RAST
  70. STI = -STI*RAST*RAST
  71. S1R = ZETA1R(J) - STR
  72. S1I = ZETA1I(J) - STI
  73. GO TO 30
  74. 20 CONTINUE
  75. S1R = ZETA1R(J) - ZETA2R(J)
  76. S1I = ZETA1I(J) - ZETA2I(J)
  77. 30 CONTINUE
  78. RS1 = S1R
  79. C-----------------------------------------------------------------------
  80. C TEST FOR UNDERFLOW AND OVERFLOW
  81. C-----------------------------------------------------------------------
  82. IF (DABS(RS1).GT.ELIM) GO TO 60
  83. IF (KDFLG.EQ.1) KFLAG = 2
  84. IF (DABS(RS1).LT.ALIM) GO TO 40
  85. C-----------------------------------------------------------------------
  86. C REFINE TEST AND SCALE
  87. C-----------------------------------------------------------------------
  88. APHI = ZABS(COMPLEX(PHIR(J),PHII(J)))
  89. RS1 = RS1 + DLOG(APHI)
  90. IF (DABS(RS1).GT.ELIM) GO TO 60
  91. IF (KDFLG.EQ.1) KFLAG = 1
  92. IF (RS1.LT.0.0D0) GO TO 40
  93. IF (KDFLG.EQ.1) KFLAG = 3
  94. 40 CONTINUE
  95. C-----------------------------------------------------------------------
  96. C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
  97. C EXPONENT EXTREMES
  98. C-----------------------------------------------------------------------
  99. S2R = PHIR(J)*SUMR(J) - PHII(J)*SUMI(J)
  100. S2I = PHIR(J)*SUMI(J) + PHII(J)*SUMR(J)
  101. STR = DEXP(S1R)*CSSR(KFLAG)
  102. S1R = STR*DCOS(S1I)
  103. S1I = STR*DSIN(S1I)
  104. STR = S2R*S1R - S2I*S1I
  105. S2I = S1R*S2I + S2R*S1I
  106. S2R = STR
  107. IF (KFLAG.NE.1) GO TO 50
  108. CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
  109. IF (NW.NE.0) GO TO 60
  110. 50 CONTINUE
  111. CYR(KDFLG) = S2R
  112. CYI(KDFLG) = S2I
  113. YR(I) = S2R*CSRR(KFLAG)
  114. YI(I) = S2I*CSRR(KFLAG)
  115. IF (KDFLG.EQ.2) GO TO 75
  116. KDFLG = 2
  117. GO TO 70
  118. 60 CONTINUE
  119. IF (RS1.GT.0.0D0) GO TO 300
  120. C-----------------------------------------------------------------------
  121. C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
  122. C-----------------------------------------------------------------------
  123. IF (ZR.LT.0.0D0) GO TO 300
  124. KDFLG = 1
  125. YR(I)=ZEROR
  126. YI(I)=ZEROI
  127. NZ=NZ+1
  128. IF (I.EQ.1) GO TO 70
  129. IF ((YR(I-1).EQ.ZEROR).AND.(YI(I-1).EQ.ZEROI)) GO TO 70
  130. YR(I-1)=ZEROR
  131. YI(I-1)=ZEROI
  132. NZ=NZ+1
  133. 70 CONTINUE
  134. I = N
  135. 75 CONTINUE
  136. RAZR = 1.0D0/ZABS(COMPLEX(ZRR,ZRI))
  137. STR = ZRR*RAZR
  138. STI = -ZRI*RAZR
  139. RZR = (STR+STR)*RAZR
  140. RZI = (STI+STI)*RAZR
  141. CKR = FN*RZR
  142. CKI = FN*RZI
  143. IB = I + 1
  144. IF (N.LT.IB) GO TO 160
  145. C-----------------------------------------------------------------------
  146. C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
  147. C ON UNDERFLOW.
  148. C-----------------------------------------------------------------------
  149. FN = FNU + DBLE(FLOAT(N-1))
  150. IPARD = 1
  151. IF (MR.NE.0) IPARD = 0
  152. INITD = 0
  153. CALL ZUNIK(ZRR, ZRI, FN, 2, IPARD, TOL, INITD, PHIDR, PHIDI,
  154. * ZET1DR, ZET1DI, ZET2DR, ZET2DI, SUMDR, SUMDI, CWRKR(1,3),
  155. * CWRKI(1,3))
  156. IF (KODE.EQ.1) GO TO 80
  157. STR = ZRR + ZET2DR
  158. STI = ZRI + ZET2DI
  159. RAST = FN/ZABS(COMPLEX(STR,STI))
  160. STR = STR*RAST*RAST
  161. STI = -STI*RAST*RAST
  162. S1R = ZET1DR - STR
  163. S1I = ZET1DI - STI
  164. GO TO 90
  165. 80 CONTINUE
  166. S1R = ZET1DR - ZET2DR
  167. S1I = ZET1DI - ZET2DI
  168. 90 CONTINUE
  169. RS1 = S1R
  170. IF (DABS(RS1).GT.ELIM) GO TO 95
  171. IF (DABS(RS1).LT.ALIM) GO TO 100
  172. C----------------------------------------------------------------------------
  173. C REFINE ESTIMATE AND TEST
  174. C-------------------------------------------------------------------------
  175. APHI = ZABS(COMPLEX(PHIDR,PHIDI))
  176. RS1 = RS1+DLOG(APHI)
  177. IF (DABS(RS1).LT.ELIM) GO TO 100
  178. 95 CONTINUE
  179. IF (DABS(RS1).GT.0.0D0) GO TO 300
  180. C-----------------------------------------------------------------------
  181. C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
  182. C-----------------------------------------------------------------------
  183. IF (ZR.LT.0.0D0) GO TO 300
  184. NZ = N
  185. DO 96 I=1,N
  186. YR(I) = ZEROR
  187. YI(I) = ZEROI
  188. 96 CONTINUE
  189. RETURN
  190. C---------------------------------------------------------------------------
  191. C FORWARD RECUR FOR REMAINDER OF THE SEQUENCE
  192. C----------------------------------------------------------------------------
  193. 100 CONTINUE
  194. S1R = CYR(1)
  195. S1I = CYI(1)
  196. S2R = CYR(2)
  197. S2I = CYI(2)
  198. C1R = CSRR(KFLAG)
  199. ASCLE = BRY(KFLAG)
  200. DO 120 I=IB,N
  201. C2R = S2R
  202. C2I = S2I
  203. S2R = CKR*C2R - CKI*C2I + S1R
  204. S2I = CKR*C2I + CKI*C2R + S1I
  205. S1R = C2R
  206. S1I = C2I
  207. CKR = CKR + RZR
  208. CKI = CKI + RZI
  209. C2R = S2R*C1R
  210. C2I = S2I*C1R
  211. YR(I) = C2R
  212. YI(I) = C2I
  213. IF (KFLAG.GE.3) GO TO 120
  214. STR = DABS(C2R)
  215. STI = DABS(C2I)
  216. C2M = DMAX1(STR,STI)
  217. IF (C2M.LE.ASCLE) GO TO 120
  218. KFLAG = KFLAG + 1
  219. ASCLE = BRY(KFLAG)
  220. S1R = S1R*C1R
  221. S1I = S1I*C1R
  222. S2R = C2R
  223. S2I = C2I
  224. S1R = S1R*CSSR(KFLAG)
  225. S1I = S1I*CSSR(KFLAG)
  226. S2R = S2R*CSSR(KFLAG)
  227. S2I = S2I*CSSR(KFLAG)
  228. C1R = CSRR(KFLAG)
  229. 120 CONTINUE
  230. 160 CONTINUE
  231. IF (MR.EQ.0) RETURN
  232. C-----------------------------------------------------------------------
  233. C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
  234. C-----------------------------------------------------------------------
  235. NZ = 0
  236. FMR = DBLE(FLOAT(MR))
  237. SGN = -DSIGN(PI,FMR)
  238. C-----------------------------------------------------------------------
  239. C CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
  240. C-----------------------------------------------------------------------
  241. CSGNI = SGN
  242. INU = INT(SNGL(FNU))
  243. FNF = FNU - DBLE(FLOAT(INU))
  244. IFN = INU + N - 1
  245. ANG = FNF*SGN
  246. CSPNR = DCOS(ANG)
  247. CSPNI = DSIN(ANG)
  248. IF (MOD(IFN,2).EQ.0) GO TO 170
  249. CSPNR = -CSPNR
  250. CSPNI = -CSPNI
  251. 170 CONTINUE
  252. ASC = BRY(1)
  253. IUF = 0
  254. KK = N
  255. KDFLG = 1
  256. IB = IB - 1
  257. IC = IB - 1
  258. DO 270 K=1,N
  259. FN = FNU + DBLE(FLOAT(KK-1))
  260. C-----------------------------------------------------------------------
  261. C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
  262. C FUNCTION ABOVE
  263. C-----------------------------------------------------------------------
  264. M=3
  265. IF (N.GT.2) GO TO 175
  266. 172 CONTINUE
  267. INITD = INIT(J)
  268. PHIDR = PHIR(J)
  269. PHIDI = PHII(J)
  270. ZET1DR = ZETA1R(J)
  271. ZET1DI = ZETA1I(J)
  272. ZET2DR = ZETA2R(J)
  273. ZET2DI = ZETA2I(J)
  274. SUMDR = SUMR(J)
  275. SUMDI = SUMI(J)
  276. M = J
  277. J = 3 - J
  278. GO TO 180
  279. 175 CONTINUE
  280. IF ((KK.EQ.N).AND.(IB.LT.N)) GO TO 180
  281. IF ((KK.EQ.IB).OR.(KK.EQ.IC)) GO TO 172
  282. INITD = 0
  283. 180 CONTINUE
  284. CALL ZUNIK(ZRR, ZRI, FN, 1, 0, TOL, INITD, PHIDR, PHIDI,
  285. * ZET1DR, ZET1DI, ZET2DR, ZET2DI, SUMDR, SUMDI,
  286. * CWRKR(1,M), CWRKI(1,M))
  287. IF (KODE.EQ.1) GO TO 200
  288. STR = ZRR + ZET2DR
  289. STI = ZRI + ZET2DI
  290. RAST = FN/ZABS(COMPLEX(STR,STI))
  291. STR = STR*RAST*RAST
  292. STI = -STI*RAST*RAST
  293. S1R = -ZET1DR + STR
  294. S1I = -ZET1DI + STI
  295. GO TO 210
  296. 200 CONTINUE
  297. S1R = -ZET1DR + ZET2DR
  298. S1I = -ZET1DI + ZET2DI
  299. 210 CONTINUE
  300. C-----------------------------------------------------------------------
  301. C TEST FOR UNDERFLOW AND OVERFLOW
  302. C-----------------------------------------------------------------------
  303. RS1 = S1R
  304. IF (DABS(RS1).GT.ELIM) GO TO 260
  305. IF (KDFLG.EQ.1) IFLAG = 2
  306. IF (DABS(RS1).LT.ALIM) GO TO 220
  307. C-----------------------------------------------------------------------
  308. C REFINE TEST AND SCALE
  309. C-----------------------------------------------------------------------
  310. APHI = ZABS(COMPLEX(PHIDR,PHIDI))
  311. RS1 = RS1 + DLOG(APHI)
  312. IF (DABS(RS1).GT.ELIM) GO TO 260
  313. IF (KDFLG.EQ.1) IFLAG = 1
  314. IF (RS1.LT.0.0D0) GO TO 220
  315. IF (KDFLG.EQ.1) IFLAG = 3
  316. 220 CONTINUE
  317. STR = PHIDR*SUMDR - PHIDI*SUMDI
  318. STI = PHIDR*SUMDI + PHIDI*SUMDR
  319. S2R = -CSGNI*STI
  320. S2I = CSGNI*STR
  321. STR = DEXP(S1R)*CSSR(IFLAG)
  322. S1R = STR*DCOS(S1I)
  323. S1I = STR*DSIN(S1I)
  324. STR = S2R*S1R - S2I*S1I
  325. S2I = S2R*S1I + S2I*S1R
  326. S2R = STR
  327. IF (IFLAG.NE.1) GO TO 230
  328. CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
  329. IF (NW.EQ.0) GO TO 230
  330. S2R = ZEROR
  331. S2I = ZEROI
  332. 230 CONTINUE
  333. CYR(KDFLG) = S2R
  334. CYI(KDFLG) = S2I
  335. C2R = S2R
  336. C2I = S2I
  337. S2R = S2R*CSRR(IFLAG)
  338. S2I = S2I*CSRR(IFLAG)
  339. C-----------------------------------------------------------------------
  340. C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
  341. C-----------------------------------------------------------------------
  342. S1R = YR(KK)
  343. S1I = YI(KK)
  344. IF (KODE.EQ.1) GO TO 250
  345. CALL ZS1S2(ZRR, ZRI, S1R, S1I, S2R, S2I, NW, ASC, ALIM, IUF)
  346. NZ = NZ + NW
  347. 250 CONTINUE
  348. YR(KK) = S1R*CSPNR - S1I*CSPNI + S2R
  349. YI(KK) = CSPNR*S1I + CSPNI*S1R + S2I
  350. KK = KK - 1
  351. CSPNR = -CSPNR
  352. CSPNI = -CSPNI
  353. IF (C2R.NE.0.0D0 .OR. C2I.NE.0.0D0) GO TO 255
  354. KDFLG = 1
  355. GO TO 270
  356. 255 CONTINUE
  357. IF (KDFLG.EQ.2) GO TO 275
  358. KDFLG = 2
  359. GO TO 270
  360. 260 CONTINUE
  361. IF (RS1.GT.0.0D0) GO TO 300
  362. S2R = ZEROR
  363. S2I = ZEROI
  364. GO TO 230
  365. 270 CONTINUE
  366. K = N
  367. 275 CONTINUE
  368. IL = N - K
  369. IF (IL.EQ.0) RETURN
  370. C-----------------------------------------------------------------------
  371. C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
  372. C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
  373. C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
  374. C-----------------------------------------------------------------------
  375. S1R = CYR(1)
  376. S1I = CYI(1)
  377. S2R = CYR(2)
  378. S2I = CYI(2)
  379. CSR = CSRR(IFLAG)
  380. ASCLE = BRY(IFLAG)
  381. FN = DBLE(FLOAT(INU+IL))
  382. DO 290 I=1,IL
  383. C2R = S2R
  384. C2I = S2I
  385. S2R = S1R + (FN+FNF)*(RZR*C2R-RZI*C2I)
  386. S2I = S1I + (FN+FNF)*(RZR*C2I+RZI*C2R)
  387. S1R = C2R
  388. S1I = C2I
  389. FN = FN - 1.0D0
  390. C2R = S2R*CSR
  391. C2I = S2I*CSR
  392. CKR = C2R
  393. CKI = C2I
  394. C1R = YR(KK)
  395. C1I = YI(KK)
  396. IF (KODE.EQ.1) GO TO 280
  397. CALL ZS1S2(ZRR, ZRI, C1R, C1I, C2R, C2I, NW, ASC, ALIM, IUF)
  398. NZ = NZ + NW
  399. 280 CONTINUE
  400. YR(KK) = C1R*CSPNR - C1I*CSPNI + C2R
  401. YI(KK) = C1R*CSPNI + C1I*CSPNR + C2I
  402. KK = KK - 1
  403. CSPNR = -CSPNR
  404. CSPNI = -CSPNI
  405. IF (IFLAG.GE.3) GO TO 290
  406. C2R = DABS(CKR)
  407. C2I = DABS(CKI)
  408. C2M = DMAX1(C2R,C2I)
  409. IF (C2M.LE.ASCLE) GO TO 290
  410. IFLAG = IFLAG + 1
  411. ASCLE = BRY(IFLAG)
  412. S1R = S1R*CSR
  413. S1I = S1I*CSR
  414. S2R = CKR
  415. S2I = CKI
  416. S1R = S1R*CSSR(IFLAG)
  417. S1I = S1I*CSSR(IFLAG)
  418. S2R = S2R*CSSR(IFLAG)
  419. S2I = S2I*CSSR(IFLAG)
  420. CSR = CSRR(IFLAG)
  421. 290 CONTINUE
  422. RETURN
  423. 300 CONTINUE
  424. NZ = -1
  425. RETURN
  426. END