asyjy.f 22 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491
  1. *DECK ASYJY
  2. SUBROUTINE ASYJY (FUNJY, X, FNU, FLGJY, IN, Y, WK, IFLW)
  3. C***BEGIN PROLOGUE ASYJY
  4. C***SUBSIDIARY
  5. C***PURPOSE Subsidiary to BESJ and BESY
  6. C***LIBRARY SLATEC
  7. C***TYPE SINGLE PRECISION (ASYJY-S, DASYJY-D)
  8. C***AUTHOR Amos, D. E., (SNLA)
  9. C***DESCRIPTION
  10. C
  11. C ASYJY computes Bessel functions J and Y
  12. C for arguments X.GT.0.0 and orders FNU.GE.35.0
  13. C on FLGJY = 1 and FLGJY = -1 respectively
  14. C
  15. C INPUT
  16. C
  17. C FUNJY - external function JAIRY or YAIRY
  18. C X - argument, X.GT.0.0E0
  19. C FNU - order of the first Bessel function
  20. C FLGJY - selection flag
  21. C FLGJY = 1.0E0 gives the J function
  22. C FLGJY = -1.0E0 gives the Y function
  23. C IN - number of functions desired, IN = 1 or 2
  24. C
  25. C OUTPUT
  26. C
  27. C Y - a vector whose first in components contain the sequence
  28. C IFLW - a flag indicating underflow or overflow
  29. C return variables for BESJ only
  30. C WK(1) = 1 - (X/FNU)**2 = W**2
  31. C WK(2) = SQRT(ABS(WK(1)))
  32. C WK(3) = ABS(WK(2) - ATAN(WK(2))) or
  33. C ABS(LN((1 + WK(2))/(X/FNU)) - WK(2))
  34. C = ABS((2/3)*ZETA**(3/2))
  35. C WK(4) = FNU*WK(3)
  36. C WK(5) = (1.5*WK(3)*FNU)**(1/3) = SQRT(ZETA)*FNU**(1/3)
  37. C WK(6) = SIGN(1.,W**2)*WK(5)**2 = SIGN(1.,W**2)*ZETA*FNU**(2/3)
  38. C WK(7) = FNU**(1/3)
  39. C
  40. C Abstract
  41. C ASYJY implements the uniform asymptotic expansion of
  42. C the J and Y Bessel functions for FNU.GE.35 and real
  43. C X.GT.0.0E0. The forms are identical except for a change
  44. C in sign of some of the terms. This change in sign is
  45. C accomplished by means of the flag FLGJY = 1 or -1. On
  46. C FLGJY = 1 the AIRY functions AI(X) and DAI(X) are
  47. C supplied by the external function JAIRY, and on
  48. C FLGJY = -1 the AIRY functions BI(X) and DBI(X) are
  49. C supplied by the external function YAIRY.
  50. C
  51. C***SEE ALSO BESJ, BESY
  52. C***ROUTINES CALLED I1MACH, R1MACH
  53. C***REVISION HISTORY (YYMMDD)
  54. C 750101 DATE WRITTEN
  55. C 890531 Changed all specific intrinsics to generic. (WRB)
  56. C 891009 Removed unreferenced variable. (WRB)
  57. C 891214 Prologue converted to Version 4.0 format. (BAB)
  58. C 900328 Added TYPE section. (WRB)
  59. C 910408 Updated the AUTHOR section. (WRB)
  60. C***END PROLOGUE ASYJY
  61. INTEGER I, IFLW, IN, J, JN,JR,JU,K, KB,KLAST,KMAX,KP1, KS, KSP1,
  62. * KSTEMP, L, LR, LRP1, ISETA, ISETB
  63. INTEGER I1MACH
  64. REAL ABW2, AKM, ALFA, ALFA1, ALFA2, AP, AR, ASUM, AZ,
  65. * BETA, BETA1, BETA2, BETA3, BR, BSUM, C, CON1, CON2,
  66. * CON548,CR,CRZ32, DFI,ELIM, DR,FI, FLGJY, FN, FNU,
  67. * FN2, GAMA, PHI, RCZ, RDEN, RELB, RFN2, RTZ, RZDEN,
  68. * SA, SB, SUMA, SUMB, S1, TA, TAU, TB, TFN, TOL, TOLS, T2, UPOL,
  69. * WK, X, XX, Y, Z, Z32
  70. REAL R1MACH
  71. DIMENSION Y(*), WK(*), C(65)
  72. DIMENSION ALFA(26,4), BETA(26,5)
  73. DIMENSION ALFA1(26,2), ALFA2(26,2)
  74. DIMENSION BETA1(26,2), BETA2(26,2), BETA3(26,1)
  75. DIMENSION GAMA(26), KMAX(5), AR(8), BR(10), UPOL(10)
  76. DIMENSION CR(10), DR(10)
  77. EQUIVALENCE (ALFA(1,1),ALFA1(1,1))
  78. EQUIVALENCE (ALFA(1,3),ALFA2(1,1))
  79. EQUIVALENCE (BETA(1,1),BETA1(1,1))
  80. EQUIVALENCE (BETA(1,3),BETA2(1,1))
  81. EQUIVALENCE (BETA(1,5),BETA3(1,1))
  82. SAVE TOLS, CON1, CON2, CON548, AR, BR, C, ALFA1, ALFA2,
  83. 1 BETA1, BETA2, BETA3, GAMA
  84. DATA TOLS /-6.90775527898214E+00/
  85. DATA CON1,CON2,CON548/
  86. 1 6.66666666666667E-01, 3.33333333333333E-01, 1.04166666666667E-01/
  87. DATA AR(1), AR(2), AR(3), AR(4), AR(5), AR(6), AR(7),
  88. A AR(8) / 8.35503472222222E-02, 1.28226574556327E-01,
  89. 1 2.91849026464140E-01, 8.81627267443758E-01, 3.32140828186277E+00,
  90. 2 1.49957629868626E+01, 7.89230130115865E+01, 4.74451538868264E+02/
  91. DATA BR(1), BR(2), BR(3), BR(4), BR(5), BR(6), BR(7), BR(8),
  92. A BR(9), BR(10) /-1.45833333333333E-01,-9.87413194444444E-02,
  93. 1-1.43312053915895E-01,-3.17227202678414E-01,-9.42429147957120E-01,
  94. 2-3.51120304082635E+00,-1.57272636203680E+01,-8.22814390971859E+01,
  95. 3-4.92355370523671E+02,-3.31621856854797E+03/
  96. DATA C(1), C(2), C(3), C(4), C(5), C(6), C(7), C(8), C(9), C(10),
  97. 1 C(11), C(12), C(13), C(14), C(15), C(16), C(17), C(18),
  98. 2 C(19), C(20), C(21), C(22), C(23), C(24)/
  99. 3 -2.08333333333333E-01, 1.25000000000000E-01,
  100. 4 3.34201388888889E-01, -4.01041666666667E-01,
  101. 5 7.03125000000000E-02, -1.02581259645062E+00,
  102. 6 1.84646267361111E+00, -8.91210937500000E-01,
  103. 7 7.32421875000000E-02, 4.66958442342625E+00,
  104. 8 -1.12070026162230E+01, 8.78912353515625E+00,
  105. 9 -2.36408691406250E+00, 1.12152099609375E-01,
  106. A -2.82120725582002E+01, 8.46362176746007E+01,
  107. B -9.18182415432400E+01, 4.25349987453885E+01,
  108. C -7.36879435947963E+00, 2.27108001708984E-01,
  109. D 2.12570130039217E+02, -7.65252468141182E+02,
  110. E 1.05999045252800E+03, -6.99579627376133E+02/
  111. DATA C(25), C(26), C(27), C(28), C(29), C(30), C(31), C(32),
  112. 1 C(33), C(34), C(35), C(36), C(37), C(38), C(39), C(40),
  113. 2 C(41), C(42), C(43), C(44), C(45), C(46), C(47), C(48)/
  114. 3 2.18190511744212E+02, -2.64914304869516E+01,
  115. 4 5.72501420974731E-01, -1.91945766231841E+03,
  116. 5 8.06172218173731E+03, -1.35865500064341E+04,
  117. 6 1.16553933368645E+04, -5.30564697861340E+03,
  118. 7 1.20090291321635E+03, -1.08090919788395E+02,
  119. 8 1.72772750258446E+00, 2.02042913309661E+04,
  120. 9 -9.69805983886375E+04, 1.92547001232532E+05,
  121. A -2.03400177280416E+05, 1.22200464983017E+05,
  122. B -4.11926549688976E+04, 7.10951430248936E+03,
  123. C -4.93915304773088E+02, 6.07404200127348E+00,
  124. D -2.42919187900551E+05, 1.31176361466298E+06,
  125. E -2.99801591853811E+06, 3.76327129765640E+06/
  126. DATA C(49), C(50), C(51), C(52), C(53), C(54), C(55), C(56),
  127. 1 C(57), C(58), C(59), C(60), C(61), C(62), C(63), C(64),
  128. 2 C(65)/
  129. 3 -2.81356322658653E+06, 1.26836527332162E+06,
  130. 4 -3.31645172484564E+05, 4.52187689813627E+04,
  131. 5 -2.49983048181121E+03, 2.43805296995561E+01,
  132. 6 3.28446985307204E+06, -1.97068191184322E+07,
  133. 7 5.09526024926646E+07, -7.41051482115327E+07,
  134. 8 6.63445122747290E+07, -3.75671766607634E+07,
  135. 9 1.32887671664218E+07, -2.78561812808645E+06,
  136. A 3.08186404612662E+05, -1.38860897537170E+04,
  137. B 1.10017140269247E+02/
  138. DATA ALFA1(1,1), ALFA1(2,1), ALFA1(3,1), ALFA1(4,1), ALFA1(5,1),
  139. 1 ALFA1(6,1), ALFA1(7,1), ALFA1(8,1), ALFA1(9,1), ALFA1(10,1),
  140. 2 ALFA1(11,1),ALFA1(12,1),ALFA1(13,1),ALFA1(14,1),ALFA1(15,1),
  141. 3 ALFA1(16,1),ALFA1(17,1),ALFA1(18,1),ALFA1(19,1),ALFA1(20,1),
  142. 4 ALFA1(21,1),ALFA1(22,1),ALFA1(23,1),ALFA1(24,1),ALFA1(25,1),
  143. 5 ALFA1(26,1) /-4.44444444444444E-03,-9.22077922077922E-04,
  144. 6-8.84892884892885E-05, 1.65927687832450E-04, 2.46691372741793E-04,
  145. 7 2.65995589346255E-04, 2.61824297061501E-04, 2.48730437344656E-04,
  146. 8 2.32721040083232E-04, 2.16362485712365E-04, 2.00738858762752E-04,
  147. 9 1.86267636637545E-04, 1.73060775917876E-04, 1.61091705929016E-04,
  148. 1 1.50274774160908E-04, 1.40503497391270E-04, 1.31668816545923E-04,
  149. 2 1.23667445598253E-04, 1.16405271474738E-04, 1.09798298372713E-04,
  150. 3 1.03772410422993E-04, 9.82626078369363E-05, 9.32120517249503E-05,
  151. 4 8.85710852478712E-05, 8.42963105715700E-05, 8.03497548407791E-05/
  152. DATA ALFA1(1,2), ALFA1(2,2), ALFA1(3,2), ALFA1(4,2), ALFA1(5,2),
  153. 1 ALFA1(6,2), ALFA1(7,2), ALFA1(8,2), ALFA1(9,2), ALFA1(10,2),
  154. 2 ALFA1(11,2),ALFA1(12,2),ALFA1(13,2),ALFA1(14,2),ALFA1(15,2),
  155. 3 ALFA1(16,2),ALFA1(17,2),ALFA1(18,2),ALFA1(19,2),ALFA1(20,2),
  156. 4 ALFA1(21,2),ALFA1(22,2),ALFA1(23,2),ALFA1(24,2),ALFA1(25,2),
  157. 5 ALFA1(26,2) / 6.93735541354589E-04, 2.32241745182922E-04,
  158. 6-1.41986273556691E-05,-1.16444931672049E-04,-1.50803558053049E-04,
  159. 7-1.55121924918096E-04,-1.46809756646466E-04,-1.33815503867491E-04,
  160. 8-1.19744975684254E-04,-1.06184319207974E-04,-9.37699549891194E-05,
  161. 9-8.26923045588193E-05,-7.29374348155221E-05,-6.44042357721016E-05,
  162. 1-5.69611566009369E-05,-5.04731044303562E-05,-4.48134868008883E-05,
  163. 2-3.98688727717599E-05,-3.55400532972042E-05,-3.17414256609022E-05,
  164. 3-2.83996793904175E-05,-2.54522720634871E-05,-2.28459297164725E-05,
  165. 4-2.05352753106481E-05,-1.84816217627666E-05,-1.66519330021394E-05/
  166. DATA ALFA2(1,1), ALFA2(2,1), ALFA2(3,1), ALFA2(4,1), ALFA2(5,1),
  167. 1 ALFA2(6,1), ALFA2(7,1), ALFA2(8,1), ALFA2(9,1), ALFA2(10,1),
  168. 2 ALFA2(11,1),ALFA2(12,1),ALFA2(13,1),ALFA2(14,1),ALFA2(15,1),
  169. 3 ALFA2(16,1),ALFA2(17,1),ALFA2(18,1),ALFA2(19,1),ALFA2(20,1),
  170. 4 ALFA2(21,1),ALFA2(22,1),ALFA2(23,1),ALFA2(24,1),ALFA2(25,1),
  171. 5 ALFA2(26,1) /-3.54211971457744E-04,-1.56161263945159E-04,
  172. 6 3.04465503594936E-05, 1.30198655773243E-04, 1.67471106699712E-04,
  173. 7 1.70222587683593E-04, 1.56501427608595E-04, 1.36339170977445E-04,
  174. 8 1.14886692029825E-04, 9.45869093034688E-05, 7.64498419250898E-05,
  175. 9 6.07570334965197E-05, 4.74394299290509E-05, 3.62757512005344E-05,
  176. 1 2.69939714979225E-05, 1.93210938247939E-05, 1.30056674793963E-05,
  177. 2 7.82620866744497E-06, 3.59257485819352E-06, 1.44040049814252E-07,
  178. 3-2.65396769697939E-06,-4.91346867098486E-06,-6.72739296091248E-06,
  179. 4-8.17269379678658E-06,-9.31304715093561E-06,-1.02011418798016E-05/
  180. DATA ALFA2(1,2), ALFA2(2,2), ALFA2(3,2), ALFA2(4,2), ALFA2(5,2),
  181. 1 ALFA2(6,2), ALFA2(7,2), ALFA2(8,2), ALFA2(9,2), ALFA2(10,2),
  182. 2 ALFA2(11,2),ALFA2(12,2),ALFA2(13,2),ALFA2(14,2),ALFA2(15,2),
  183. 3 ALFA2(16,2),ALFA2(17,2),ALFA2(18,2),ALFA2(19,2),ALFA2(20,2),
  184. 4 ALFA2(21,2),ALFA2(22,2),ALFA2(23,2),ALFA2(24,2),ALFA2(25,2),
  185. 5 ALFA2(26,2) / 3.78194199201773E-04, 2.02471952761816E-04,
  186. 6-6.37938506318862E-05,-2.38598230603006E-04,-3.10916256027362E-04,
  187. 7-3.13680115247576E-04,-2.78950273791323E-04,-2.28564082619141E-04,
  188. 8-1.75245280340847E-04,-1.25544063060690E-04,-8.22982872820208E-05,
  189. 9-4.62860730588116E-05,-1.72334302366962E-05, 5.60690482304602E-06,
  190. 1 2.31395443148287E-05, 3.62642745856794E-05, 4.58006124490189E-05,
  191. 2 5.24595294959114E-05, 5.68396208545815E-05, 5.94349820393104E-05,
  192. 3 6.06478527578422E-05, 6.08023907788436E-05, 6.01577894539460E-05,
  193. 4 5.89199657344698E-05, 5.72515823777593E-05, 5.52804375585853E-05/
  194. DATA BETA1(1,1), BETA1(2,1), BETA1(3,1), BETA1(4,1), BETA1(5,1),
  195. 1 BETA1(6,1), BETA1(7,1), BETA1(8,1), BETA1(9,1), BETA1(10,1),
  196. 2 BETA1(11,1),BETA1(12,1),BETA1(13,1),BETA1(14,1),BETA1(15,1),
  197. 3 BETA1(16,1),BETA1(17,1),BETA1(18,1),BETA1(19,1),BETA1(20,1),
  198. 4 BETA1(21,1),BETA1(22,1),BETA1(23,1),BETA1(24,1),BETA1(25,1),
  199. 5 BETA1(26,1) / 1.79988721413553E-02, 5.59964911064388E-03,
  200. 6 2.88501402231133E-03, 1.80096606761054E-03, 1.24753110589199E-03,
  201. 7 9.22878876572938E-04, 7.14430421727287E-04, 5.71787281789705E-04,
  202. 8 4.69431007606482E-04, 3.93232835462917E-04, 3.34818889318298E-04,
  203. 9 2.88952148495752E-04, 2.52211615549573E-04, 2.22280580798883E-04,
  204. 1 1.97541838033063E-04, 1.76836855019718E-04, 1.59316899661821E-04,
  205. 2 1.44347930197334E-04, 1.31448068119965E-04, 1.20245444949303E-04,
  206. 3 1.10449144504599E-04, 1.01828770740567E-04, 9.41998224204238E-05,
  207. 4 8.74130545753834E-05, 8.13466262162801E-05, 7.59002269646219E-05/
  208. DATA BETA1(1,2), BETA1(2,2), BETA1(3,2), BETA1(4,2), BETA1(5,2),
  209. 1 BETA1(6,2), BETA1(7,2), BETA1(8,2), BETA1(9,2), BETA1(10,2),
  210. 2 BETA1(11,2),BETA1(12,2),BETA1(13,2),BETA1(14,2),BETA1(15,2),
  211. 3 BETA1(16,2),BETA1(17,2),BETA1(18,2),BETA1(19,2),BETA1(20,2),
  212. 4 BETA1(21,2),BETA1(22,2),BETA1(23,2),BETA1(24,2),BETA1(25,2),
  213. 5 BETA1(26,2) /-1.49282953213429E-03,-8.78204709546389E-04,
  214. 6-5.02916549572035E-04,-2.94822138512746E-04,-1.75463996970783E-04,
  215. 7-1.04008550460816E-04,-5.96141953046458E-05,-3.12038929076098E-05,
  216. 8-1.26089735980230E-05,-2.42892608575730E-07, 8.05996165414274E-06,
  217. 9 1.36507009262147E-05, 1.73964125472926E-05, 1.98672978842134E-05,
  218. 1 2.14463263790823E-05, 2.23954659232457E-05, 2.28967783814713E-05,
  219. 2 2.30785389811178E-05, 2.30321976080909E-05, 2.28236073720349E-05,
  220. 3 2.25005881105292E-05, 2.20981015361991E-05, 2.16418427448104E-05,
  221. 4 2.11507649256221E-05, 2.06388749782171E-05, 2.01165241997082E-05/
  222. DATA BETA2(1,1), BETA2(2,1), BETA2(3,1), BETA2(4,1), BETA2(5,1),
  223. 1 BETA2(6,1), BETA2(7,1), BETA2(8,1), BETA2(9,1), BETA2(10,1),
  224. 2 BETA2(11,1),BETA2(12,1),BETA2(13,1),BETA2(14,1),BETA2(15,1),
  225. 3 BETA2(16,1),BETA2(17,1),BETA2(18,1),BETA2(19,1),BETA2(20,1),
  226. 4 BETA2(21,1),BETA2(22,1),BETA2(23,1),BETA2(24,1),BETA2(25,1),
  227. 5 BETA2(26,1) / 5.52213076721293E-04, 4.47932581552385E-04,
  228. 6 2.79520653992021E-04, 1.52468156198447E-04, 6.93271105657044E-05,
  229. 7 1.76258683069991E-05,-1.35744996343269E-05,-3.17972413350427E-05,
  230. 8-4.18861861696693E-05,-4.69004889379141E-05,-4.87665447413787E-05,
  231. 9-4.87010031186735E-05,-4.74755620890087E-05,-4.55813058138628E-05,
  232. 1-4.33309644511266E-05,-4.09230193157750E-05,-3.84822638603221E-05,
  233. 2-3.60857167535411E-05,-3.37793306123367E-05,-3.15888560772110E-05,
  234. 3-2.95269561750807E-05,-2.75978914828336E-05,-2.58006174666884E-05,
  235. 4-2.41308356761280E-05,-2.25823509518346E-05,-2.11479656768913E-05/
  236. DATA BETA2(1,2), BETA2(2,2), BETA2(3,2), BETA2(4,2), BETA2(5,2),
  237. 1 BETA2(6,2), BETA2(7,2), BETA2(8,2), BETA2(9,2), BETA2(10,2),
  238. 2 BETA2(11,2),BETA2(12,2),BETA2(13,2),BETA2(14,2),BETA2(15,2),
  239. 3 BETA2(16,2),BETA2(17,2),BETA2(18,2),BETA2(19,2),BETA2(20,2),
  240. 4 BETA2(21,2),BETA2(22,2),BETA2(23,2),BETA2(24,2),BETA2(25,2),
  241. 5 BETA2(26,2) /-4.74617796559960E-04,-4.77864567147321E-04,
  242. 6-3.20390228067038E-04,-1.61105016119962E-04,-4.25778101285435E-05,
  243. 7 3.44571294294968E-05, 7.97092684075675E-05, 1.03138236708272E-04,
  244. 8 1.12466775262204E-04, 1.13103642108481E-04, 1.08651634848774E-04,
  245. 9 1.01437951597662E-04, 9.29298396593364E-05, 8.40293133016090E-05,
  246. 1 7.52727991349134E-05, 6.69632521975731E-05, 5.92564547323195E-05,
  247. 2 5.22169308826976E-05, 4.58539485165361E-05, 4.01445513891487E-05,
  248. 3 3.50481730031328E-05, 3.05157995034347E-05, 2.64956119950516E-05,
  249. 4 2.29363633690998E-05, 1.97893056664022E-05, 1.70091984636413E-05/
  250. DATA BETA3(1,1), BETA3(2,1), BETA3(3,1), BETA3(4,1), BETA3(5,1),
  251. 1 BETA3(6,1), BETA3(7,1), BETA3(8,1), BETA3(9,1), BETA3(10,1),
  252. 2 BETA3(11,1),BETA3(12,1),BETA3(13,1),BETA3(14,1),BETA3(15,1),
  253. 3 BETA3(16,1),BETA3(17,1),BETA3(18,1),BETA3(19,1),BETA3(20,1),
  254. 4 BETA3(21,1),BETA3(22,1),BETA3(23,1),BETA3(24,1),BETA3(25,1),
  255. 5 BETA3(26,1) / 7.36465810572578E-04, 8.72790805146194E-04,
  256. 6 6.22614862573135E-04, 2.85998154194304E-04, 3.84737672879366E-06,
  257. 7-1.87906003636972E-04,-2.97603646594555E-04,-3.45998126832656E-04,
  258. 8-3.53382470916038E-04,-3.35715635775049E-04,-3.04321124789040E-04,
  259. 9-2.66722723047613E-04,-2.27654214122820E-04,-1.89922611854562E-04,
  260. 1-1.55058918599094E-04,-1.23778240761874E-04,-9.62926147717644E-05,
  261. 2-7.25178327714425E-05,-5.22070028895634E-05,-3.50347750511901E-05,
  262. 3-2.06489761035552E-05,-8.70106096849767E-06, 1.13698686675100E-06,
  263. 4 9.16426474122779E-06, 1.56477785428873E-05, 2.08223629482467E-05/
  264. DATA GAMA(1), GAMA(2), GAMA(3), GAMA(4), GAMA(5),
  265. 1 GAMA(6), GAMA(7), GAMA(8), GAMA(9), GAMA(10),
  266. 2 GAMA(11), GAMA(12), GAMA(13), GAMA(14), GAMA(15),
  267. 3 GAMA(16), GAMA(17), GAMA(18), GAMA(19), GAMA(20),
  268. 4 GAMA(21), GAMA(22), GAMA(23), GAMA(24), GAMA(25),
  269. 5 GAMA(26) / 6.29960524947437E-01, 2.51984209978975E-01,
  270. 6 1.54790300415656E-01, 1.10713062416159E-01, 8.57309395527395E-02,
  271. 7 6.97161316958684E-02, 5.86085671893714E-02, 5.04698873536311E-02,
  272. 8 4.42600580689155E-02, 3.93720661543510E-02, 3.54283195924455E-02,
  273. 9 3.21818857502098E-02, 2.94646240791158E-02, 2.71581677112934E-02,
  274. 1 2.51768272973862E-02, 2.34570755306079E-02, 2.19508390134907E-02,
  275. 2 2.06210828235646E-02, 1.94388240897881E-02, 1.83810633800683E-02,
  276. 3 1.74293213231963E-02, 1.65685837786612E-02, 1.57865285987918E-02,
  277. 4 1.50729501494096E-02, 1.44193250839955E-02, 1.38184805735342E-02/
  278. C***FIRST EXECUTABLE STATEMENT ASYJY
  279. TA = R1MACH(3)
  280. TOL = MAX(TA,1.0E-15)
  281. TB = R1MACH(5)
  282. JU = I1MACH(12)
  283. IF(FLGJY.EQ.1.0E0) GO TO 6
  284. JR = I1MACH(11)
  285. ELIM = -2.303E0*TB*(JU+JR)
  286. GO TO 7
  287. 6 CONTINUE
  288. ELIM = -2.303E0*(TB*JU+3.0E0)
  289. 7 CONTINUE
  290. FN = FNU
  291. IFLW = 0
  292. DO 170 JN=1,IN
  293. XX = X/FN
  294. WK(1) = 1.0E0 - XX*XX
  295. ABW2 = ABS(WK(1))
  296. WK(2) = SQRT(ABW2)
  297. WK(7) = FN**CON2
  298. IF (ABW2.GT.0.27750E0) GO TO 80
  299. C
  300. C ASYMPTOTIC EXPANSION
  301. C CASES NEAR X=FN, ABS(1.-(X/FN)**2).LE.0.2775
  302. C COEFFICIENTS OF ASYMPTOTIC EXPANSION BY SERIES
  303. C
  304. C ZETA AND TRUNCATION FOR A(ZETA) AND B(ZETA) SERIES
  305. C
  306. C KMAX IS TRUNCATION INDEX FOR A(ZETA) AND B(ZETA) SERIES=MAX(2,SA)
  307. C
  308. SA = 0.0E0
  309. IF (ABW2.EQ.0.0E0) GO TO 10
  310. SA = TOLS/LOG(ABW2)
  311. 10 SB = SA
  312. DO 20 I=1,5
  313. AKM = MAX(SA,2.0E0)
  314. KMAX(I) = INT(AKM)
  315. SA = SA + SB
  316. 20 CONTINUE
  317. KB = KMAX(5)
  318. KLAST = KB - 1
  319. SA = GAMA(KB)
  320. DO 30 K=1,KLAST
  321. KB = KB - 1
  322. SA = SA*WK(1) + GAMA(KB)
  323. 30 CONTINUE
  324. Z = WK(1)*SA
  325. AZ = ABS(Z)
  326. RTZ = SQRT(AZ)
  327. WK(3) = CON1*AZ*RTZ
  328. WK(4) = WK(3)*FN
  329. WK(5) = RTZ*WK(7)
  330. WK(6) = -WK(5)*WK(5)
  331. IF(Z.LE.0.0E0) GO TO 35
  332. IF(WK(4).GT.ELIM) GO TO 75
  333. WK(6) = -WK(6)
  334. 35 CONTINUE
  335. PHI = SQRT(SQRT(SA+SA+SA+SA))
  336. C
  337. C B(ZETA) FOR S=0
  338. C
  339. KB = KMAX(5)
  340. KLAST = KB - 1
  341. SB = BETA(KB,1)
  342. DO 40 K=1,KLAST
  343. KB = KB - 1
  344. SB = SB*WK(1) + BETA(KB,1)
  345. 40 CONTINUE
  346. KSP1 = 1
  347. FN2 = FN*FN
  348. RFN2 = 1.0E0/FN2
  349. RDEN = 1.0E0
  350. ASUM = 1.0E0
  351. RELB = TOL*ABS(SB)
  352. BSUM = SB
  353. DO 60 KS=1,4
  354. KSP1 = KSP1 + 1
  355. RDEN = RDEN*RFN2
  356. C
  357. C A(ZETA) AND B(ZETA) FOR S=1,2,3,4
  358. C
  359. KSTEMP = 5 - KS
  360. KB = KMAX(KSTEMP)
  361. KLAST = KB - 1
  362. SA = ALFA(KB,KS)
  363. SB = BETA(KB,KSP1)
  364. DO 50 K=1,KLAST
  365. KB = KB - 1
  366. SA = SA*WK(1) + ALFA(KB,KS)
  367. SB = SB*WK(1) + BETA(KB,KSP1)
  368. 50 CONTINUE
  369. TA = SA*RDEN
  370. TB = SB*RDEN
  371. ASUM = ASUM + TA
  372. BSUM = BSUM + TB
  373. IF (ABS(TA).LE.TOL .AND. ABS(TB).LE.RELB) GO TO 70
  374. 60 CONTINUE
  375. 70 CONTINUE
  376. BSUM = BSUM/(FN*WK(7))
  377. GO TO 160
  378. C
  379. 75 CONTINUE
  380. IFLW = 1
  381. RETURN
  382. C
  383. 80 CONTINUE
  384. UPOL(1) = 1.0E0
  385. TAU = 1.0E0/WK(2)
  386. T2 = 1.0E0/WK(1)
  387. IF (WK(1).GE.0.0E0) GO TO 90
  388. C
  389. C CASES FOR (X/FN).GT.SQRT(1.2775)
  390. C
  391. WK(3) = ABS(WK(2)-ATAN(WK(2)))
  392. WK(4) = WK(3)*FN
  393. RCZ = -CON1/WK(4)
  394. Z32 = 1.5E0*WK(3)
  395. RTZ = Z32**CON2
  396. WK(5) = RTZ*WK(7)
  397. WK(6) = -WK(5)*WK(5)
  398. GO TO 100
  399. 90 CONTINUE
  400. C
  401. C CASES FOR (X/FN).LT.SQRT(0.7225)
  402. C
  403. WK(3) = ABS(LOG((1.0E0+WK(2))/XX)-WK(2))
  404. WK(4) = WK(3)*FN
  405. RCZ = CON1/WK(4)
  406. IF(WK(4).GT.ELIM) GO TO 75
  407. Z32 = 1.5E0*WK(3)
  408. RTZ = Z32**CON2
  409. WK(7) = FN**CON2
  410. WK(5) = RTZ*WK(7)
  411. WK(6) = WK(5)*WK(5)
  412. 100 CONTINUE
  413. PHI = SQRT((RTZ+RTZ)*TAU)
  414. TB = 1.0E0
  415. ASUM = 1.0E0
  416. TFN = TAU/FN
  417. RDEN=1.0E0/FN
  418. RFN2=RDEN*RDEN
  419. RDEN=1.0E0
  420. UPOL(2) = (C(1)*T2+C(2))*TFN
  421. CRZ32 = CON548*RCZ
  422. BSUM = UPOL(2) + CRZ32
  423. RELB = TOL*ABS(BSUM)
  424. AP = TFN
  425. KS = 0
  426. KP1 = 2
  427. RZDEN = RCZ
  428. L = 2
  429. ISETA=0
  430. ISETB=0
  431. DO 140 LR=2,8,2
  432. C
  433. C COMPUTE TWO U POLYNOMIALS FOR NEXT A(ZETA) AND B(ZETA)
  434. C
  435. LRP1 = LR + 1
  436. DO 120 K=LR,LRP1
  437. KS = KS + 1
  438. KP1 = KP1 + 1
  439. L = L + 1
  440. S1 = C(L)
  441. DO 110 J=2,KP1
  442. L = L + 1
  443. S1 = S1*T2 + C(L)
  444. 110 CONTINUE
  445. AP = AP*TFN
  446. UPOL(KP1) = AP*S1
  447. CR(KS) = BR(KS)*RZDEN
  448. RZDEN = RZDEN*RCZ
  449. DR(KS) = AR(KS)*RZDEN
  450. 120 CONTINUE
  451. SUMA = UPOL(LRP1)
  452. SUMB = UPOL(LR+2) + UPOL(LRP1)*CRZ32
  453. JU = LRP1
  454. DO 130 JR=1,LR
  455. JU = JU - 1
  456. SUMA = SUMA + CR(JR)*UPOL(JU)
  457. SUMB = SUMB + DR(JR)*UPOL(JU)
  458. 130 CONTINUE
  459. RDEN=RDEN*RFN2
  460. TB = -TB
  461. IF (WK(1).GT.0.0E0) TB = ABS(TB)
  462. IF (RDEN.LT.TOL) GO TO 131
  463. ASUM = ASUM + SUMA*TB
  464. BSUM = BSUM + SUMB*TB
  465. GO TO 140
  466. 131 IF(ISETA.EQ.1) GO TO 132
  467. IF(ABS(SUMA).LT.TOL) ISETA=1
  468. ASUM=ASUM+SUMA*TB
  469. 132 IF(ISETB.EQ.1) GO TO 133
  470. IF(ABS(SUMB).LT.RELB) ISETB=1
  471. BSUM=BSUM+SUMB*TB
  472. 133 IF(ISETA.EQ.1 .AND. ISETB.EQ.1) GO TO 150
  473. 140 CONTINUE
  474. 150 TB = WK(5)
  475. IF (WK(1).GT.0.0E0) TB = -TB
  476. BSUM = BSUM/TB
  477. C
  478. 160 CONTINUE
  479. CALL FUNJY(WK(6), WK(5), WK(4), FI, DFI)
  480. TA=1.0E0/TOL
  481. TB=R1MACH(1)*TA*1.0E+3
  482. IF(ABS(FI).GT.TB) GO TO 165
  483. FI=FI*TA
  484. DFI=DFI*TA
  485. PHI=PHI*TOL
  486. 165 CONTINUE
  487. Y(JN) = FLGJY*PHI*(FI*ASUM+DFI*BSUM)/WK(7)
  488. FN = FN - FLGJY
  489. 170 CONTINUE
  490. RETURN
  491. END