zunk2.f 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505
  1. SUBROUTINE ZUNK2(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM,
  2. * ALIM)
  3. C***BEGIN PROLOGUE ZUNK2
  4. C***REFER TO ZBESK
  5. C
  6. C ZUNK2 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 EXPANSIONS FOR H(KIND,FNU,ZN) AND J(FNU,ZN)
  9. C WHERE ZN IS IN THE RIGHT HALF PLANE, KIND=(3-MR)/2, MR=+1 OR
  10. C -1. HERE ZN=ZR*I OR -ZR*I WHERE ZR=Z IF Z IS IN THE RIGHT
  11. C HALF PLANE OR ZR=-Z IF Z IS IN THE LEFT HALF PLANE. MR INDIC-
  12. C ATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
  13. C NZ=-1 MEANS AN OVERFLOW WILL OCCUR
  14. C
  15. C***ROUTINES CALLED ZAIRY,ZKSCL,ZS1S2,ZUCHK,ZUNHJ,D1MACH,ZABS
  16. C***END PROLOGUE ZUNK2
  17. C COMPLEX AI,ARG,ARGD,ASUM,ASUMD,BSUM,BSUMD,CFN,CI,CIP,CK,CONE,CRSC,
  18. C *CR1,CR2,CS,CSCL,CSGN,CSPN,CSR,CSS,CY,CZERO,C1,C2,DAI,PHI,PHID,RZ,
  19. C *S1,S2,Y,Z,ZB,ZETA1,ZETA1D,ZETA2,ZETA2D,ZN,ZR
  20. DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGDI,
  21. * ARGDR, ARGI, ARGR, ASC, ASCLE, ASUMDI, ASUMDR, ASUMI, ASUMR,
  22. * BRY, BSUMDI, BSUMDR, BSUMI, BSUMR, CAR, CIPI, CIPR, CKI, CKR,
  23. * CONER, CRSC, CR1I, CR1R, CR2I, CR2R, CSCL, CSGNI, CSI,
  24. * CSPNI, CSPNR, CSR, CSRR, CSSR, CYI, CYR, C1I, C1R, C2I, C2M,
  25. * C2R, DAII, DAIR, ELIM, FMR, FN, FNF, FNU, HPI, PHIDI, PHIDR,
  26. * PHII, PHIR, PI, PTI, PTR, RAST, RAZR, RS1, RZI, RZR, SAR, SGN,
  27. * STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, YY, ZBI, ZBR, ZEROI,
  28. * ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZET1DI, ZET1DR, ZET2DI,
  29. * ZET2DR, ZI, ZNI, ZNR, ZR, ZRI, ZRR, D1MACH, ZABS
  30. INTEGER I, IB, IFLAG, IFN, IL, IN, INU, IUF, K, KDFLG, KFLAG, KK,
  31. * KODE, MR, N, NAI, NDAI, NW, NZ, IDUM, J, IPARD, IC
  32. DIMENSION BRY(3), YR(N), YI(N), ASUMR(2), ASUMI(2), BSUMR(2),
  33. * BSUMI(2), PHIR(2), PHII(2), ARGR(2), ARGI(2), ZETA1R(2),
  34. * ZETA1I(2), ZETA2R(2), ZETA2I(2), CYR(2), CYI(2), CIPR(4),
  35. * CIPI(4), CSSR(3), CSRR(3)
  36. DATA ZEROR,ZEROI,CONER,CR1R,CR1I,CR2R,CR2I /
  37. 1 0.0D0, 0.0D0, 1.0D0,
  38. 1 1.0D0,1.73205080756887729D0 , -0.5D0,-8.66025403784438647D-01 /
  39. DATA HPI, PI, AIC /
  40. 1 1.57079632679489662D+00, 3.14159265358979324D+00,
  41. 1 1.26551212348464539D+00/
  42. DATA CIPR(1),CIPI(1),CIPR(2),CIPI(2),CIPR(3),CIPI(3),CIPR(4),
  43. * CIPI(4) /
  44. 1 1.0D0,0.0D0 , 0.0D0,-1.0D0 , -1.0D0,0.0D0 , 0.0D0,1.0D0 /
  45. C
  46. KDFLG = 1
  47. NZ = 0
  48. C-----------------------------------------------------------------------
  49. C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
  50. C THE UNDERFLOW LIMIT
  51. C-----------------------------------------------------------------------
  52. CSCL = 1.0D0/TOL
  53. CRSC = TOL
  54. CSSR(1) = CSCL
  55. CSSR(2) = CONER
  56. CSSR(3) = CRSC
  57. CSRR(1) = CRSC
  58. CSRR(2) = CONER
  59. CSRR(3) = CSCL
  60. BRY(1) = 1.0D+3*D1MACH(1)/TOL
  61. BRY(2) = 1.0D0/BRY(1)
  62. BRY(3) = D1MACH(2)
  63. ZRR = ZR
  64. ZRI = ZI
  65. IF (ZR.GE.0.0D0) GO TO 10
  66. ZRR = -ZR
  67. ZRI = -ZI
  68. 10 CONTINUE
  69. YY = ZRI
  70. ZNR = ZRI
  71. ZNI = -ZRR
  72. ZBR = ZRR
  73. ZBI = ZRI
  74. INU = INT(SNGL(FNU))
  75. FNF = FNU - DBLE(FLOAT(INU))
  76. ANG = -HPI*FNF
  77. CAR = DCOS(ANG)
  78. SAR = DSIN(ANG)
  79. C2R = HPI*SAR
  80. C2I = -HPI*CAR
  81. KK = MOD(INU,4) + 1
  82. STR = C2R*CIPR(KK) - C2I*CIPI(KK)
  83. STI = C2R*CIPI(KK) + C2I*CIPR(KK)
  84. CSR = CR1R*STR - CR1I*STI
  85. CSI = CR1R*STI + CR1I*STR
  86. IF (YY.GT.0.0D0) GO TO 20
  87. ZNR = -ZNR
  88. ZBI = -ZBI
  89. 20 CONTINUE
  90. C-----------------------------------------------------------------------
  91. C K(FNU,Z) IS COMPUTED FROM H(2,FNU,-I*Z) WHERE Z IS IN THE FIRST
  92. C QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
  93. C CONJUGATION SINCE THE K FUNCTION IS REAL ON THE POSITIVE REAL AXIS
  94. C-----------------------------------------------------------------------
  95. J = 2
  96. DO 80 I=1,N
  97. C-----------------------------------------------------------------------
  98. C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
  99. C-----------------------------------------------------------------------
  100. J = 3 - J
  101. FN = FNU + DBLE(FLOAT(I-1))
  102. CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIR(J), PHII(J), ARGR(J),
  103. * ARGI(J), ZETA1R(J), ZETA1I(J), ZETA2R(J), ZETA2I(J), ASUMR(J),
  104. * ASUMI(J), BSUMR(J), BSUMI(J))
  105. IF (KODE.EQ.1) GO TO 30
  106. STR = ZBR + ZETA2R(J)
  107. STI = ZBI + ZETA2I(J)
  108. RAST = FN/ZABS(COMPLEX(STR,STI))
  109. STR = STR*RAST*RAST
  110. STI = -STI*RAST*RAST
  111. S1R = ZETA1R(J) - STR
  112. S1I = ZETA1I(J) - STI
  113. GO TO 40
  114. 30 CONTINUE
  115. S1R = ZETA1R(J) - ZETA2R(J)
  116. S1I = ZETA1I(J) - ZETA2I(J)
  117. 40 CONTINUE
  118. C-----------------------------------------------------------------------
  119. C TEST FOR UNDERFLOW AND OVERFLOW
  120. C-----------------------------------------------------------------------
  121. RS1 = S1R
  122. IF (DABS(RS1).GT.ELIM) GO TO 70
  123. IF (KDFLG.EQ.1) KFLAG = 2
  124. IF (DABS(RS1).LT.ALIM) GO TO 50
  125. C-----------------------------------------------------------------------
  126. C REFINE TEST AND SCALE
  127. C-----------------------------------------------------------------------
  128. APHI = ZABS(COMPLEX(PHIR(J),PHII(J)))
  129. AARG = ZABS(COMPLEX(ARGR(J),ARGI(J)))
  130. RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC
  131. IF (DABS(RS1).GT.ELIM) GO TO 70
  132. IF (KDFLG.EQ.1) KFLAG = 1
  133. IF (RS1.LT.0.0D0) GO TO 50
  134. IF (KDFLG.EQ.1) KFLAG = 3
  135. 50 CONTINUE
  136. C-----------------------------------------------------------------------
  137. C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
  138. C EXPONENT EXTREMES
  139. C-----------------------------------------------------------------------
  140. C2R = ARGR(J)*CR2R - ARGI(J)*CR2I
  141. C2I = ARGR(J)*CR2I + ARGI(J)*CR2R
  142. CALL ZAIRY(C2R, C2I, 0, 2, AIR, AII, NAI, IDUM)
  143. CALL ZAIRY(C2R, C2I, 1, 2, DAIR, DAII, NDAI, IDUM)
  144. STR = DAIR*BSUMR(J) - DAII*BSUMI(J)
  145. STI = DAIR*BSUMI(J) + DAII*BSUMR(J)
  146. PTR = STR*CR2R - STI*CR2I
  147. PTI = STR*CR2I + STI*CR2R
  148. STR = PTR + (AIR*ASUMR(J)-AII*ASUMI(J))
  149. STI = PTI + (AIR*ASUMI(J)+AII*ASUMR(J))
  150. PTR = STR*PHIR(J) - STI*PHII(J)
  151. PTI = STR*PHII(J) + STI*PHIR(J)
  152. S2R = PTR*CSR - PTI*CSI
  153. S2I = PTR*CSI + PTI*CSR
  154. STR = DEXP(S1R)*CSSR(KFLAG)
  155. S1R = STR*DCOS(S1I)
  156. S1I = STR*DSIN(S1I)
  157. STR = S2R*S1R - S2I*S1I
  158. S2I = S1R*S2I + S2R*S1I
  159. S2R = STR
  160. IF (KFLAG.NE.1) GO TO 60
  161. CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
  162. IF (NW.NE.0) GO TO 70
  163. 60 CONTINUE
  164. IF (YY.LE.0.0D0) S2I = -S2I
  165. CYR(KDFLG) = S2R
  166. CYI(KDFLG) = S2I
  167. YR(I) = S2R*CSRR(KFLAG)
  168. YI(I) = S2I*CSRR(KFLAG)
  169. STR = CSI
  170. CSI = -CSR
  171. CSR = STR
  172. IF (KDFLG.EQ.2) GO TO 85
  173. KDFLG = 2
  174. GO TO 80
  175. 70 CONTINUE
  176. IF (RS1.GT.0.0D0) GO TO 320
  177. C-----------------------------------------------------------------------
  178. C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
  179. C-----------------------------------------------------------------------
  180. IF (ZR.LT.0.0D0) GO TO 320
  181. KDFLG = 1
  182. YR(I)=ZEROR
  183. YI(I)=ZEROI
  184. NZ=NZ+1
  185. STR = CSI
  186. CSI =-CSR
  187. CSR = STR
  188. IF (I.EQ.1) GO TO 80
  189. IF ((YR(I-1).EQ.ZEROR).AND.(YI(I-1).EQ.ZEROI)) GO TO 80
  190. YR(I-1)=ZEROR
  191. YI(I-1)=ZEROI
  192. NZ=NZ+1
  193. 80 CONTINUE
  194. I = N
  195. 85 CONTINUE
  196. RAZR = 1.0D0/ZABS(COMPLEX(ZRR,ZRI))
  197. STR = ZRR*RAZR
  198. STI = -ZRI*RAZR
  199. RZR = (STR+STR)*RAZR
  200. RZI = (STI+STI)*RAZR
  201. CKR = FN*RZR
  202. CKI = FN*RZI
  203. IB = I + 1
  204. IF (N.LT.IB) GO TO 180
  205. C-----------------------------------------------------------------------
  206. C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
  207. C ON UNDERFLOW.
  208. C-----------------------------------------------------------------------
  209. FN = FNU + DBLE(FLOAT(N-1))
  210. IPARD = 1
  211. IF (MR.NE.0) IPARD = 0
  212. CALL ZUNHJ(ZNR, ZNI, FN, IPARD, TOL, PHIDR, PHIDI, ARGDR, ARGDI,
  213. * ZET1DR, ZET1DI, ZET2DR, ZET2DI, ASUMDR, ASUMDI, BSUMDR, BSUMDI)
  214. IF (KODE.EQ.1) GO TO 90
  215. STR = ZBR + ZET2DR
  216. STI = ZBI + ZET2DI
  217. RAST = FN/ZABS(COMPLEX(STR,STI))
  218. STR = STR*RAST*RAST
  219. STI = -STI*RAST*RAST
  220. S1R = ZET1DR - STR
  221. S1I = ZET1DI - STI
  222. GO TO 100
  223. 90 CONTINUE
  224. S1R = ZET1DR - ZET2DR
  225. S1I = ZET1DI - ZET2DI
  226. 100 CONTINUE
  227. RS1 = S1R
  228. IF (DABS(RS1).GT.ELIM) GO TO 105
  229. IF (DABS(RS1).LT.ALIM) GO TO 120
  230. C----------------------------------------------------------------------------
  231. C REFINE ESTIMATE AND TEST
  232. C-------------------------------------------------------------------------
  233. APHI = ZABS(COMPLEX(PHIDR,PHIDI))
  234. RS1 = RS1+DLOG(APHI)
  235. IF (DABS(RS1).LT.ELIM) GO TO 120
  236. 105 CONTINUE
  237. IF (RS1.GT.0.0D0) GO TO 320
  238. C-----------------------------------------------------------------------
  239. C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
  240. C-----------------------------------------------------------------------
  241. IF (ZR.LT.0.0D0) GO TO 320
  242. NZ = N
  243. DO 106 I=1,N
  244. YR(I) = ZEROR
  245. YI(I) = ZEROI
  246. 106 CONTINUE
  247. RETURN
  248. 120 CONTINUE
  249. S1R = CYR(1)
  250. S1I = CYI(1)
  251. S2R = CYR(2)
  252. S2I = CYI(2)
  253. C1R = CSRR(KFLAG)
  254. ASCLE = BRY(KFLAG)
  255. DO 130 I=IB,N
  256. C2R = S2R
  257. C2I = S2I
  258. S2R = CKR*C2R - CKI*C2I + S1R
  259. S2I = CKR*C2I + CKI*C2R + S1I
  260. S1R = C2R
  261. S1I = C2I
  262. CKR = CKR + RZR
  263. CKI = CKI + RZI
  264. C2R = S2R*C1R
  265. C2I = S2I*C1R
  266. YR(I) = C2R
  267. YI(I) = C2I
  268. IF (KFLAG.GE.3) GO TO 130
  269. STR = DABS(C2R)
  270. STI = DABS(C2I)
  271. C2M = DMAX1(STR,STI)
  272. IF (C2M.LE.ASCLE) GO TO 130
  273. KFLAG = KFLAG + 1
  274. ASCLE = BRY(KFLAG)
  275. S1R = S1R*C1R
  276. S1I = S1I*C1R
  277. S2R = C2R
  278. S2I = C2I
  279. S1R = S1R*CSSR(KFLAG)
  280. S1I = S1I*CSSR(KFLAG)
  281. S2R = S2R*CSSR(KFLAG)
  282. S2I = S2I*CSSR(KFLAG)
  283. C1R = CSRR(KFLAG)
  284. 130 CONTINUE
  285. 180 CONTINUE
  286. IF (MR.EQ.0) RETURN
  287. C-----------------------------------------------------------------------
  288. C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
  289. C-----------------------------------------------------------------------
  290. NZ = 0
  291. FMR = DBLE(FLOAT(MR))
  292. SGN = -DSIGN(PI,FMR)
  293. C-----------------------------------------------------------------------
  294. C CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP.
  295. C-----------------------------------------------------------------------
  296. CSGNI = SGN
  297. IF (YY.LE.0.0D0) CSGNI = -CSGNI
  298. IFN = INU + N - 1
  299. ANG = FNF*SGN
  300. CSPNR = DCOS(ANG)
  301. CSPNI = DSIN(ANG)
  302. IF (MOD(IFN,2).EQ.0) GO TO 190
  303. CSPNR = -CSPNR
  304. CSPNI = -CSPNI
  305. 190 CONTINUE
  306. C-----------------------------------------------------------------------
  307. C CS=COEFF OF THE J FUNCTION TO GET THE I FUNCTION. I(FNU,Z) IS
  308. C COMPUTED FROM EXP(I*FNU*HPI)*J(FNU,-I*Z) WHERE Z IS IN THE FIRST
  309. C QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
  310. C CONJUGATION SINCE THE I FUNCTION IS REAL ON THE POSITIVE REAL AXIS
  311. C-----------------------------------------------------------------------
  312. CSR = SAR*CSGNI
  313. CSI = CAR*CSGNI
  314. IN = MOD(IFN,4) + 1
  315. C2R = CIPR(IN)
  316. C2I = CIPI(IN)
  317. STR = CSR*C2R + CSI*C2I
  318. CSI = -CSR*C2I + CSI*C2R
  319. CSR = STR
  320. ASC = BRY(1)
  321. IUF = 0
  322. KK = N
  323. KDFLG = 1
  324. IB = IB - 1
  325. IC = IB - 1
  326. DO 290 K=1,N
  327. FN = FNU + DBLE(FLOAT(KK-1))
  328. C-----------------------------------------------------------------------
  329. C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
  330. C FUNCTION ABOVE
  331. C-----------------------------------------------------------------------
  332. IF (N.GT.2) GO TO 175
  333. 172 CONTINUE
  334. PHIDR = PHIR(J)
  335. PHIDI = PHII(J)
  336. ARGDR = ARGR(J)
  337. ARGDI = ARGI(J)
  338. ZET1DR = ZETA1R(J)
  339. ZET1DI = ZETA1I(J)
  340. ZET2DR = ZETA2R(J)
  341. ZET2DI = ZETA2I(J)
  342. ASUMDR = ASUMR(J)
  343. ASUMDI = ASUMI(J)
  344. BSUMDR = BSUMR(J)
  345. BSUMDI = BSUMI(J)
  346. J = 3 - J
  347. GO TO 210
  348. 175 CONTINUE
  349. IF ((KK.EQ.N).AND.(IB.LT.N)) GO TO 210
  350. IF ((KK.EQ.IB).OR.(KK.EQ.IC)) GO TO 172
  351. CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIDR, PHIDI, ARGDR,
  352. * ARGDI, ZET1DR, ZET1DI, ZET2DR, ZET2DI, ASUMDR,
  353. * ASUMDI, BSUMDR, BSUMDI)
  354. 210 CONTINUE
  355. IF (KODE.EQ.1) GO TO 220
  356. STR = ZBR + ZET2DR
  357. STI = ZBI + ZET2DI
  358. RAST = FN/ZABS(COMPLEX(STR,STI))
  359. STR = STR*RAST*RAST
  360. STI = -STI*RAST*RAST
  361. S1R = -ZET1DR + STR
  362. S1I = -ZET1DI + STI
  363. GO TO 230
  364. 220 CONTINUE
  365. S1R = -ZET1DR + ZET2DR
  366. S1I = -ZET1DI + ZET2DI
  367. 230 CONTINUE
  368. C-----------------------------------------------------------------------
  369. C TEST FOR UNDERFLOW AND OVERFLOW
  370. C-----------------------------------------------------------------------
  371. RS1 = S1R
  372. IF (DABS(RS1).GT.ELIM) GO TO 280
  373. IF (KDFLG.EQ.1) IFLAG = 2
  374. IF (DABS(RS1).LT.ALIM) GO TO 240
  375. C-----------------------------------------------------------------------
  376. C REFINE TEST AND SCALE
  377. C-----------------------------------------------------------------------
  378. APHI = ZABS(COMPLEX(PHIDR,PHIDI))
  379. AARG = ZABS(COMPLEX(ARGDR,ARGDI))
  380. RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC
  381. IF (DABS(RS1).GT.ELIM) GO TO 280
  382. IF (KDFLG.EQ.1) IFLAG = 1
  383. IF (RS1.LT.0.0D0) GO TO 240
  384. IF (KDFLG.EQ.1) IFLAG = 3
  385. 240 CONTINUE
  386. CALL ZAIRY(ARGDR, ARGDI, 0, 2, AIR, AII, NAI, IDUM)
  387. CALL ZAIRY(ARGDR, ARGDI, 1, 2, DAIR, DAII, NDAI, IDUM)
  388. STR = DAIR*BSUMDR - DAII*BSUMDI
  389. STI = DAIR*BSUMDI + DAII*BSUMDR
  390. STR = STR + (AIR*ASUMDR-AII*ASUMDI)
  391. STI = STI + (AIR*ASUMDI+AII*ASUMDR)
  392. PTR = STR*PHIDR - STI*PHIDI
  393. PTI = STR*PHIDI + STI*PHIDR
  394. S2R = PTR*CSR - PTI*CSI
  395. S2I = PTR*CSI + PTI*CSR
  396. STR = DEXP(S1R)*CSSR(IFLAG)
  397. S1R = STR*DCOS(S1I)
  398. S1I = STR*DSIN(S1I)
  399. STR = S2R*S1R - S2I*S1I
  400. S2I = S2R*S1I + S2I*S1R
  401. S2R = STR
  402. IF (IFLAG.NE.1) GO TO 250
  403. CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
  404. IF (NW.EQ.0) GO TO 250
  405. S2R = ZEROR
  406. S2I = ZEROI
  407. 250 CONTINUE
  408. IF (YY.LE.0.0D0) S2I = -S2I
  409. CYR(KDFLG) = S2R
  410. CYI(KDFLG) = S2I
  411. C2R = S2R
  412. C2I = S2I
  413. S2R = S2R*CSRR(IFLAG)
  414. S2I = S2I*CSRR(IFLAG)
  415. C-----------------------------------------------------------------------
  416. C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
  417. C-----------------------------------------------------------------------
  418. S1R = YR(KK)
  419. S1I = YI(KK)
  420. IF (KODE.EQ.1) GO TO 270
  421. CALL ZS1S2(ZRR, ZRI, S1R, S1I, S2R, S2I, NW, ASC, ALIM, IUF)
  422. NZ = NZ + NW
  423. 270 CONTINUE
  424. YR(KK) = S1R*CSPNR - S1I*CSPNI + S2R
  425. YI(KK) = S1R*CSPNI + S1I*CSPNR + S2I
  426. KK = KK - 1
  427. CSPNR = -CSPNR
  428. CSPNI = -CSPNI
  429. STR = CSI
  430. CSI = -CSR
  431. CSR = STR
  432. IF (C2R.NE.0.0D0 .OR. C2I.NE.0.0D0) GO TO 255
  433. KDFLG = 1
  434. GO TO 290
  435. 255 CONTINUE
  436. IF (KDFLG.EQ.2) GO TO 295
  437. KDFLG = 2
  438. GO TO 290
  439. 280 CONTINUE
  440. IF (RS1.GT.0.0D0) GO TO 320
  441. S2R = ZEROR
  442. S2I = ZEROI
  443. GO TO 250
  444. 290 CONTINUE
  445. K = N
  446. 295 CONTINUE
  447. IL = N - K
  448. IF (IL.EQ.0) RETURN
  449. C-----------------------------------------------------------------------
  450. C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
  451. C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
  452. C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
  453. C-----------------------------------------------------------------------
  454. S1R = CYR(1)
  455. S1I = CYI(1)
  456. S2R = CYR(2)
  457. S2I = CYI(2)
  458. CSR = CSRR(IFLAG)
  459. ASCLE = BRY(IFLAG)
  460. FN = DBLE(FLOAT(INU+IL))
  461. DO 310 I=1,IL
  462. C2R = S2R
  463. C2I = S2I
  464. S2R = S1R + (FN+FNF)*(RZR*C2R-RZI*C2I)
  465. S2I = S1I + (FN+FNF)*(RZR*C2I+RZI*C2R)
  466. S1R = C2R
  467. S1I = C2I
  468. FN = FN - 1.0D0
  469. C2R = S2R*CSR
  470. C2I = S2I*CSR
  471. CKR = C2R
  472. CKI = C2I
  473. C1R = YR(KK)
  474. C1I = YI(KK)
  475. IF (KODE.EQ.1) GO TO 300
  476. CALL ZS1S2(ZRR, ZRI, C1R, C1I, C2R, C2I, NW, ASC, ALIM, IUF)
  477. NZ = NZ + NW
  478. 300 CONTINUE
  479. YR(KK) = C1R*CSPNR - C1I*CSPNI + C2R
  480. YI(KK) = C1R*CSPNI + C1I*CSPNR + C2I
  481. KK = KK - 1
  482. CSPNR = -CSPNR
  483. CSPNI = -CSPNI
  484. IF (IFLAG.GE.3) GO TO 310
  485. C2R = DABS(CKR)
  486. C2I = DABS(CKI)
  487. C2M = DMAX1(C2R,C2I)
  488. IF (C2M.LE.ASCLE) GO TO 310
  489. IFLAG = IFLAG + 1
  490. ASCLE = BRY(IFLAG)
  491. S1R = S1R*CSR
  492. S1I = S1I*CSR
  493. S2R = CKR
  494. S2I = CKI
  495. S1R = S1R*CSSR(IFLAG)
  496. S1I = S1I*CSSR(IFLAG)
  497. S2R = S2R*CSSR(IFLAG)
  498. S2I = S2I*CSSR(IFLAG)
  499. CSR = CSRR(IFLAG)
  500. 310 CONTINUE
  501. RETURN
  502. 320 CONTINUE
  503. NZ = -1
  504. RETURN
  505. END