zbknu.f 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581
  1. *DECK ZBKNU
  2. SUBROUTINE ZBKNU (ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM,
  3. + ALIM)
  4. C***BEGIN PROLOGUE ZBKNU
  5. C***SUBSIDIARY
  6. C***PURPOSE Subsidiary to ZAIRY, ZBESH, ZBESI and ZBESK
  7. C***LIBRARY SLATEC
  8. C***TYPE ALL (CBKNU-A, ZBKNU-A)
  9. C***AUTHOR Amos, D. E., (SNL)
  10. C***DESCRIPTION
  11. C
  12. C ZBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE.
  13. C
  14. C***SEE ALSO ZAIRY, ZBESH, ZBESI, ZBESK
  15. C***ROUTINES CALLED D1MACH, DGAMLN, I1MACH, ZABS, ZDIV, ZEXP, ZKSCL,
  16. C ZLOG, ZMLT, ZSHCH, ZSQRT, ZUCHK
  17. C***REVISION HISTORY (YYMMDD)
  18. C 830501 DATE WRITTEN
  19. C 910415 Prologue converted to Version 4.0 format. (BAB)
  20. C 930122 Added ZEXP, ZLOG and ZSQRT to EXTERNAL statement. (RWC)
  21. C***END PROLOGUE ZBKNU
  22. C
  23. DOUBLE PRECISION AA, AK, ALIM, ASCLE, A1, A2, BB, BK, BRY, CAZ,
  24. * CBI, CBR, CC, CCHI, CCHR, CKI, CKR, COEFI, COEFR, CONEI, CONER,
  25. * CRSCR, CSCLR, CSHI, CSHR, CSI, CSR, CSRR, CSSR, CTWOR,
  26. * CZEROI, CZEROR, CZI, CZR, DNU, DNU2, DPI, ELIM, ETEST, FC, FHS,
  27. * FI, FK, FKS, FMUI, FMUR, FNU, FPI, FR, G1, G2, HPI, PI, PR, PTI,
  28. * PTR, P1I, P1R, P2I, P2M, P2R, QI, QR, RAK, RCAZ, RTHPI, RZI,
  29. * RZR, R1, S, SMUI, SMUR, SPI, STI, STR, S1I, S1R, S2I, S2R, TM,
  30. * TOL, TTH, T1, T2, YI, YR, ZI, ZR, DGAMLN, D1MACH, ZABS, ELM,
  31. * CELMR, ZDR, ZDI, AS, ALAS, HELIM, CYR, CYI
  32. INTEGER I, IFLAG, INU, K, KFLAG, KK, KMAX, KODE, KODED, N, NZ,
  33. * IDUM, I1MACH, J, IC, INUB, NW
  34. DIMENSION YR(N), YI(N), CC(8), CSSR(3), CSRR(3), BRY(3), CYR(2),
  35. * CYI(2)
  36. EXTERNAL ZABS, ZEXP, ZLOG, ZSQRT
  37. C COMPLEX Z,Y,A,B,RZ,SMU,FU,FMU,F,FLRZ,CZ,S1,S2,CSH,CCH
  38. C COMPLEX CK,P,Q,COEF,P1,P2,CBK,PT,CZERO,CONE,CTWO,ST,EZ,CS,DK
  39. C
  40. DATA KMAX / 30 /
  41. DATA CZEROR,CZEROI,CONER,CONEI,CTWOR,R1/
  42. 1 0.0D0 , 0.0D0 , 1.0D0 , 0.0D0 , 2.0D0 , 2.0D0 /
  43. DATA DPI, RTHPI, SPI ,HPI, FPI, TTH /
  44. 1 3.14159265358979324D0, 1.25331413731550025D0,
  45. 2 1.90985931710274403D0, 1.57079632679489662D0,
  46. 3 1.89769999331517738D0, 6.66666666666666666D-01/
  47. DATA CC(1), CC(2), CC(3), CC(4), CC(5), CC(6), CC(7), CC(8)/
  48. 1 5.77215664901532861D-01, -4.20026350340952355D-02,
  49. 2 -4.21977345555443367D-02, 7.21894324666309954D-03,
  50. 3 -2.15241674114950973D-04, -2.01348547807882387D-05,
  51. 4 1.13302723198169588D-06, 6.11609510448141582D-09/
  52. C***FIRST EXECUTABLE STATEMENT ZBKNU
  53. CAZ = ZABS(ZR,ZI)
  54. CSCLR = 1.0D0/TOL
  55. CRSCR = TOL
  56. CSSR(1) = CSCLR
  57. CSSR(2) = 1.0D0
  58. CSSR(3) = CRSCR
  59. CSRR(1) = CRSCR
  60. CSRR(2) = 1.0D0
  61. CSRR(3) = CSCLR
  62. BRY(1) = 1.0D+3*D1MACH(1)/TOL
  63. BRY(2) = 1.0D0/BRY(1)
  64. BRY(3) = D1MACH(2)
  65. NZ = 0
  66. IFLAG = 0
  67. KODED = KODE
  68. RCAZ = 1.0D0/CAZ
  69. STR = ZR*RCAZ
  70. STI = -ZI*RCAZ
  71. RZR = (STR+STR)*RCAZ
  72. RZI = (STI+STI)*RCAZ
  73. INU = FNU+0.5D0
  74. DNU = FNU - INU
  75. IF (ABS(DNU).EQ.0.5D0) GO TO 110
  76. DNU2 = 0.0D0
  77. IF (ABS(DNU).GT.TOL) DNU2 = DNU*DNU
  78. IF (CAZ.GT.R1) GO TO 110
  79. C-----------------------------------------------------------------------
  80. C SERIES FOR ABS(Z).LE.R1
  81. C-----------------------------------------------------------------------
  82. FC = 1.0D0
  83. CALL ZLOG(RZR, RZI, SMUR, SMUI, IDUM)
  84. FMUR = SMUR*DNU
  85. FMUI = SMUI*DNU
  86. CALL ZSHCH(FMUR, FMUI, CSHR, CSHI, CCHR, CCHI)
  87. IF (DNU.EQ.0.0D0) GO TO 10
  88. FC = DNU*DPI
  89. FC = FC/SIN(FC)
  90. SMUR = CSHR/DNU
  91. SMUI = CSHI/DNU
  92. 10 CONTINUE
  93. A2 = 1.0D0 + DNU
  94. C-----------------------------------------------------------------------
  95. C GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU)
  96. C-----------------------------------------------------------------------
  97. T2 = EXP(-DGAMLN(A2,IDUM))
  98. T1 = 1.0D0/(T2*FC)
  99. IF (ABS(DNU).GT.0.1D0) GO TO 40
  100. C-----------------------------------------------------------------------
  101. C SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
  102. C-----------------------------------------------------------------------
  103. AK = 1.0D0
  104. S = CC(1)
  105. DO 20 K=2,8
  106. AK = AK*DNU2
  107. TM = CC(K)*AK
  108. S = S + TM
  109. IF (ABS(TM).LT.TOL) GO TO 30
  110. 20 CONTINUE
  111. 30 G1 = -S
  112. GO TO 50
  113. 40 CONTINUE
  114. G1 = (T1-T2)/(DNU+DNU)
  115. 50 CONTINUE
  116. G2 = (T1+T2)*0.5D0
  117. FR = FC*(CCHR*G1+SMUR*G2)
  118. FI = FC*(CCHI*G1+SMUI*G2)
  119. CALL ZEXP(FMUR, FMUI, STR, STI)
  120. PR = 0.5D0*STR/T2
  121. PI = 0.5D0*STI/T2
  122. CALL ZDIV(0.5D0, 0.0D0, STR, STI, PTR, PTI)
  123. QR = PTR/T1
  124. QI = PTI/T1
  125. S1R = FR
  126. S1I = FI
  127. S2R = PR
  128. S2I = PI
  129. AK = 1.0D0
  130. A1 = 1.0D0
  131. CKR = CONER
  132. CKI = CONEI
  133. BK = 1.0D0 - DNU2
  134. IF (INU.GT.0 .OR. N.GT.1) GO TO 80
  135. C-----------------------------------------------------------------------
  136. C GENERATE K(FNU,Z), 0.0D0 .LE. FNU .LT. 0.5D0 AND N=1
  137. C-----------------------------------------------------------------------
  138. IF (CAZ.LT.TOL) GO TO 70
  139. CALL ZMLT(ZR, ZI, ZR, ZI, CZR, CZI)
  140. CZR = 0.25D0*CZR
  141. CZI = 0.25D0*CZI
  142. T1 = 0.25D0*CAZ*CAZ
  143. 60 CONTINUE
  144. FR = (FR*AK+PR+QR)/BK
  145. FI = (FI*AK+PI+QI)/BK
  146. STR = 1.0D0/(AK-DNU)
  147. PR = PR*STR
  148. PI = PI*STR
  149. STR = 1.0D0/(AK+DNU)
  150. QR = QR*STR
  151. QI = QI*STR
  152. STR = CKR*CZR - CKI*CZI
  153. RAK = 1.0D0/AK
  154. CKI = (CKR*CZI+CKI*CZR)*RAK
  155. CKR = STR*RAK
  156. S1R = CKR*FR - CKI*FI + S1R
  157. S1I = CKR*FI + CKI*FR + S1I
  158. A1 = A1*T1*RAK
  159. BK = BK + AK + AK + 1.0D0
  160. AK = AK + 1.0D0
  161. IF (A1.GT.TOL) GO TO 60
  162. 70 CONTINUE
  163. YR(1) = S1R
  164. YI(1) = S1I
  165. IF (KODED.EQ.1) RETURN
  166. CALL ZEXP(ZR, ZI, STR, STI)
  167. CALL ZMLT(S1R, S1I, STR, STI, YR(1), YI(1))
  168. RETURN
  169. C-----------------------------------------------------------------------
  170. C GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE
  171. C-----------------------------------------------------------------------
  172. 80 CONTINUE
  173. IF (CAZ.LT.TOL) GO TO 100
  174. CALL ZMLT(ZR, ZI, ZR, ZI, CZR, CZI)
  175. CZR = 0.25D0*CZR
  176. CZI = 0.25D0*CZI
  177. T1 = 0.25D0*CAZ*CAZ
  178. 90 CONTINUE
  179. FR = (FR*AK+PR+QR)/BK
  180. FI = (FI*AK+PI+QI)/BK
  181. STR = 1.0D0/(AK-DNU)
  182. PR = PR*STR
  183. PI = PI*STR
  184. STR = 1.0D0/(AK+DNU)
  185. QR = QR*STR
  186. QI = QI*STR
  187. STR = CKR*CZR - CKI*CZI
  188. RAK = 1.0D0/AK
  189. CKI = (CKR*CZI+CKI*CZR)*RAK
  190. CKR = STR*RAK
  191. S1R = CKR*FR - CKI*FI + S1R
  192. S1I = CKR*FI + CKI*FR + S1I
  193. STR = PR - FR*AK
  194. STI = PI - FI*AK
  195. S2R = CKR*STR - CKI*STI + S2R
  196. S2I = CKR*STI + CKI*STR + S2I
  197. A1 = A1*T1*RAK
  198. BK = BK + AK + AK + 1.0D0
  199. AK = AK + 1.0D0
  200. IF (A1.GT.TOL) GO TO 90
  201. 100 CONTINUE
  202. KFLAG = 2
  203. A1 = FNU + 1.0D0
  204. AK = A1*ABS(SMUR)
  205. IF (AK.GT.ALIM) KFLAG = 3
  206. STR = CSSR(KFLAG)
  207. P2R = S2R*STR
  208. P2I = S2I*STR
  209. CALL ZMLT(P2R, P2I, RZR, RZI, S2R, S2I)
  210. S1R = S1R*STR
  211. S1I = S1I*STR
  212. IF (KODED.EQ.1) GO TO 210
  213. CALL ZEXP(ZR, ZI, FR, FI)
  214. CALL ZMLT(S1R, S1I, FR, FI, S1R, S1I)
  215. CALL ZMLT(S2R, S2I, FR, FI, S2R, S2I)
  216. GO TO 210
  217. C-----------------------------------------------------------------------
  218. C IFLAG=0 MEANS NO UNDERFLOW OCCURRED
  219. C IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH
  220. C KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD
  221. C RECURSION
  222. C-----------------------------------------------------------------------
  223. 110 CONTINUE
  224. CALL ZSQRT(ZR, ZI, STR, STI)
  225. CALL ZDIV(RTHPI, CZEROI, STR, STI, COEFR, COEFI)
  226. KFLAG = 2
  227. IF (KODED.EQ.2) GO TO 120
  228. IF (ZR.GT.ALIM) GO TO 290
  229. C BLANK LINE
  230. STR = EXP(-ZR)*CSSR(KFLAG)
  231. STI = -STR*SIN(ZI)
  232. STR = STR*COS(ZI)
  233. CALL ZMLT(COEFR, COEFI, STR, STI, COEFR, COEFI)
  234. 120 CONTINUE
  235. IF (ABS(DNU).EQ.0.5D0) GO TO 300
  236. C-----------------------------------------------------------------------
  237. C MILLER ALGORITHM FOR ABS(Z).GT.R1
  238. C-----------------------------------------------------------------------
  239. AK = COS(DPI*DNU)
  240. AK = ABS(AK)
  241. IF (AK.EQ.CZEROR) GO TO 300
  242. FHS = ABS(0.25D0-DNU2)
  243. IF (FHS.EQ.CZEROR) GO TO 300
  244. C-----------------------------------------------------------------------
  245. C COMPUTE R2=F(E). IF ABS(Z).GE.R2, USE FORWARD RECURRENCE TO
  246. C DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON
  247. C 12.LE.E.LE.60. E IS COMPUTED FROM 2**(-E)=B**(1-I1MACH(14))=
  248. C TOL WHERE B IS THE BASE OF THE ARITHMETIC.
  249. C-----------------------------------------------------------------------
  250. T1 = I1MACH(14)-1
  251. T1 = T1*D1MACH(5)*3.321928094D0
  252. T1 = MAX(T1,12.0D0)
  253. T1 = MIN(T1,60.0D0)
  254. T2 = TTH*T1 - 6.0D0
  255. IF (ZR.NE.0.0D0) GO TO 130
  256. T1 = HPI
  257. GO TO 140
  258. 130 CONTINUE
  259. T1 = DATAN(ZI/ZR)
  260. T1 = ABS(T1)
  261. 140 CONTINUE
  262. IF (T2.GT.CAZ) GO TO 170
  263. C-----------------------------------------------------------------------
  264. C FORWARD RECURRENCE LOOP WHEN ABS(Z).GE.R2
  265. C-----------------------------------------------------------------------
  266. ETEST = AK/(DPI*CAZ*TOL)
  267. FK = CONER
  268. IF (ETEST.LT.CONER) GO TO 180
  269. FKS = CTWOR
  270. CKR = CAZ + CAZ + CTWOR
  271. P1R = CZEROR
  272. P2R = CONER
  273. DO 150 I=1,KMAX
  274. AK = FHS/FKS
  275. CBR = CKR/(FK+CONER)
  276. PTR = P2R
  277. P2R = CBR*P2R - P1R*AK
  278. P1R = PTR
  279. CKR = CKR + CTWOR
  280. FKS = FKS + FK + FK + CTWOR
  281. FHS = FHS + FK + FK
  282. FK = FK + CONER
  283. STR = ABS(P2R)*FK
  284. IF (ETEST.LT.STR) GO TO 160
  285. 150 CONTINUE
  286. GO TO 310
  287. 160 CONTINUE
  288. FK = FK + SPI*T1*SQRT(T2/CAZ)
  289. FHS = ABS(0.25D0-DNU2)
  290. GO TO 180
  291. 170 CONTINUE
  292. C-----------------------------------------------------------------------
  293. C COMPUTE BACKWARD INDEX K FOR ABS(Z).LT.R2
  294. C-----------------------------------------------------------------------
  295. A2 = SQRT(CAZ)
  296. AK = FPI*AK/(TOL*SQRT(A2))
  297. AA = 3.0D0*T1/(1.0D0+CAZ)
  298. BB = 14.7D0*T1/(28.0D0+CAZ)
  299. AK = (LOG(AK)+CAZ*COS(AA)/(1.0D0+0.008D0*CAZ))/COS(BB)
  300. FK = 0.12125D0*AK*AK/CAZ + 1.5D0
  301. 180 CONTINUE
  302. C-----------------------------------------------------------------------
  303. C BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM
  304. C-----------------------------------------------------------------------
  305. K = FK
  306. FK = K
  307. FKS = FK*FK
  308. P1R = CZEROR
  309. P1I = CZEROI
  310. P2R = TOL
  311. P2I = CZEROI
  312. CSR = P2R
  313. CSI = P2I
  314. DO 190 I=1,K
  315. A1 = FKS - FK
  316. AK = (FKS+FK)/(A1+FHS)
  317. RAK = 2.0D0/(FK+CONER)
  318. CBR = (FK+ZR)*RAK
  319. CBI = ZI*RAK
  320. PTR = P2R
  321. PTI = P2I
  322. P2R = (PTR*CBR-PTI*CBI-P1R)*AK
  323. P2I = (PTI*CBR+PTR*CBI-P1I)*AK
  324. P1R = PTR
  325. P1I = PTI
  326. CSR = CSR + P2R
  327. CSI = CSI + P2I
  328. FKS = A1 - FK + CONER
  329. FK = FK - CONER
  330. 190 CONTINUE
  331. C-----------------------------------------------------------------------
  332. C COMPUTE (P2/CS)=(P2/ABS(CS))*(CONJG(CS)/ABS(CS)) FOR BETTER
  333. C SCALING
  334. C-----------------------------------------------------------------------
  335. TM = ZABS(CSR,CSI)
  336. PTR = 1.0D0/TM
  337. S1R = P2R*PTR
  338. S1I = P2I*PTR
  339. CSR = CSR*PTR
  340. CSI = -CSI*PTR
  341. CALL ZMLT(COEFR, COEFI, S1R, S1I, STR, STI)
  342. CALL ZMLT(STR, STI, CSR, CSI, S1R, S1I)
  343. IF (INU.GT.0 .OR. N.GT.1) GO TO 200
  344. ZDR = ZR
  345. ZDI = ZI
  346. IF(IFLAG.EQ.1) GO TO 270
  347. GO TO 240
  348. 200 CONTINUE
  349. C-----------------------------------------------------------------------
  350. C COMPUTE P1/P2=(P1/ABS(P2)*CONJG(P2)/ABS(P2) FOR SCALING
  351. C-----------------------------------------------------------------------
  352. TM = ZABS(P2R,P2I)
  353. PTR = 1.0D0/TM
  354. P1R = P1R*PTR
  355. P1I = P1I*PTR
  356. P2R = P2R*PTR
  357. P2I = -P2I*PTR
  358. CALL ZMLT(P1R, P1I, P2R, P2I, PTR, PTI)
  359. STR = DNU + 0.5D0 - PTR
  360. STI = -PTI
  361. CALL ZDIV(STR, STI, ZR, ZI, STR, STI)
  362. STR = STR + 1.0D0
  363. CALL ZMLT(STR, STI, S1R, S1I, S2R, S2I)
  364. C-----------------------------------------------------------------------
  365. C FORWARD RECURSION ON THE THREE TERM RECURSION WITH RELATION WITH
  366. C SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3
  367. C-----------------------------------------------------------------------
  368. 210 CONTINUE
  369. STR = DNU + 1.0D0
  370. CKR = STR*RZR
  371. CKI = STR*RZI
  372. IF (N.EQ.1) INU = INU - 1
  373. IF (INU.GT.0) GO TO 220
  374. IF (N.GT.1) GO TO 215
  375. S1R = S2R
  376. S1I = S2I
  377. 215 CONTINUE
  378. ZDR = ZR
  379. ZDI = ZI
  380. IF(IFLAG.EQ.1) GO TO 270
  381. GO TO 240
  382. 220 CONTINUE
  383. INUB = 1
  384. IF(IFLAG.EQ.1) GO TO 261
  385. 225 CONTINUE
  386. P1R = CSRR(KFLAG)
  387. ASCLE = BRY(KFLAG)
  388. DO 230 I=INUB,INU
  389. STR = S2R
  390. STI = S2I
  391. S2R = CKR*STR - CKI*STI + S1R
  392. S2I = CKR*STI + CKI*STR + S1I
  393. S1R = STR
  394. S1I = STI
  395. CKR = CKR + RZR
  396. CKI = CKI + RZI
  397. IF (KFLAG.GE.3) GO TO 230
  398. P2R = S2R*P1R
  399. P2I = S2I*P1R
  400. STR = ABS(P2R)
  401. STI = ABS(P2I)
  402. P2M = MAX(STR,STI)
  403. IF (P2M.LE.ASCLE) GO TO 230
  404. KFLAG = KFLAG + 1
  405. ASCLE = BRY(KFLAG)
  406. S1R = S1R*P1R
  407. S1I = S1I*P1R
  408. S2R = P2R
  409. S2I = P2I
  410. STR = CSSR(KFLAG)
  411. S1R = S1R*STR
  412. S1I = S1I*STR
  413. S2R = S2R*STR
  414. S2I = S2I*STR
  415. P1R = CSRR(KFLAG)
  416. 230 CONTINUE
  417. IF (N.NE.1) GO TO 240
  418. S1R = S2R
  419. S1I = S2I
  420. 240 CONTINUE
  421. STR = CSRR(KFLAG)
  422. YR(1) = S1R*STR
  423. YI(1) = S1I*STR
  424. IF (N.EQ.1) RETURN
  425. YR(2) = S2R*STR
  426. YI(2) = S2I*STR
  427. IF (N.EQ.2) RETURN
  428. KK = 2
  429. 250 CONTINUE
  430. KK = KK + 1
  431. IF (KK.GT.N) RETURN
  432. P1R = CSRR(KFLAG)
  433. ASCLE = BRY(KFLAG)
  434. DO 260 I=KK,N
  435. P2R = S2R
  436. P2I = S2I
  437. S2R = CKR*P2R - CKI*P2I + S1R
  438. S2I = CKI*P2R + CKR*P2I + S1I
  439. S1R = P2R
  440. S1I = P2I
  441. CKR = CKR + RZR
  442. CKI = CKI + RZI
  443. P2R = S2R*P1R
  444. P2I = S2I*P1R
  445. YR(I) = P2R
  446. YI(I) = P2I
  447. IF (KFLAG.GE.3) GO TO 260
  448. STR = ABS(P2R)
  449. STI = ABS(P2I)
  450. P2M = MAX(STR,STI)
  451. IF (P2M.LE.ASCLE) GO TO 260
  452. KFLAG = KFLAG + 1
  453. ASCLE = BRY(KFLAG)
  454. S1R = S1R*P1R
  455. S1I = S1I*P1R
  456. S2R = P2R
  457. S2I = P2I
  458. STR = CSSR(KFLAG)
  459. S1R = S1R*STR
  460. S1I = S1I*STR
  461. S2R = S2R*STR
  462. S2I = S2I*STR
  463. P1R = CSRR(KFLAG)
  464. 260 CONTINUE
  465. RETURN
  466. C-----------------------------------------------------------------------
  467. C IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW
  468. C-----------------------------------------------------------------------
  469. 261 CONTINUE
  470. HELIM = 0.5D0*ELIM
  471. ELM = EXP(-ELIM)
  472. CELMR = ELM
  473. ASCLE = BRY(1)
  474. ZDR = ZR
  475. ZDI = ZI
  476. IC = -1
  477. J = 2
  478. DO 262 I=1,INU
  479. STR = S2R
  480. STI = S2I
  481. S2R = STR*CKR-STI*CKI+S1R
  482. S2I = STI*CKR+STR*CKI+S1I
  483. S1R = STR
  484. S1I = STI
  485. CKR = CKR+RZR
  486. CKI = CKI+RZI
  487. AS = ZABS(S2R,S2I)
  488. ALAS = LOG(AS)
  489. P2R = -ZDR+ALAS
  490. IF(P2R.LT.(-ELIM)) GO TO 263
  491. CALL ZLOG(S2R,S2I,STR,STI,IDUM)
  492. P2R = -ZDR+STR
  493. P2I = -ZDI+STI
  494. P2M = EXP(P2R)/TOL
  495. P1R = P2M*COS(P2I)
  496. P1I = P2M*SIN(P2I)
  497. CALL ZUCHK(P1R,P1I,NW,ASCLE,TOL)
  498. IF(NW.NE.0) GO TO 263
  499. J = 3 - J
  500. CYR(J) = P1R
  501. CYI(J) = P1I
  502. IF(IC.EQ.(I-1)) GO TO 264
  503. IC = I
  504. GO TO 262
  505. 263 CONTINUE
  506. IF(ALAS.LT.HELIM) GO TO 262
  507. ZDR = ZDR-ELIM
  508. S1R = S1R*CELMR
  509. S1I = S1I*CELMR
  510. S2R = S2R*CELMR
  511. S2I = S2I*CELMR
  512. 262 CONTINUE
  513. IF(N.NE.1) GO TO 270
  514. S1R = S2R
  515. S1I = S2I
  516. GO TO 270
  517. 264 CONTINUE
  518. KFLAG = 1
  519. INUB = I+1
  520. S2R = CYR(J)
  521. S2I = CYI(J)
  522. J = 3 - J
  523. S1R = CYR(J)
  524. S1I = CYI(J)
  525. IF(INUB.LE.INU) GO TO 225
  526. IF(N.NE.1) GO TO 240
  527. S1R = S2R
  528. S1I = S2I
  529. GO TO 240
  530. 270 CONTINUE
  531. YR(1) = S1R
  532. YI(1) = S1I
  533. IF(N.EQ.1) GO TO 280
  534. YR(2) = S2R
  535. YI(2) = S2I
  536. 280 CONTINUE
  537. ASCLE = BRY(1)
  538. CALL ZKSCL(ZDR,ZDI,FNU,N,YR,YI,NZ,RZR,RZI,ASCLE,TOL,ELIM)
  539. INU = N - NZ
  540. IF (INU.LE.0) RETURN
  541. KK = NZ + 1
  542. S1R = YR(KK)
  543. S1I = YI(KK)
  544. YR(KK) = S1R*CSRR(1)
  545. YI(KK) = S1I*CSRR(1)
  546. IF (INU.EQ.1) RETURN
  547. KK = NZ + 2
  548. S2R = YR(KK)
  549. S2I = YI(KK)
  550. YR(KK) = S2R*CSRR(1)
  551. YI(KK) = S2I*CSRR(1)
  552. IF (INU.EQ.2) RETURN
  553. T2 = FNU + (KK-1)
  554. CKR = T2*RZR
  555. CKI = T2*RZI
  556. KFLAG = 1
  557. GO TO 250
  558. 290 CONTINUE
  559. C-----------------------------------------------------------------------
  560. C SCALE BY EXP(Z), IFLAG = 1 CASES
  561. C-----------------------------------------------------------------------
  562. KODED = 2
  563. IFLAG = 1
  564. KFLAG = 2
  565. GO TO 120
  566. C-----------------------------------------------------------------------
  567. C FNU=HALF ODD INTEGER CASE, DNU=-0.5
  568. C-----------------------------------------------------------------------
  569. 300 CONTINUE
  570. S1R = COEFR
  571. S1I = COEFI
  572. S2R = COEFR
  573. S2I = COEFI
  574. GO TO 210
  575. C
  576. C
  577. 310 CONTINUE
  578. NZ=-2
  579. RETURN
  580. END