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