|
@@ -1,491 +0,0 @@
|
|
|
-*DECK ASYJY
|
|
|
- SUBROUTINE ASYJY (FUNJY, X, FNU, FLGJY, IN, Y, WK, IFLW)
|
|
|
-C***BEGIN PROLOGUE ASYJY
|
|
|
-C***SUBSIDIARY
|
|
|
-C***PURPOSE Subsidiary to BESJ and BESY
|
|
|
-C***LIBRARY SLATEC
|
|
|
-C***TYPE SINGLE PRECISION (ASYJY-S, DASYJY-D)
|
|
|
-C***AUTHOR Amos, D. E., (SNLA)
|
|
|
-C***DESCRIPTION
|
|
|
-C
|
|
|
-C ASYJY computes Bessel functions J and Y
|
|
|
-C for arguments X.GT.0.0 and orders FNU.GE.35.0
|
|
|
-C on FLGJY = 1 and FLGJY = -1 respectively
|
|
|
-C
|
|
|
-C INPUT
|
|
|
-C
|
|
|
-C FUNJY - external function JAIRY or YAIRY
|
|
|
-C X - argument, X.GT.0.0E0
|
|
|
-C FNU - order of the first Bessel function
|
|
|
-C FLGJY - selection flag
|
|
|
-C FLGJY = 1.0E0 gives the J function
|
|
|
-C FLGJY = -1.0E0 gives the Y function
|
|
|
-C IN - number of functions desired, IN = 1 or 2
|
|
|
-C
|
|
|
-C OUTPUT
|
|
|
-C
|
|
|
-C Y - a vector whose first in components contain the sequence
|
|
|
-C IFLW - a flag indicating underflow or overflow
|
|
|
-C return variables for BESJ only
|
|
|
-C WK(1) = 1 - (X/FNU)**2 = W**2
|
|
|
-C WK(2) = SQRT(ABS(WK(1)))
|
|
|
-C WK(3) = ABS(WK(2) - ATAN(WK(2))) or
|
|
|
-C ABS(LN((1 + WK(2))/(X/FNU)) - WK(2))
|
|
|
-C = ABS((2/3)*ZETA**(3/2))
|
|
|
-C WK(4) = FNU*WK(3)
|
|
|
-C WK(5) = (1.5*WK(3)*FNU)**(1/3) = SQRT(ZETA)*FNU**(1/3)
|
|
|
-C WK(6) = SIGN(1.,W**2)*WK(5)**2 = SIGN(1.,W**2)*ZETA*FNU**(2/3)
|
|
|
-C WK(7) = FNU**(1/3)
|
|
|
-C
|
|
|
-C Abstract
|
|
|
-C ASYJY implements the uniform asymptotic expansion of
|
|
|
-C the J and Y Bessel functions for FNU.GE.35 and real
|
|
|
-C X.GT.0.0E0. The forms are identical except for a change
|
|
|
-C in sign of some of the terms. This change in sign is
|
|
|
-C accomplished by means of the flag FLGJY = 1 or -1. On
|
|
|
-C FLGJY = 1 the AIRY functions AI(X) and DAI(X) are
|
|
|
-C supplied by the external function JAIRY, and on
|
|
|
-C FLGJY = -1 the AIRY functions BI(X) and DBI(X) are
|
|
|
-C supplied by the external function YAIRY.
|
|
|
-C
|
|
|
-C***SEE ALSO BESJ, BESY
|
|
|
-C***ROUTINES CALLED I1MACH, R1MACH
|
|
|
-C***REVISION HISTORY (YYMMDD)
|
|
|
-C 750101 DATE WRITTEN
|
|
|
-C 890531 Changed all specific intrinsics to generic. (WRB)
|
|
|
-C 891009 Removed unreferenced variable. (WRB)
|
|
|
-C 891214 Prologue converted to Version 4.0 format. (BAB)
|
|
|
-C 900328 Added TYPE section. (WRB)
|
|
|
-C 910408 Updated the AUTHOR section. (WRB)
|
|
|
-C***END PROLOGUE ASYJY
|
|
|
- INTEGER I, IFLW, IN, J, JN,JR,JU,K, KB,KLAST,KMAX,KP1, KS, KSP1,
|
|
|
- * KSTEMP, L, LR, LRP1, ISETA, ISETB
|
|
|
- INTEGER I1MACH
|
|
|
- REAL ABW2, AKM, ALFA, ALFA1, ALFA2, AP, AR, ASUM, AZ,
|
|
|
- * BETA, BETA1, BETA2, BETA3, BR, BSUM, C, CON1, CON2,
|
|
|
- * CON548,CR,CRZ32, DFI,ELIM, DR,FI, FLGJY, FN, FNU,
|
|
|
- * FN2, GAMA, PHI, RCZ, RDEN, RELB, RFN2, RTZ, RZDEN,
|
|
|
- * SA, SB, SUMA, SUMB, S1, TA, TAU, TB, TFN, TOL, TOLS, T2, UPOL,
|
|
|
- * WK, X, XX, Y, Z, Z32
|
|
|
- REAL R1MACH
|
|
|
- DIMENSION Y(*), WK(*), C(65)
|
|
|
- DIMENSION ALFA(26,4), BETA(26,5)
|
|
|
- DIMENSION ALFA1(26,2), ALFA2(26,2)
|
|
|
- DIMENSION BETA1(26,2), BETA2(26,2), BETA3(26,1)
|
|
|
- DIMENSION GAMA(26), KMAX(5), AR(8), BR(10), UPOL(10)
|
|
|
- DIMENSION CR(10), DR(10)
|
|
|
- EQUIVALENCE (ALFA(1,1),ALFA1(1,1))
|
|
|
- EQUIVALENCE (ALFA(1,3),ALFA2(1,1))
|
|
|
- EQUIVALENCE (BETA(1,1),BETA1(1,1))
|
|
|
- EQUIVALENCE (BETA(1,3),BETA2(1,1))
|
|
|
- EQUIVALENCE (BETA(1,5),BETA3(1,1))
|
|
|
- SAVE TOLS, CON1, CON2, CON548, AR, BR, C, ALFA1, ALFA2,
|
|
|
- 1 BETA1, BETA2, BETA3, GAMA
|
|
|
- DATA TOLS /-6.90775527898214E+00/
|
|
|
- DATA CON1,CON2,CON548/
|
|
|
- 1 6.66666666666667E-01, 3.33333333333333E-01, 1.04166666666667E-01/
|
|
|
- DATA AR(1), AR(2), AR(3), AR(4), AR(5), AR(6), AR(7),
|
|
|
- A AR(8) / 8.35503472222222E-02, 1.28226574556327E-01,
|
|
|
- 1 2.91849026464140E-01, 8.81627267443758E-01, 3.32140828186277E+00,
|
|
|
- 2 1.49957629868626E+01, 7.89230130115865E+01, 4.74451538868264E+02/
|
|
|
- DATA BR(1), BR(2), BR(3), BR(4), BR(5), BR(6), BR(7), BR(8),
|
|
|
- A BR(9), BR(10) /-1.45833333333333E-01,-9.87413194444444E-02,
|
|
|
- 1-1.43312053915895E-01,-3.17227202678414E-01,-9.42429147957120E-01,
|
|
|
- 2-3.51120304082635E+00,-1.57272636203680E+01,-8.22814390971859E+01,
|
|
|
- 3-4.92355370523671E+02,-3.31621856854797E+03/
|
|
|
- DATA C(1), C(2), C(3), C(4), C(5), C(6), C(7), C(8), C(9), C(10),
|
|
|
- 1 C(11), C(12), C(13), C(14), C(15), C(16), C(17), C(18),
|
|
|
- 2 C(19), C(20), C(21), C(22), C(23), C(24)/
|
|
|
- 3 -2.08333333333333E-01, 1.25000000000000E-01,
|
|
|
- 4 3.34201388888889E-01, -4.01041666666667E-01,
|
|
|
- 5 7.03125000000000E-02, -1.02581259645062E+00,
|
|
|
- 6 1.84646267361111E+00, -8.91210937500000E-01,
|
|
|
- 7 7.32421875000000E-02, 4.66958442342625E+00,
|
|
|
- 8 -1.12070026162230E+01, 8.78912353515625E+00,
|
|
|
- 9 -2.36408691406250E+00, 1.12152099609375E-01,
|
|
|
- A -2.82120725582002E+01, 8.46362176746007E+01,
|
|
|
- B -9.18182415432400E+01, 4.25349987453885E+01,
|
|
|
- C -7.36879435947963E+00, 2.27108001708984E-01,
|
|
|
- D 2.12570130039217E+02, -7.65252468141182E+02,
|
|
|
- E 1.05999045252800E+03, -6.99579627376133E+02/
|
|
|
- DATA C(25), C(26), C(27), C(28), C(29), C(30), C(31), C(32),
|
|
|
- 1 C(33), C(34), C(35), C(36), C(37), C(38), C(39), C(40),
|
|
|
- 2 C(41), C(42), C(43), C(44), C(45), C(46), C(47), C(48)/
|
|
|
- 3 2.18190511744212E+02, -2.64914304869516E+01,
|
|
|
- 4 5.72501420974731E-01, -1.91945766231841E+03,
|
|
|
- 5 8.06172218173731E+03, -1.35865500064341E+04,
|
|
|
- 6 1.16553933368645E+04, -5.30564697861340E+03,
|
|
|
- 7 1.20090291321635E+03, -1.08090919788395E+02,
|
|
|
- 8 1.72772750258446E+00, 2.02042913309661E+04,
|
|
|
- 9 -9.69805983886375E+04, 1.92547001232532E+05,
|
|
|
- A -2.03400177280416E+05, 1.22200464983017E+05,
|
|
|
- B -4.11926549688976E+04, 7.10951430248936E+03,
|
|
|
- C -4.93915304773088E+02, 6.07404200127348E+00,
|
|
|
- D -2.42919187900551E+05, 1.31176361466298E+06,
|
|
|
- E -2.99801591853811E+06, 3.76327129765640E+06/
|
|
|
- DATA C(49), C(50), C(51), C(52), C(53), C(54), C(55), C(56),
|
|
|
- 1 C(57), C(58), C(59), C(60), C(61), C(62), C(63), C(64),
|
|
|
- 2 C(65)/
|
|
|
- 3 -2.81356322658653E+06, 1.26836527332162E+06,
|
|
|
- 4 -3.31645172484564E+05, 4.52187689813627E+04,
|
|
|
- 5 -2.49983048181121E+03, 2.43805296995561E+01,
|
|
|
- 6 3.28446985307204E+06, -1.97068191184322E+07,
|
|
|
- 7 5.09526024926646E+07, -7.41051482115327E+07,
|
|
|
- 8 6.63445122747290E+07, -3.75671766607634E+07,
|
|
|
- 9 1.32887671664218E+07, -2.78561812808645E+06,
|
|
|
- A 3.08186404612662E+05, -1.38860897537170E+04,
|
|
|
- B 1.10017140269247E+02/
|
|
|
- DATA ALFA1(1,1), ALFA1(2,1), ALFA1(3,1), ALFA1(4,1), ALFA1(5,1),
|
|
|
- 1 ALFA1(6,1), ALFA1(7,1), ALFA1(8,1), ALFA1(9,1), ALFA1(10,1),
|
|
|
- 2 ALFA1(11,1),ALFA1(12,1),ALFA1(13,1),ALFA1(14,1),ALFA1(15,1),
|
|
|
- 3 ALFA1(16,1),ALFA1(17,1),ALFA1(18,1),ALFA1(19,1),ALFA1(20,1),
|
|
|
- 4 ALFA1(21,1),ALFA1(22,1),ALFA1(23,1),ALFA1(24,1),ALFA1(25,1),
|
|
|
- 5 ALFA1(26,1) /-4.44444444444444E-03,-9.22077922077922E-04,
|
|
|
- 6-8.84892884892885E-05, 1.65927687832450E-04, 2.46691372741793E-04,
|
|
|
- 7 2.65995589346255E-04, 2.61824297061501E-04, 2.48730437344656E-04,
|
|
|
- 8 2.32721040083232E-04, 2.16362485712365E-04, 2.00738858762752E-04,
|
|
|
- 9 1.86267636637545E-04, 1.73060775917876E-04, 1.61091705929016E-04,
|
|
|
- 1 1.50274774160908E-04, 1.40503497391270E-04, 1.31668816545923E-04,
|
|
|
- 2 1.23667445598253E-04, 1.16405271474738E-04, 1.09798298372713E-04,
|
|
|
- 3 1.03772410422993E-04, 9.82626078369363E-05, 9.32120517249503E-05,
|
|
|
- 4 8.85710852478712E-05, 8.42963105715700E-05, 8.03497548407791E-05/
|
|
|
- DATA ALFA1(1,2), ALFA1(2,2), ALFA1(3,2), ALFA1(4,2), ALFA1(5,2),
|
|
|
- 1 ALFA1(6,2), ALFA1(7,2), ALFA1(8,2), ALFA1(9,2), ALFA1(10,2),
|
|
|
- 2 ALFA1(11,2),ALFA1(12,2),ALFA1(13,2),ALFA1(14,2),ALFA1(15,2),
|
|
|
- 3 ALFA1(16,2),ALFA1(17,2),ALFA1(18,2),ALFA1(19,2),ALFA1(20,2),
|
|
|
- 4 ALFA1(21,2),ALFA1(22,2),ALFA1(23,2),ALFA1(24,2),ALFA1(25,2),
|
|
|
- 5 ALFA1(26,2) / 6.93735541354589E-04, 2.32241745182922E-04,
|
|
|
- 6-1.41986273556691E-05,-1.16444931672049E-04,-1.50803558053049E-04,
|
|
|
- 7-1.55121924918096E-04,-1.46809756646466E-04,-1.33815503867491E-04,
|
|
|
- 8-1.19744975684254E-04,-1.06184319207974E-04,-9.37699549891194E-05,
|
|
|
- 9-8.26923045588193E-05,-7.29374348155221E-05,-6.44042357721016E-05,
|
|
|
- 1-5.69611566009369E-05,-5.04731044303562E-05,-4.48134868008883E-05,
|
|
|
- 2-3.98688727717599E-05,-3.55400532972042E-05,-3.17414256609022E-05,
|
|
|
- 3-2.83996793904175E-05,-2.54522720634871E-05,-2.28459297164725E-05,
|
|
|
- 4-2.05352753106481E-05,-1.84816217627666E-05,-1.66519330021394E-05/
|
|
|
- DATA ALFA2(1,1), ALFA2(2,1), ALFA2(3,1), ALFA2(4,1), ALFA2(5,1),
|
|
|
- 1 ALFA2(6,1), ALFA2(7,1), ALFA2(8,1), ALFA2(9,1), ALFA2(10,1),
|
|
|
- 2 ALFA2(11,1),ALFA2(12,1),ALFA2(13,1),ALFA2(14,1),ALFA2(15,1),
|
|
|
- 3 ALFA2(16,1),ALFA2(17,1),ALFA2(18,1),ALFA2(19,1),ALFA2(20,1),
|
|
|
- 4 ALFA2(21,1),ALFA2(22,1),ALFA2(23,1),ALFA2(24,1),ALFA2(25,1),
|
|
|
- 5 ALFA2(26,1) /-3.54211971457744E-04,-1.56161263945159E-04,
|
|
|
- 6 3.04465503594936E-05, 1.30198655773243E-04, 1.67471106699712E-04,
|
|
|
- 7 1.70222587683593E-04, 1.56501427608595E-04, 1.36339170977445E-04,
|
|
|
- 8 1.14886692029825E-04, 9.45869093034688E-05, 7.64498419250898E-05,
|
|
|
- 9 6.07570334965197E-05, 4.74394299290509E-05, 3.62757512005344E-05,
|
|
|
- 1 2.69939714979225E-05, 1.93210938247939E-05, 1.30056674793963E-05,
|
|
|
- 2 7.82620866744497E-06, 3.59257485819352E-06, 1.44040049814252E-07,
|
|
|
- 3-2.65396769697939E-06,-4.91346867098486E-06,-6.72739296091248E-06,
|
|
|
- 4-8.17269379678658E-06,-9.31304715093561E-06,-1.02011418798016E-05/
|
|
|
- DATA ALFA2(1,2), ALFA2(2,2), ALFA2(3,2), ALFA2(4,2), ALFA2(5,2),
|
|
|
- 1 ALFA2(6,2), ALFA2(7,2), ALFA2(8,2), ALFA2(9,2), ALFA2(10,2),
|
|
|
- 2 ALFA2(11,2),ALFA2(12,2),ALFA2(13,2),ALFA2(14,2),ALFA2(15,2),
|
|
|
- 3 ALFA2(16,2),ALFA2(17,2),ALFA2(18,2),ALFA2(19,2),ALFA2(20,2),
|
|
|
- 4 ALFA2(21,2),ALFA2(22,2),ALFA2(23,2),ALFA2(24,2),ALFA2(25,2),
|
|
|
- 5 ALFA2(26,2) / 3.78194199201773E-04, 2.02471952761816E-04,
|
|
|
- 6-6.37938506318862E-05,-2.38598230603006E-04,-3.10916256027362E-04,
|
|
|
- 7-3.13680115247576E-04,-2.78950273791323E-04,-2.28564082619141E-04,
|
|
|
- 8-1.75245280340847E-04,-1.25544063060690E-04,-8.22982872820208E-05,
|
|
|
- 9-4.62860730588116E-05,-1.72334302366962E-05, 5.60690482304602E-06,
|
|
|
- 1 2.31395443148287E-05, 3.62642745856794E-05, 4.58006124490189E-05,
|
|
|
- 2 5.24595294959114E-05, 5.68396208545815E-05, 5.94349820393104E-05,
|
|
|
- 3 6.06478527578422E-05, 6.08023907788436E-05, 6.01577894539460E-05,
|
|
|
- 4 5.89199657344698E-05, 5.72515823777593E-05, 5.52804375585853E-05/
|
|
|
- DATA BETA1(1,1), BETA1(2,1), BETA1(3,1), BETA1(4,1), BETA1(5,1),
|
|
|
- 1 BETA1(6,1), BETA1(7,1), BETA1(8,1), BETA1(9,1), BETA1(10,1),
|
|
|
- 2 BETA1(11,1),BETA1(12,1),BETA1(13,1),BETA1(14,1),BETA1(15,1),
|
|
|
- 3 BETA1(16,1),BETA1(17,1),BETA1(18,1),BETA1(19,1),BETA1(20,1),
|
|
|
- 4 BETA1(21,1),BETA1(22,1),BETA1(23,1),BETA1(24,1),BETA1(25,1),
|
|
|
- 5 BETA1(26,1) / 1.79988721413553E-02, 5.59964911064388E-03,
|
|
|
- 6 2.88501402231133E-03, 1.80096606761054E-03, 1.24753110589199E-03,
|
|
|
- 7 9.22878876572938E-04, 7.14430421727287E-04, 5.71787281789705E-04,
|
|
|
- 8 4.69431007606482E-04, 3.93232835462917E-04, 3.34818889318298E-04,
|
|
|
- 9 2.88952148495752E-04, 2.52211615549573E-04, 2.22280580798883E-04,
|
|
|
- 1 1.97541838033063E-04, 1.76836855019718E-04, 1.59316899661821E-04,
|
|
|
- 2 1.44347930197334E-04, 1.31448068119965E-04, 1.20245444949303E-04,
|
|
|
- 3 1.10449144504599E-04, 1.01828770740567E-04, 9.41998224204238E-05,
|
|
|
- 4 8.74130545753834E-05, 8.13466262162801E-05, 7.59002269646219E-05/
|
|
|
- DATA BETA1(1,2), BETA1(2,2), BETA1(3,2), BETA1(4,2), BETA1(5,2),
|
|
|
- 1 BETA1(6,2), BETA1(7,2), BETA1(8,2), BETA1(9,2), BETA1(10,2),
|
|
|
- 2 BETA1(11,2),BETA1(12,2),BETA1(13,2),BETA1(14,2),BETA1(15,2),
|
|
|
- 3 BETA1(16,2),BETA1(17,2),BETA1(18,2),BETA1(19,2),BETA1(20,2),
|
|
|
- 4 BETA1(21,2),BETA1(22,2),BETA1(23,2),BETA1(24,2),BETA1(25,2),
|
|
|
- 5 BETA1(26,2) /-1.49282953213429E-03,-8.78204709546389E-04,
|
|
|
- 6-5.02916549572035E-04,-2.94822138512746E-04,-1.75463996970783E-04,
|
|
|
- 7-1.04008550460816E-04,-5.96141953046458E-05,-3.12038929076098E-05,
|
|
|
- 8-1.26089735980230E-05,-2.42892608575730E-07, 8.05996165414274E-06,
|
|
|
- 9 1.36507009262147E-05, 1.73964125472926E-05, 1.98672978842134E-05,
|
|
|
- 1 2.14463263790823E-05, 2.23954659232457E-05, 2.28967783814713E-05,
|
|
|
- 2 2.30785389811178E-05, 2.30321976080909E-05, 2.28236073720349E-05,
|
|
|
- 3 2.25005881105292E-05, 2.20981015361991E-05, 2.16418427448104E-05,
|
|
|
- 4 2.11507649256221E-05, 2.06388749782171E-05, 2.01165241997082E-05/
|
|
|
- DATA BETA2(1,1), BETA2(2,1), BETA2(3,1), BETA2(4,1), BETA2(5,1),
|
|
|
- 1 BETA2(6,1), BETA2(7,1), BETA2(8,1), BETA2(9,1), BETA2(10,1),
|
|
|
- 2 BETA2(11,1),BETA2(12,1),BETA2(13,1),BETA2(14,1),BETA2(15,1),
|
|
|
- 3 BETA2(16,1),BETA2(17,1),BETA2(18,1),BETA2(19,1),BETA2(20,1),
|
|
|
- 4 BETA2(21,1),BETA2(22,1),BETA2(23,1),BETA2(24,1),BETA2(25,1),
|
|
|
- 5 BETA2(26,1) / 5.52213076721293E-04, 4.47932581552385E-04,
|
|
|
- 6 2.79520653992021E-04, 1.52468156198447E-04, 6.93271105657044E-05,
|
|
|
- 7 1.76258683069991E-05,-1.35744996343269E-05,-3.17972413350427E-05,
|
|
|
- 8-4.18861861696693E-05,-4.69004889379141E-05,-4.87665447413787E-05,
|
|
|
- 9-4.87010031186735E-05,-4.74755620890087E-05,-4.55813058138628E-05,
|
|
|
- 1-4.33309644511266E-05,-4.09230193157750E-05,-3.84822638603221E-05,
|
|
|
- 2-3.60857167535411E-05,-3.37793306123367E-05,-3.15888560772110E-05,
|
|
|
- 3-2.95269561750807E-05,-2.75978914828336E-05,-2.58006174666884E-05,
|
|
|
- 4-2.41308356761280E-05,-2.25823509518346E-05,-2.11479656768913E-05/
|
|
|
- DATA BETA2(1,2), BETA2(2,2), BETA2(3,2), BETA2(4,2), BETA2(5,2),
|
|
|
- 1 BETA2(6,2), BETA2(7,2), BETA2(8,2), BETA2(9,2), BETA2(10,2),
|
|
|
- 2 BETA2(11,2),BETA2(12,2),BETA2(13,2),BETA2(14,2),BETA2(15,2),
|
|
|
- 3 BETA2(16,2),BETA2(17,2),BETA2(18,2),BETA2(19,2),BETA2(20,2),
|
|
|
- 4 BETA2(21,2),BETA2(22,2),BETA2(23,2),BETA2(24,2),BETA2(25,2),
|
|
|
- 5 BETA2(26,2) /-4.74617796559960E-04,-4.77864567147321E-04,
|
|
|
- 6-3.20390228067038E-04,-1.61105016119962E-04,-4.25778101285435E-05,
|
|
|
- 7 3.44571294294968E-05, 7.97092684075675E-05, 1.03138236708272E-04,
|
|
|
- 8 1.12466775262204E-04, 1.13103642108481E-04, 1.08651634848774E-04,
|
|
|
- 9 1.01437951597662E-04, 9.29298396593364E-05, 8.40293133016090E-05,
|
|
|
- 1 7.52727991349134E-05, 6.69632521975731E-05, 5.92564547323195E-05,
|
|
|
- 2 5.22169308826976E-05, 4.58539485165361E-05, 4.01445513891487E-05,
|
|
|
- 3 3.50481730031328E-05, 3.05157995034347E-05, 2.64956119950516E-05,
|
|
|
- 4 2.29363633690998E-05, 1.97893056664022E-05, 1.70091984636413E-05/
|
|
|
- DATA BETA3(1,1), BETA3(2,1), BETA3(3,1), BETA3(4,1), BETA3(5,1),
|
|
|
- 1 BETA3(6,1), BETA3(7,1), BETA3(8,1), BETA3(9,1), BETA3(10,1),
|
|
|
- 2 BETA3(11,1),BETA3(12,1),BETA3(13,1),BETA3(14,1),BETA3(15,1),
|
|
|
- 3 BETA3(16,1),BETA3(17,1),BETA3(18,1),BETA3(19,1),BETA3(20,1),
|
|
|
- 4 BETA3(21,1),BETA3(22,1),BETA3(23,1),BETA3(24,1),BETA3(25,1),
|
|
|
- 5 BETA3(26,1) / 7.36465810572578E-04, 8.72790805146194E-04,
|
|
|
- 6 6.22614862573135E-04, 2.85998154194304E-04, 3.84737672879366E-06,
|
|
|
- 7-1.87906003636972E-04,-2.97603646594555E-04,-3.45998126832656E-04,
|
|
|
- 8-3.53382470916038E-04,-3.35715635775049E-04,-3.04321124789040E-04,
|
|
|
- 9-2.66722723047613E-04,-2.27654214122820E-04,-1.89922611854562E-04,
|
|
|
- 1-1.55058918599094E-04,-1.23778240761874E-04,-9.62926147717644E-05,
|
|
|
- 2-7.25178327714425E-05,-5.22070028895634E-05,-3.50347750511901E-05,
|
|
|
- 3-2.06489761035552E-05,-8.70106096849767E-06, 1.13698686675100E-06,
|
|
|
- 4 9.16426474122779E-06, 1.56477785428873E-05, 2.08223629482467E-05/
|
|
|
- DATA GAMA(1), GAMA(2), GAMA(3), GAMA(4), GAMA(5),
|
|
|
- 1 GAMA(6), GAMA(7), GAMA(8), GAMA(9), GAMA(10),
|
|
|
- 2 GAMA(11), GAMA(12), GAMA(13), GAMA(14), GAMA(15),
|
|
|
- 3 GAMA(16), GAMA(17), GAMA(18), GAMA(19), GAMA(20),
|
|
|
- 4 GAMA(21), GAMA(22), GAMA(23), GAMA(24), GAMA(25),
|
|
|
- 5 GAMA(26) / 6.29960524947437E-01, 2.51984209978975E-01,
|
|
|
- 6 1.54790300415656E-01, 1.10713062416159E-01, 8.57309395527395E-02,
|
|
|
- 7 6.97161316958684E-02, 5.86085671893714E-02, 5.04698873536311E-02,
|
|
|
- 8 4.42600580689155E-02, 3.93720661543510E-02, 3.54283195924455E-02,
|
|
|
- 9 3.21818857502098E-02, 2.94646240791158E-02, 2.71581677112934E-02,
|
|
|
- 1 2.51768272973862E-02, 2.34570755306079E-02, 2.19508390134907E-02,
|
|
|
- 2 2.06210828235646E-02, 1.94388240897881E-02, 1.83810633800683E-02,
|
|
|
- 3 1.74293213231963E-02, 1.65685837786612E-02, 1.57865285987918E-02,
|
|
|
- 4 1.50729501494096E-02, 1.44193250839955E-02, 1.38184805735342E-02/
|
|
|
-C***FIRST EXECUTABLE STATEMENT ASYJY
|
|
|
- TA = R1MACH(3)
|
|
|
- TOL = MAX(TA,1.0E-15)
|
|
|
- TB = R1MACH(5)
|
|
|
- JU = I1MACH(12)
|
|
|
- IF(FLGJY.EQ.1.0E0) GO TO 6
|
|
|
- JR = I1MACH(11)
|
|
|
- ELIM = -2.303E0*TB*(JU+JR)
|
|
|
- GO TO 7
|
|
|
- 6 CONTINUE
|
|
|
- ELIM = -2.303E0*(TB*JU+3.0E0)
|
|
|
- 7 CONTINUE
|
|
|
- FN = FNU
|
|
|
- IFLW = 0
|
|
|
- DO 170 JN=1,IN
|
|
|
- XX = X/FN
|
|
|
- WK(1) = 1.0E0 - XX*XX
|
|
|
- ABW2 = ABS(WK(1))
|
|
|
- WK(2) = SQRT(ABW2)
|
|
|
- WK(7) = FN**CON2
|
|
|
- IF (ABW2.GT.0.27750E0) GO TO 80
|
|
|
-C
|
|
|
-C ASYMPTOTIC EXPANSION
|
|
|
-C CASES NEAR X=FN, ABS(1.-(X/FN)**2).LE.0.2775
|
|
|
-C COEFFICIENTS OF ASYMPTOTIC EXPANSION BY SERIES
|
|
|
-C
|
|
|
-C ZETA AND TRUNCATION FOR A(ZETA) AND B(ZETA) SERIES
|
|
|
-C
|
|
|
-C KMAX IS TRUNCATION INDEX FOR A(ZETA) AND B(ZETA) SERIES=MAX(2,SA)
|
|
|
-C
|
|
|
- SA = 0.0E0
|
|
|
- IF (ABW2.EQ.0.0E0) GO TO 10
|
|
|
- SA = TOLS/LOG(ABW2)
|
|
|
- 10 SB = SA
|
|
|
- DO 20 I=1,5
|
|
|
- AKM = MAX(SA,2.0E0)
|
|
|
- KMAX(I) = INT(AKM)
|
|
|
- SA = SA + SB
|
|
|
- 20 CONTINUE
|
|
|
- KB = KMAX(5)
|
|
|
- KLAST = KB - 1
|
|
|
- SA = GAMA(KB)
|
|
|
- DO 30 K=1,KLAST
|
|
|
- KB = KB - 1
|
|
|
- SA = SA*WK(1) + GAMA(KB)
|
|
|
- 30 CONTINUE
|
|
|
- Z = WK(1)*SA
|
|
|
- AZ = ABS(Z)
|
|
|
- RTZ = SQRT(AZ)
|
|
|
- WK(3) = CON1*AZ*RTZ
|
|
|
- WK(4) = WK(3)*FN
|
|
|
- WK(5) = RTZ*WK(7)
|
|
|
- WK(6) = -WK(5)*WK(5)
|
|
|
- IF(Z.LE.0.0E0) GO TO 35
|
|
|
- IF(WK(4).GT.ELIM) GO TO 75
|
|
|
- WK(6) = -WK(6)
|
|
|
- 35 CONTINUE
|
|
|
- PHI = SQRT(SQRT(SA+SA+SA+SA))
|
|
|
-C
|
|
|
-C B(ZETA) FOR S=0
|
|
|
-C
|
|
|
- KB = KMAX(5)
|
|
|
- KLAST = KB - 1
|
|
|
- SB = BETA(KB,1)
|
|
|
- DO 40 K=1,KLAST
|
|
|
- KB = KB - 1
|
|
|
- SB = SB*WK(1) + BETA(KB,1)
|
|
|
- 40 CONTINUE
|
|
|
- KSP1 = 1
|
|
|
- FN2 = FN*FN
|
|
|
- RFN2 = 1.0E0/FN2
|
|
|
- RDEN = 1.0E0
|
|
|
- ASUM = 1.0E0
|
|
|
- RELB = TOL*ABS(SB)
|
|
|
- BSUM = SB
|
|
|
- DO 60 KS=1,4
|
|
|
- KSP1 = KSP1 + 1
|
|
|
- RDEN = RDEN*RFN2
|
|
|
-C
|
|
|
-C A(ZETA) AND B(ZETA) FOR S=1,2,3,4
|
|
|
-C
|
|
|
- KSTEMP = 5 - KS
|
|
|
- KB = KMAX(KSTEMP)
|
|
|
- KLAST = KB - 1
|
|
|
- SA = ALFA(KB,KS)
|
|
|
- SB = BETA(KB,KSP1)
|
|
|
- DO 50 K=1,KLAST
|
|
|
- KB = KB - 1
|
|
|
- SA = SA*WK(1) + ALFA(KB,KS)
|
|
|
- SB = SB*WK(1) + BETA(KB,KSP1)
|
|
|
- 50 CONTINUE
|
|
|
- TA = SA*RDEN
|
|
|
- TB = SB*RDEN
|
|
|
- ASUM = ASUM + TA
|
|
|
- BSUM = BSUM + TB
|
|
|
- IF (ABS(TA).LE.TOL .AND. ABS(TB).LE.RELB) GO TO 70
|
|
|
- 60 CONTINUE
|
|
|
- 70 CONTINUE
|
|
|
- BSUM = BSUM/(FN*WK(7))
|
|
|
- GO TO 160
|
|
|
-C
|
|
|
- 75 CONTINUE
|
|
|
- IFLW = 1
|
|
|
- RETURN
|
|
|
-C
|
|
|
- 80 CONTINUE
|
|
|
- UPOL(1) = 1.0E0
|
|
|
- TAU = 1.0E0/WK(2)
|
|
|
- T2 = 1.0E0/WK(1)
|
|
|
- IF (WK(1).GE.0.0E0) GO TO 90
|
|
|
-C
|
|
|
-C CASES FOR (X/FN).GT.SQRT(1.2775)
|
|
|
-C
|
|
|
- WK(3) = ABS(WK(2)-ATAN(WK(2)))
|
|
|
- WK(4) = WK(3)*FN
|
|
|
- RCZ = -CON1/WK(4)
|
|
|
- Z32 = 1.5E0*WK(3)
|
|
|
- RTZ = Z32**CON2
|
|
|
- WK(5) = RTZ*WK(7)
|
|
|
- WK(6) = -WK(5)*WK(5)
|
|
|
- GO TO 100
|
|
|
- 90 CONTINUE
|
|
|
-C
|
|
|
-C CASES FOR (X/FN).LT.SQRT(0.7225)
|
|
|
-C
|
|
|
- WK(3) = ABS(LOG((1.0E0+WK(2))/XX)-WK(2))
|
|
|
- WK(4) = WK(3)*FN
|
|
|
- RCZ = CON1/WK(4)
|
|
|
- IF(WK(4).GT.ELIM) GO TO 75
|
|
|
- Z32 = 1.5E0*WK(3)
|
|
|
- RTZ = Z32**CON2
|
|
|
- WK(7) = FN**CON2
|
|
|
- WK(5) = RTZ*WK(7)
|
|
|
- WK(6) = WK(5)*WK(5)
|
|
|
- 100 CONTINUE
|
|
|
- PHI = SQRT((RTZ+RTZ)*TAU)
|
|
|
- TB = 1.0E0
|
|
|
- ASUM = 1.0E0
|
|
|
- TFN = TAU/FN
|
|
|
- RDEN=1.0E0/FN
|
|
|
- RFN2=RDEN*RDEN
|
|
|
- RDEN=1.0E0
|
|
|
- UPOL(2) = (C(1)*T2+C(2))*TFN
|
|
|
- CRZ32 = CON548*RCZ
|
|
|
- BSUM = UPOL(2) + CRZ32
|
|
|
- RELB = TOL*ABS(BSUM)
|
|
|
- AP = TFN
|
|
|
- KS = 0
|
|
|
- KP1 = 2
|
|
|
- RZDEN = RCZ
|
|
|
- L = 2
|
|
|
- ISETA=0
|
|
|
- ISETB=0
|
|
|
- DO 140 LR=2,8,2
|
|
|
-C
|
|
|
-C COMPUTE TWO U POLYNOMIALS FOR NEXT A(ZETA) AND B(ZETA)
|
|
|
-C
|
|
|
- LRP1 = LR + 1
|
|
|
- DO 120 K=LR,LRP1
|
|
|
- KS = KS + 1
|
|
|
- KP1 = KP1 + 1
|
|
|
- L = L + 1
|
|
|
- S1 = C(L)
|
|
|
- DO 110 J=2,KP1
|
|
|
- L = L + 1
|
|
|
- S1 = S1*T2 + C(L)
|
|
|
- 110 CONTINUE
|
|
|
- AP = AP*TFN
|
|
|
- UPOL(KP1) = AP*S1
|
|
|
- CR(KS) = BR(KS)*RZDEN
|
|
|
- RZDEN = RZDEN*RCZ
|
|
|
- DR(KS) = AR(KS)*RZDEN
|
|
|
- 120 CONTINUE
|
|
|
- SUMA = UPOL(LRP1)
|
|
|
- SUMB = UPOL(LR+2) + UPOL(LRP1)*CRZ32
|
|
|
- JU = LRP1
|
|
|
- DO 130 JR=1,LR
|
|
|
- JU = JU - 1
|
|
|
- SUMA = SUMA + CR(JR)*UPOL(JU)
|
|
|
- SUMB = SUMB + DR(JR)*UPOL(JU)
|
|
|
- 130 CONTINUE
|
|
|
- RDEN=RDEN*RFN2
|
|
|
- TB = -TB
|
|
|
- IF (WK(1).GT.0.0E0) TB = ABS(TB)
|
|
|
- IF (RDEN.LT.TOL) GO TO 131
|
|
|
- ASUM = ASUM + SUMA*TB
|
|
|
- BSUM = BSUM + SUMB*TB
|
|
|
- GO TO 140
|
|
|
- 131 IF(ISETA.EQ.1) GO TO 132
|
|
|
- IF(ABS(SUMA).LT.TOL) ISETA=1
|
|
|
- ASUM=ASUM+SUMA*TB
|
|
|
- 132 IF(ISETB.EQ.1) GO TO 133
|
|
|
- IF(ABS(SUMB).LT.RELB) ISETB=1
|
|
|
- BSUM=BSUM+SUMB*TB
|
|
|
- 133 IF(ISETA.EQ.1 .AND. ISETB.EQ.1) GO TO 150
|
|
|
- 140 CONTINUE
|
|
|
- 150 TB = WK(5)
|
|
|
- IF (WK(1).GT.0.0E0) TB = -TB
|
|
|
- BSUM = BSUM/TB
|
|
|
-C
|
|
|
- 160 CONTINUE
|
|
|
- CALL FUNJY(WK(6), WK(5), WK(4), FI, DFI)
|
|
|
- TA=1.0E0/TOL
|
|
|
- TB=R1MACH(1)*TA*1.0E+3
|
|
|
- IF(ABS(FI).GT.TB) GO TO 165
|
|
|
- FI=FI*TA
|
|
|
- DFI=DFI*TA
|
|
|
- PHI=PHI*TOL
|
|
|
- 165 CONTINUE
|
|
|
- Y(JN) = FLGJY*PHI*(FI*ASUM+DFI*BSUM)/WK(7)
|
|
|
- FN = FN - FLGJY
|
|
|
- 170 CONTINUE
|
|
|
- RETURN
|
|
|
- END
|