DOUBLE PRECISION FUNCTION DGAMLN (Z, IERR) C***BEGIN PROLOGUE DGAMLN C***SUBSIDIARY C***PURPOSE Compute the logarithm of the Gamma function C***LIBRARY SLATEC C***CATEGORY C7A C***TYPE DOUBLE PRECISION (GAMLN-S, DGAMLN-D) C***KEYWORDS LOGARITHM OF GAMMA FUNCTION C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C **** A DOUBLE PRECISION ROUTINE **** C DGAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR C Z.GT.0. THE ASYMPTOTIC EXPANSION IS USED TO GENERATE VALUES C GREATER THAN ZMIN WHICH ARE ADJUSTED BY THE RECURSION C G(Z+1)=Z*G(Z) FOR Z.LE.ZMIN. THE FUNCTION WAS MADE AS C PORTABLE AS POSSIBLE BY COMPUTIMG ZMIN FROM THE NUMBER OF BASE C 10 DIGITS IN A WORD, RLN=MAX(-ALOG10(R1MACH(4)),0.5E-18) C LIMITED TO 18 DIGITS OF (RELATIVE) ACCURACY. C C SINCE INTEGER ARGUMENTS ARE COMMON, A TABLE LOOK UP ON 100 C VALUES IS USED FOR SPEED OF EXECUTION. C C DESCRIPTION OF ARGUMENTS C C INPUT Z IS D0UBLE PRECISION C Z - ARGUMENT, Z.GT.0.0D0 C C OUTPUT DGAMLN IS DOUBLE PRECISION C DGAMLN - NATURAL LOG OF THE GAMMA FUNCTION AT Z.NE.0.0D0 C IERR - ERROR FLAG C IERR=0, NORMAL RETURN, COMPUTATION COMPLETED C IERR=1, Z.LE.0.0D0, NO COMPUTATION C C C***REFERENCES COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT C BY D. E. AMOS, SAND83-0083, MAY, 1983. C***ROUTINES CALLED D1MACH, I1MACH C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 830501 REVISION DATE from Version 3.2 C 910415 Prologue converted to Version 4.0 format. (BAB) C 920128 Category corrected. (WRB) C 921215 DGAMLN defined for Z negative. (WRB) C***END PROLOGUE DGAMLN DOUBLE PRECISION CF, CON, FLN, FZ, GLN, RLN, S, TLG, TRM, TST, * T1, WDTOL, Z, ZDMY, ZINC, ZM, ZMIN, ZP, ZSQ, D1MACH INTEGER I, IERR, I1M, K, MZ, NZ, I1MACH DIMENSION CF(22), GLN(100) C LNGAMMA(N), N=1,100 DATA GLN(1), GLN(2), GLN(3), GLN(4), GLN(5), GLN(6), GLN(7), 1 GLN(8), GLN(9), GLN(10), GLN(11), GLN(12), GLN(13), GLN(14), 2 GLN(15), GLN(16), GLN(17), GLN(18), GLN(19), GLN(20), 3 GLN(21), GLN(22)/ 4 0.00000000000000000D+00, 0.00000000000000000D+00, 5 6.93147180559945309D-01, 1.79175946922805500D+00, 6 3.17805383034794562D+00, 4.78749174278204599D+00, 7 6.57925121201010100D+00, 8.52516136106541430D+00, 8 1.06046029027452502D+01, 1.28018274800814696D+01, 9 1.51044125730755153D+01, 1.75023078458738858D+01, A 1.99872144956618861D+01, 2.25521638531234229D+01, B 2.51912211827386815D+01, 2.78992713838408916D+01, C 3.06718601060806728D+01, 3.35050734501368889D+01, D 3.63954452080330536D+01, 3.93398841871994940D+01, E 4.23356164607534850D+01, 4.53801388984769080D+01/ DATA GLN(23), GLN(24), GLN(25), GLN(26), GLN(27), GLN(28), 1 GLN(29), GLN(30), GLN(31), GLN(32), GLN(33), GLN(34), 2 GLN(35), GLN(36), GLN(37), GLN(38), GLN(39), GLN(40), 3 GLN(41), GLN(42), GLN(43), GLN(44)/ 4 4.84711813518352239D+01, 5.16066755677643736D+01, 5 5.47847293981123192D+01, 5.80036052229805199D+01, 6 6.12617017610020020D+01, 6.45575386270063311D+01, 7 6.78897431371815350D+01, 7.12570389671680090D+01, 8 7.46582363488301644D+01, 7.80922235533153106D+01, 9 8.15579594561150372D+01, 8.50544670175815174D+01, A 8.85808275421976788D+01, 9.21361756036870925D+01, B 9.57196945421432025D+01, 9.93306124547874269D+01, C 1.02968198614513813D+02, 1.06631760260643459D+02, D 1.10320639714757395D+02, 1.14034211781461703D+02, E 1.17771881399745072D+02, 1.21533081515438634D+02/ DATA GLN(45), GLN(46), GLN(47), GLN(48), GLN(49), GLN(50), 1 GLN(51), GLN(52), GLN(53), GLN(54), GLN(55), GLN(56), 2 GLN(57), GLN(58), GLN(59), GLN(60), GLN(61), GLN(62), 3 GLN(63), GLN(64), GLN(65), GLN(66)/ 4 1.25317271149356895D+02, 1.29123933639127215D+02, 5 1.32952575035616310D+02, 1.36802722637326368D+02, 6 1.40673923648234259D+02, 1.44565743946344886D+02, 7 1.48477766951773032D+02, 1.52409592584497358D+02, 8 1.56360836303078785D+02, 1.60331128216630907D+02, 9 1.64320112263195181D+02, 1.68327445448427652D+02, A 1.72352797139162802D+02, 1.76395848406997352D+02, B 1.80456291417543771D+02, 1.84533828861449491D+02, C 1.88628173423671591D+02, 1.92739047287844902D+02, D 1.96866181672889994D+02, 2.01009316399281527D+02, E 2.05168199482641199D+02, 2.09342586752536836D+02/ DATA GLN(67), GLN(68), GLN(69), GLN(70), GLN(71), GLN(72), 1 GLN(73), GLN(74), GLN(75), GLN(76), GLN(77), GLN(78), 2 GLN(79), GLN(80), GLN(81), GLN(82), GLN(83), GLN(84), 3 GLN(85), GLN(86), GLN(87), GLN(88)/ 4 2.13532241494563261D+02, 2.17736934113954227D+02, 5 2.21956441819130334D+02, 2.26190548323727593D+02, 6 2.30439043565776952D+02, 2.34701723442818268D+02, 7 2.38978389561834323D+02, 2.43268849002982714D+02, 8 2.47572914096186884D+02, 2.51890402209723194D+02, 9 2.56221135550009525D+02, 2.60564940971863209D+02, A 2.64921649798552801D+02, 2.69291097651019823D+02, B 2.73673124285693704D+02, 2.78067573440366143D+02, C 2.82474292687630396D+02, 2.86893133295426994D+02, D 2.91323950094270308D+02, 2.95766601350760624D+02, E 3.00220948647014132D+02, 3.04686856765668715D+02/ DATA GLN(89), GLN(90), GLN(91), GLN(92), GLN(93), GLN(94), 1 GLN(95), GLN(96), GLN(97), GLN(98), GLN(99), GLN(100)/ 2 3.09164193580146922D+02, 3.13652829949879062D+02, 3 3.18152639620209327D+02, 3.22663499126726177D+02, 4 3.27185287703775217D+02, 3.31717887196928473D+02, 5 3.36261181979198477D+02, 3.40815058870799018D+02, 6 3.45379407062266854D+02, 3.49954118040770237D+02, 7 3.54539085519440809D+02, 3.59134205369575399D+02/ C COEFFICIENTS OF ASYMPTOTIC EXPANSION DATA CF(1), CF(2), CF(3), CF(4), CF(5), CF(6), CF(7), CF(8), 1 CF(9), CF(10), CF(11), CF(12), CF(13), CF(14), CF(15), 2 CF(16), CF(17), CF(18), CF(19), CF(20), CF(21), CF(22)/ 3 8.33333333333333333D-02, -2.77777777777777778D-03, 4 7.93650793650793651D-04, -5.95238095238095238D-04, 5 8.41750841750841751D-04, -1.91752691752691753D-03, 6 6.41025641025641026D-03, -2.95506535947712418D-02, 7 1.79644372368830573D-01, -1.39243221690590112D+00, 8 1.34028640441683920D+01, -1.56848284626002017D+02, 9 2.19310333333333333D+03, -3.61087712537249894D+04, A 6.91472268851313067D+05, -1.52382215394074162D+07, B 3.82900751391414141D+08, -1.08822660357843911D+10, C 3.47320283765002252D+11, -1.23696021422692745D+13, D 4.88788064793079335D+14, -2.13203339609193739D+16/ C C LN(2*PI) DATA CON / 1.83787706640934548D+00/ C C***FIRST EXECUTABLE STATEMENT DGAMLN IERR=0 IF (Z.LE.0.0D0) GO TO 70 IF (Z.GT.101.0D0) GO TO 10 NZ = Z FZ = Z - NZ IF (FZ.GT.0.0D0) GO TO 10 IF (NZ.GT.100) GO TO 10 DGAMLN = GLN(NZ) RETURN 10 CONTINUE WDTOL = D1MACH(4) WDTOL = MAX(WDTOL,0.5D-18) I1M = I1MACH(14) RLN = D1MACH(5)*I1M FLN = MIN(RLN,20.0D0) FLN = MAX(FLN,3.0D0) FLN = FLN - 3.0D0 ZM = 1.8000D0 + 0.3875D0*FLN MZ = ZM + 1 ZMIN = MZ ZDMY = Z ZINC = 0.0D0 IF (Z.GE.ZMIN) GO TO 20 ZINC = ZMIN - NZ ZDMY = Z + ZINC 20 CONTINUE ZP = 1.0D0/ZDMY T1 = CF(1)*ZP S = T1 IF (ZP.LT.WDTOL) GO TO 40 ZSQ = ZP*ZP TST = T1*WDTOL DO 30 K=2,22 ZP = ZP*ZSQ TRM = CF(K)*ZP IF (ABS(TRM).LT.TST) GO TO 40 S = S + TRM 30 CONTINUE 40 CONTINUE IF (ZINC.NE.0.0D0) GO TO 50 TLG = LOG(Z) DGAMLN = Z*(TLG-1.0D0) + 0.5D0*(CON-TLG) + S RETURN 50 CONTINUE ZP = 1.0D0 NZ = ZINC DO 60 I=1,NZ ZP = ZP*(Z+(I-1)) 60 CONTINUE TLG = LOG(ZDMY) DGAMLN = ZDMY*(TLG-1.0D0) - LOG(ZP) + 0.5D0*(CON-TLG) + S RETURN C C 70 CONTINUE DGAMLN = D1MACH(2) IERR=1 RETURN END DOUBLE PRECISION FUNCTION ZABS (ZR, ZI) C***BEGIN PROLOGUE ZABS C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESH, ZBESI, ZBESJ, ZBESK, ZBESY, ZAIRY and C ZBIRY C***LIBRARY SLATEC C***TYPE ALL (ZABS-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZABS COMPUTES THE ABSOLUTE VALUE OR MAGNITUDE OF A DOUBLE C PRECISION COMPLEX VARIABLE CMPLX(ZR,ZI) C C***SEE ALSO ZAIRY, ZBESH, ZBESI, ZBESJ, ZBESK, ZBESY, ZBIRY C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZABS DOUBLE PRECISION ZR, ZI, U, V, Q, S C***FIRST EXECUTABLE STATEMENT ZABS U = ABS(ZR) V = ABS(ZI) S = U + V C----------------------------------------------------------------------- C S*1.0D0 MAKES AN UNNORMALIZED UNDERFLOW ON CDC MACHINES INTO A C TRUE FLOATING ZERO C----------------------------------------------------------------------- S = S*1.0D+0 IF (S.EQ.0.0D+0) GO TO 20 IF (U.GT.V) GO TO 10 Q = U/V ZABS = V*SQRT(1.D+0+Q*Q) RETURN 10 Q = V/U ZABS = U*SQRT(1.D+0+Q*Q) RETURN 20 ZABS = 0.0D+0 RETURN END SUBROUTINE ZACAI (ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, TOL, + ELIM, ALIM) C***BEGIN PROLOGUE ZACAI C***SUBSIDIARY C***PURPOSE Subsidiary to ZAIRY C***LIBRARY SLATEC C***TYPE ALL (CACAI-A, ZACAI-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZACAI APPLIES THE ANALYTIC CONTINUATION FORMULA C C K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN) C MP=PI*MR*CMPLX(0.0,1.0) C C TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT C HALF Z PLANE FOR USE WITH ZAIRY WHERE FNU=1/3 OR 2/3 AND N=1. C ZACAI IS THE SAME AS ZACON WITH THE PARTS FOR LARGER ORDERS AND C RECURRENCE REMOVED. A RECURSIVE CALL TO ZACON CAN RESULT IF ZACON C IS CALLED FROM ZAIRY. C C***SEE ALSO ZAIRY C***ROUTINES CALLED D1MACH, ZABS, ZASYI, ZBKNU, ZMLRI, ZS1S2, ZSERI C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZACAI C COMPLEX CSGN,CSPN,C1,C2,Y,Z,ZN,CY DOUBLE PRECISION ALIM, ARG, ASCLE, AZ, CSGNR, CSGNI, CSPNR, * CSPNI, C1R, C1I, C2R, C2I, CYR, CYI, DFNU, ELIM, FMR, FNU, PI, * RL, SGN, TOL, YY, YR, YI, ZR, ZI, ZNR, ZNI, D1MACH, ZABS INTEGER INU, IUF, KODE, MR, N, NN, NW, NZ DIMENSION YR(N), YI(N), CYR(2), CYI(2) EXTERNAL ZABS DATA PI / 3.14159265358979324D0 / C***FIRST EXECUTABLE STATEMENT ZACAI NZ = 0 ZNR = -ZR ZNI = -ZI AZ = ZABS(ZR,ZI) NN = N DFNU = FNU + (N-1) IF (AZ.LE.2.0D0) GO TO 10 IF (AZ*AZ*0.25D0.GT.DFNU+1.0D0) GO TO 20 10 CONTINUE C----------------------------------------------------------------------- C POWER SERIES FOR THE I FUNCTION C----------------------------------------------------------------------- CALL ZSERI(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL, ELIM, ALIM) GO TO 40 20 CONTINUE IF (AZ.LT.RL) GO TO 30 C----------------------------------------------------------------------- C ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION C----------------------------------------------------------------------- CALL ZASYI(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, TOL, ELIM, * ALIM) IF (NW.LT.0) GO TO 80 GO TO 40 30 CONTINUE C----------------------------------------------------------------------- C MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION C----------------------------------------------------------------------- CALL ZMLRI(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL) IF(NW.LT.0) GO TO 80 40 CONTINUE C----------------------------------------------------------------------- C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION C----------------------------------------------------------------------- CALL ZBKNU(ZNR, ZNI, FNU, KODE, 1, CYR, CYI, NW, TOL, ELIM, ALIM) IF (NW.NE.0) GO TO 80 FMR = MR SGN = -DSIGN(PI,FMR) CSGNR = 0.0D0 CSGNI = SGN IF (KODE.EQ.1) GO TO 50 YY = -ZNI CSGNR = -CSGNI*SIN(YY) CSGNI = CSGNI*COS(YY) 50 CONTINUE C----------------------------------------------------------------------- C CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE C WHEN FNU IS LARGE C----------------------------------------------------------------------- INU = FNU ARG = (FNU-INU)*SGN CSPNR = COS(ARG) CSPNI = SIN(ARG) IF (MOD(INU,2).EQ.0) GO TO 60 CSPNR = -CSPNR CSPNI = -CSPNI 60 CONTINUE C1R = CYR(1) C1I = CYI(1) C2R = YR(1) C2I = YI(1) IF (KODE.EQ.1) GO TO 70 IUF = 0 ASCLE = 1.0D+3*D1MACH(1)/TOL CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF) NZ = NZ + NW 70 CONTINUE YR(1) = CSPNR*C1R - CSPNI*C1I + CSGNR*C2R - CSGNI*C2I YI(1) = CSPNR*C1I + CSPNI*C1R + CSGNR*C2I + CSGNI*C2R RETURN 80 CONTINUE NZ = -1 IF(NW.EQ.(-2)) NZ=-2 RETURN END SUBROUTINE ZAIRY (ZR, ZI, ID, KODE, AIR, AII, NZ, IERR) C***BEGIN PROLOGUE ZAIRY C***PURPOSE Compute the Airy function Ai(z) or its derivative dAi/dz C for complex argument z. A scaling option is available C to help avoid underflow and overflow. C***LIBRARY SLATEC C***CATEGORY C10D C***TYPE COMPLEX (CAIRY-C, ZAIRY-C) C***KEYWORDS AIRY FUNCTION, BESSEL FUNCTION OF ORDER ONE THIRD, C BESSEL FUNCTION OF ORDER TWO THIRDS C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ***A DOUBLE PRECISION ROUTINE*** C On KODE=1, ZAIRY computes the complex Airy function Ai(z) C or its derivative dAi/dz on ID=0 or ID=1 respectively. On C KODE=2, a scaling option exp(zeta)*Ai(z) or exp(zeta)*dAi/dz C is provided to remove the exponential decay in -pi/31 and from power series when abs(z)<=1. C C In most complex variable computation, one must evaluate ele- C mentary functions. When the magnitude of Z is large, losses C of significance by argument reduction occur. Consequently, if C the magnitude of ZETA=(2/3)*Z**(3/2) exceeds U1=SQRT(0.5/UR), C then losses exceeding half precision are likely and an error C flag IERR=3 is triggered where UR=MAX(D1MACH(4),1.0D-18) is C double precision unit roundoff limited to 18 digits precision. C Also, if the magnitude of ZETA is larger than U2=0.5/UR, then C all significance is lost and IERR=4. In order to use the INT C function, ZETA must be further restricted not to exceed C U3=I1MACH(9)=LARGEST INTEGER. Thus, the magnitude of ZETA C must be restricted by MIN(U2,U3). In IEEE arithmetic, U1,U2, C and U3 are approximately 2.0E+3, 4.2E+6, 2.1E+9 in single C precision and 4.7E+7, 2.3E+15, 2.1E+9 in double precision. C This makes U2 limiting is single precision and U3 limiting C in double precision. This means that the magnitude of Z C cannot exceed approximately 3.4E+4 in single precision and C 2.1E+6 in double precision. This also means that one can C expect to retain, in the worst cases on 32-bit machines, C no digits in single precision and only 6 digits in double C precision. C C The approximate relative error in the magnitude of a complex C Bessel function can be expressed as P*10**S where P=MAX(UNIT C ROUNDOFF,1.0E-18) is the nominal precision and 10**S repre- C sents the increase in error due to argument reduction in the C elementary functions. Here, S=MAX(1,ABS(LOG10(ABS(Z))), C ABS(LOG10(FNU))) approximately (i.e., S=MAX(1,ABS(EXPONENT OF C ABS(Z),ABS(EXPONENT OF FNU)) ). However, the phase angle may C have only absolute accuracy. This is most likely to occur C when one component (in magnitude) is larger than the other by C several orders of magnitude. If one component is 10**K larger C than the other, then one can expect only MAX(ABS(LOG10(P))-K, C 0) significant digits; or, stated another way, when K exceeds C the exponent of P, no significant digits remain in the smaller C component. However, the phase angle retains absolute accuracy C because, in complex arithmetic with precision P, the smaller C component will not (as a rule) decrease below P times the C magnitude of the larger component. In these extreme cases, C the principal phase angle is on the order of +P, -P, PI/2-P, C or -PI/2+P. C C***REFERENCES 1. M. Abramowitz and I. A. Stegun, Handbook of Mathe- C matical Functions, National Bureau of Standards C Applied Mathematics Series 55, U. S. Department C of Commerce, Tenth Printing (1972) or later. C 2. D. E. Amos, Computation of Bessel Functions of C Complex Argument and Large Order, Report SAND83-0643, C Sandia National Laboratories, Albuquerque, NM, May C 1983. C 3. D. E. Amos, A Subroutine Package for Bessel Functions C of a Complex Argument and Nonnegative Order, Report C SAND85-1018, Sandia National Laboratory, Albuquerque, C NM, May 1985. C 4. D. E. Amos, A portable package for Bessel functions C of a complex argument and nonnegative order, ACM C Transactions on Mathematical Software, 12 (September C 1986), pp. 265-273. C C***ROUTINES CALLED D1MACH, I1MACH, ZABS, ZACAI, ZBKNU, ZEXP, ZSQRT C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 890801 REVISION DATE from Version 3.2 C 910415 Prologue converted to Version 4.0 format. (BAB) C 920128 Category corrected. (WRB) C 920811 Prologue revised. (DWL) C***END PROLOGUE ZAIRY C COMPLEX AI,CONE,CSQ,CY,S1,S2,TRM1,TRM2,Z,ZTA,Z3 DOUBLE PRECISION AA, AD, AII, AIR, AK, ALIM, ATRM, AZ, AZ3, BK, * CC, CK, COEF, CONEI, CONER, CSQI, CSQR, CYI, CYR, C1, C2, DIG, * DK, D1, D2, ELIM, FID, FNU, PTR, RL, R1M5, SFAC, STI, STR, * S1I, S1R, S2I, S2R, TOL, TRM1I, TRM1R, TRM2I, TRM2R, TTH, ZEROI, * ZEROR, ZI, ZR, ZTAI, ZTAR, Z3I, Z3R, D1MACH, ZABS, ALAZ, BB INTEGER ID, IERR, IFLAG, K, KODE, K1, K2, MR, NN, NZ, I1MACH DIMENSION CYR(1), CYI(1) EXTERNAL ZABS DATA TTH, C1, C2, COEF /6.66666666666666667D-01, * 3.55028053887817240D-01,2.58819403792806799D-01, * 1.83776298473930683D-01/ DATA ZEROR, ZEROI, CONER, CONEI /0.0D0,0.0D0,1.0D0,0.0D0/ C***FIRST EXECUTABLE STATEMENT ZAIRY IERR = 0 NZ=0 IF (ID.LT.0 .OR. ID.GT.1) IERR=1 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 IF (IERR.NE.0) RETURN AZ = ZABS(ZR,ZI) TOL = MAX(D1MACH(4),1.0D-18) FID = ID IF (AZ.GT.1.0D0) GO TO 70 C----------------------------------------------------------------------- C POWER SERIES FOR ABS(Z).LE.1. C----------------------------------------------------------------------- S1R = CONER S1I = CONEI S2R = CONER S2I = CONEI IF (AZ.LT.TOL) GO TO 170 AA = AZ*AZ IF (AA.LT.TOL/AZ) GO TO 40 TRM1R = CONER TRM1I = CONEI TRM2R = CONER TRM2I = CONEI ATRM = 1.0D0 STR = ZR*ZR - ZI*ZI STI = ZR*ZI + ZI*ZR Z3R = STR*ZR - STI*ZI Z3I = STR*ZI + STI*ZR AZ3 = AZ*AA AK = 2.0D0 + FID BK = 3.0D0 - FID - FID CK = 4.0D0 - FID DK = 3.0D0 + FID + FID D1 = AK*DK D2 = BK*CK AD = MIN(D1,D2) AK = 24.0D0 + 9.0D0*FID BK = 30.0D0 - 9.0D0*FID DO 30 K=1,25 STR = (TRM1R*Z3R-TRM1I*Z3I)/D1 TRM1I = (TRM1R*Z3I+TRM1I*Z3R)/D1 TRM1R = STR S1R = S1R + TRM1R S1I = S1I + TRM1I STR = (TRM2R*Z3R-TRM2I*Z3I)/D2 TRM2I = (TRM2R*Z3I+TRM2I*Z3R)/D2 TRM2R = STR S2R = S2R + TRM2R S2I = S2I + TRM2I ATRM = ATRM*AZ3/AD D1 = D1 + AK D2 = D2 + BK AD = MIN(D1,D2) IF (ATRM.LT.TOL*AD) GO TO 40 AK = AK + 18.0D0 BK = BK + 18.0D0 30 CONTINUE 40 CONTINUE IF (ID.EQ.1) GO TO 50 AIR = S1R*C1 - C2*(ZR*S2R-ZI*S2I) AII = S1I*C1 - C2*(ZR*S2I+ZI*S2R) IF (KODE.EQ.1) RETURN CALL ZSQRT(ZR, ZI, STR, STI) ZTAR = TTH*(ZR*STR-ZI*STI) ZTAI = TTH*(ZR*STI+ZI*STR) CALL ZEXP(ZTAR, ZTAI, STR, STI) PTR = AIR*STR - AII*STI AII = AIR*STI + AII*STR AIR = PTR RETURN 50 CONTINUE AIR = -S2R*C2 AII = -S2I*C2 IF (AZ.LE.TOL) GO TO 60 STR = ZR*S1R - ZI*S1I STI = ZR*S1I + ZI*S1R CC = C1/(1.0D0+FID) AIR = AIR + CC*(STR*ZR-STI*ZI) AII = AII + CC*(STR*ZI+STI*ZR) 60 CONTINUE IF (KODE.EQ.1) RETURN CALL ZSQRT(ZR, ZI, STR, STI) ZTAR = TTH*(ZR*STR-ZI*STI) ZTAI = TTH*(ZR*STI+ZI*STR) CALL ZEXP(ZTAR, ZTAI, STR, STI) PTR = STR*AIR - STI*AII AII = STR*AII + STI*AIR AIR = PTR RETURN C----------------------------------------------------------------------- C CASE FOR ABS(Z).GT.1.0 C----------------------------------------------------------------------- 70 CONTINUE FNU = (1.0D0+FID)/3.0D0 C----------------------------------------------------------------------- C SET PARAMETERS RELATED TO MACHINE CONSTANTS. C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0D-18. C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). C----------------------------------------------------------------------- K1 = I1MACH(15) K2 = I1MACH(16) R1M5 = D1MACH(5) K = MIN(ABS(K1),ABS(K2)) ELIM = 2.303D0*(K*R1M5-3.0D0) K1 = I1MACH(14) - 1 AA = R1M5*K1 DIG = MIN(AA,18.0D0) AA = AA*2.303D0 ALIM = ELIM + MAX(-AA,-41.45D0) RL = 1.2D0*DIG + 3.0D0 ALAZ = LOG(AZ) C----------------------------------------------------------------------- C TEST FOR PROPER RANGE C----------------------------------------------------------------------- AA=0.5D0/TOL BB=I1MACH(9)*0.5D0 AA=MIN(AA,BB) AA=AA**TTH IF (AZ.GT.AA) GO TO 260 AA=SQRT(AA) IF (AZ.GT.AA) IERR=3 CALL ZSQRT(ZR, ZI, CSQR, CSQI) ZTAR = TTH*(ZR*CSQR-ZI*CSQI) ZTAI = TTH*(ZR*CSQI+ZI*CSQR) C----------------------------------------------------------------------- C RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL C----------------------------------------------------------------------- IFLAG = 0 SFAC = 1.0D0 AK = ZTAI IF (ZR.GE.0.0D0) GO TO 80 BK = ZTAR CK = -ABS(BK) ZTAR = CK ZTAI = AK 80 CONTINUE IF (ZI.NE.0.0D0) GO TO 90 IF (ZR.GT.0.0D0) GO TO 90 ZTAR = 0.0D0 ZTAI = AK 90 CONTINUE AA = ZTAR IF (AA.GE.0.0D0 .AND. ZR.GT.0.0D0) GO TO 110 IF (KODE.EQ.2) GO TO 100 C----------------------------------------------------------------------- C OVERFLOW TEST C----------------------------------------------------------------------- IF (AA.GT.(-ALIM)) GO TO 100 AA = -AA + 0.25D0*ALAZ IFLAG = 1 SFAC = TOL IF (AA.GT.ELIM) GO TO 270 100 CONTINUE C----------------------------------------------------------------------- C CBKNU AND CACON RETURN EXP(ZTA)*K(FNU,ZTA) ON KODE=2 C----------------------------------------------------------------------- MR = 1 IF (ZI.LT.0.0D0) MR = -1 CALL ZACAI(ZTAR, ZTAI, FNU, KODE, MR, 1, CYR, CYI, NN, RL, TOL, * ELIM, ALIM) IF (NN.LT.0) GO TO 280 NZ = NZ + NN GO TO 130 110 CONTINUE IF (KODE.EQ.2) GO TO 120 C----------------------------------------------------------------------- C UNDERFLOW TEST C----------------------------------------------------------------------- IF (AA.LT.ALIM) GO TO 120 AA = -AA - 0.25D0*ALAZ IFLAG = 2 SFAC = 1.0D0/TOL IF (AA.LT.(-ELIM)) GO TO 210 120 CONTINUE CALL ZBKNU(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, NZ, TOL, ELIM, * ALIM) 130 CONTINUE S1R = CYR(1)*COEF S1I = CYI(1)*COEF IF (IFLAG.NE.0) GO TO 150 IF (ID.EQ.1) GO TO 140 AIR = CSQR*S1R - CSQI*S1I AII = CSQR*S1I + CSQI*S1R RETURN 140 CONTINUE AIR = -(ZR*S1R-ZI*S1I) AII = -(ZR*S1I+ZI*S1R) RETURN 150 CONTINUE S1R = S1R*SFAC S1I = S1I*SFAC IF (ID.EQ.1) GO TO 160 STR = S1R*CSQR - S1I*CSQI S1I = S1R*CSQI + S1I*CSQR S1R = STR AIR = S1R/SFAC AII = S1I/SFAC RETURN 160 CONTINUE STR = -(S1R*ZR-S1I*ZI) S1I = -(S1R*ZI+S1I*ZR) S1R = STR AIR = S1R/SFAC AII = S1I/SFAC RETURN 170 CONTINUE AA = 1.0D+3*D1MACH(1) S1R = ZEROR S1I = ZEROI IF (ID.EQ.1) GO TO 190 IF (AZ.LE.AA) GO TO 180 S1R = C2*ZR S1I = C2*ZI 180 CONTINUE AIR = C1 - S1R AII = -S1I RETURN 190 CONTINUE AIR = -C2 AII = 0.0D0 AA = SQRT(AA) IF (AZ.LE.AA) GO TO 200 S1R = 0.5D0*(ZR*ZR-ZI*ZI) S1I = ZR*ZI 200 CONTINUE AIR = AIR + C1*S1R AII = AII + C1*S1I RETURN 210 CONTINUE NZ = 1 AIR = ZEROR AII = ZEROI RETURN 270 CONTINUE NZ = 0 IERR=2 RETURN 280 CONTINUE IF(NN.EQ.(-1)) GO TO 270 NZ=0 IERR=5 RETURN 260 CONTINUE IERR=4 NZ=0 RETURN END SUBROUTINE ZASYI (ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM, + ALIM) C***BEGIN PROLOGUE ZASYI C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESI and ZBESK C***LIBRARY SLATEC C***TYPE ALL (CASYI-A, ZASYI-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY C MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE ABS(Z) IN THE C REGION ABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN. C NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1. C C***SEE ALSO ZBESI, ZBESK C***ROUTINES CALLED D1MACH, ZABS, ZDIV, ZEXP, ZMLT, ZSQRT C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZASYI C COMPLEX AK1,CK,CONE,CS1,CS2,CZ,CZERO,DK,EZ,P1,RZ,S2,Y,Z DOUBLE PRECISION AA, AEZ, AK, AK1I, AK1R, ALIM, ARG, ARM, ATOL, * AZ, BB, BK, CKI, CKR, CONEI, CONER, CS1I, CS1R, CS2I, CS2R, CZI, * CZR, DFNU, DKI, DKR, DNU2, ELIM, EZI, EZR, FDN, FNU, PI, P1I, * P1R, RAZ, RL, RTPI, RTR1, RZI, RZR, S, SGN, SQK, STI, STR, S2I, * S2R, TOL, TZI, TZR, YI, YR, ZEROI, ZEROR, ZI, ZR, D1MACH, ZABS INTEGER I, IB, IL, INU, J, JL, K, KODE, KODED, M, N, NN, NZ DIMENSION YR(N), YI(N) EXTERNAL ZABS DATA PI, RTPI /3.14159265358979324D0 , 0.159154943091895336D0 / DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 / C***FIRST EXECUTABLE STATEMENT ZASYI NZ = 0 AZ = ZABS(ZR,ZI) ARM = 1.0D+3*D1MACH(1) RTR1 = SQRT(ARM) IL = MIN(2,N) DFNU = FNU + (N-IL) C----------------------------------------------------------------------- C OVERFLOW TEST C----------------------------------------------------------------------- RAZ = 1.0D0/AZ STR = ZR*RAZ STI = -ZI*RAZ AK1R = RTPI*STR*RAZ AK1I = RTPI*STI*RAZ CALL ZSQRT(AK1R, AK1I, AK1R, AK1I) CZR = ZR CZI = ZI IF (KODE.NE.2) GO TO 10 CZR = ZEROR CZI = ZI 10 CONTINUE IF (ABS(CZR).GT.ELIM) GO TO 100 DNU2 = DFNU + DFNU KODED = 1 IF ((ABS(CZR).GT.ALIM) .AND. (N.GT.2)) GO TO 20 KODED = 0 CALL ZEXP(CZR, CZI, STR, STI) CALL ZMLT(AK1R, AK1I, STR, STI, AK1R, AK1I) 20 CONTINUE FDN = 0.0D0 IF (DNU2.GT.RTR1) FDN = DNU2*DNU2 EZR = ZR*8.0D0 EZI = ZI*8.0D0 C----------------------------------------------------------------------- C WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE C FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE C EXPANSION FOR THE IMAGINARY PART. C----------------------------------------------------------------------- AEZ = 8.0D0*AZ S = TOL/AEZ JL = RL+RL + 2 P1R = ZEROR P1I = ZEROI IF (ZI.EQ.0.0D0) GO TO 30 C----------------------------------------------------------------------- C CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF C SIGNIFICANCE WHEN FNU OR N IS LARGE C----------------------------------------------------------------------- INU = FNU ARG = (FNU-INU)*PI INU = INU + N - IL AK = -SIN(ARG) BK = COS(ARG) IF (ZI.LT.0.0D0) BK = -BK P1R = AK P1I = BK IF (MOD(INU,2).EQ.0) GO TO 30 P1R = -P1R P1I = -P1I 30 CONTINUE DO 70 K=1,IL SQK = FDN - 1.0D0 ATOL = S*ABS(SQK) SGN = 1.0D0 CS1R = CONER CS1I = CONEI CS2R = CONER CS2I = CONEI CKR = CONER CKI = CONEI AK = 0.0D0 AA = 1.0D0 BB = AEZ DKR = EZR DKI = EZI DO 40 J=1,JL CALL ZDIV(CKR, CKI, DKR, DKI, STR, STI) CKR = STR*SQK CKI = STI*SQK CS2R = CS2R + CKR CS2I = CS2I + CKI SGN = -SGN CS1R = CS1R + CKR*SGN CS1I = CS1I + CKI*SGN DKR = DKR + EZR DKI = DKI + EZI AA = AA*ABS(SQK)/BB BB = BB + AEZ AK = AK + 8.0D0 SQK = SQK - AK IF (AA.LE.ATOL) GO TO 50 40 CONTINUE GO TO 110 50 CONTINUE S2R = CS1R S2I = CS1I IF (ZR+ZR.GE.ELIM) GO TO 60 TZR = ZR + ZR TZI = ZI + ZI CALL ZEXP(-TZR, -TZI, STR, STI) CALL ZMLT(STR, STI, P1R, P1I, STR, STI) CALL ZMLT(STR, STI, CS2R, CS2I, STR, STI) S2R = S2R + STR S2I = S2I + STI 60 CONTINUE FDN = FDN + 8.0D0*DFNU + 4.0D0 P1R = -P1R P1I = -P1I M = N - IL + K YR(M) = S2R*AK1R - S2I*AK1I YI(M) = S2R*AK1I + S2I*AK1R 70 CONTINUE IF (N.LE.2) RETURN NN = N K = NN - 2 AK = K STR = ZR*RAZ STI = -ZI*RAZ RZR = (STR+STR)*RAZ RZI = (STI+STI)*RAZ IB = 3 DO 80 I=IB,NN YR(K) = (AK+FNU)*(RZR*YR(K+1)-RZI*YI(K+1)) + YR(K+2) YI(K) = (AK+FNU)*(RZR*YI(K+1)+RZI*YR(K+1)) + YI(K+2) AK = AK - 1.0D0 K = K - 1 80 CONTINUE IF (KODED.EQ.0) RETURN CALL ZEXP(CZR, CZI, CKR, CKI) DO 90 I=1,NN STR = YR(I)*CKR - YI(I)*CKI YI(I) = YR(I)*CKI + YI(I)*CKR YR(I) = STR 90 CONTINUE RETURN 100 CONTINUE NZ = -1 RETURN 110 CONTINUE NZ=-2 RETURN END SUBROUTINE ZBINU (ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL, + TOL, ELIM, ALIM) C***BEGIN PROLOGUE ZBINU C***SUBSIDIARY C***PURPOSE Subsidiary to ZAIRY, ZBESH, ZBESI, ZBESJ, ZBESK and ZBIRY C***LIBRARY SLATEC C***TYPE ALL (CBINU-A, ZBINU-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZBINU COMPUTES THE I FUNCTION IN THE RIGHT HALF Z PLANE C C***SEE ALSO ZAIRY, ZBESH, ZBESI, ZBESJ, ZBESK, ZBIRY C***ROUTINES CALLED ZABS, ZASYI, ZBUNI, ZMLRI, ZSERI, ZUOIK, ZWRSK C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZBINU DOUBLE PRECISION ALIM, AZ, CWI, CWR, CYI, CYR, DFNU, ELIM, FNU, * FNUL, RL, TOL, ZEROI, ZEROR, ZI, ZR, ZABS INTEGER I, INW, KODE, N, NLAST, NN, NUI, NW, NZ DIMENSION CYR(N), CYI(N), CWR(2), CWI(2) EXTERNAL ZABS DATA ZEROR,ZEROI / 0.0D0, 0.0D0 / C***FIRST EXECUTABLE STATEMENT ZBINU NZ = 0 AZ = ZABS(ZR,ZI) NN = N DFNU = FNU + (N-1) IF (AZ.LE.2.0D0) GO TO 10 IF (AZ*AZ*0.25D0.GT.DFNU+1.0D0) GO TO 20 10 CONTINUE C----------------------------------------------------------------------- C POWER SERIES C----------------------------------------------------------------------- CALL ZSERI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, TOL, ELIM, ALIM) INW = ABS(NW) NZ = NZ + INW NN = NN - INW IF (NN.EQ.0) RETURN IF (NW.GE.0) GO TO 120 DFNU = FNU + (NN-1) 20 CONTINUE IF (AZ.LT.RL) GO TO 40 IF (DFNU.LE.1.0D0) GO TO 30 IF (AZ+AZ.LT.DFNU*DFNU) GO TO 50 C----------------------------------------------------------------------- C ASYMPTOTIC EXPANSION FOR LARGE Z C----------------------------------------------------------------------- 30 CONTINUE CALL ZASYI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, RL, TOL, ELIM, * ALIM) IF (NW.LT.0) GO TO 130 GO TO 120 40 CONTINUE IF (DFNU.LE.1.0D0) GO TO 70 50 CONTINUE C----------------------------------------------------------------------- C OVERFLOW AND UNDERFLOW TEST ON I SEQUENCE FOR MILLER ALGORITHM C----------------------------------------------------------------------- CALL ZUOIK(ZR, ZI, FNU, KODE, 1, NN, CYR, CYI, NW, TOL, ELIM, * ALIM) IF (NW.LT.0) GO TO 130 NZ = NZ + NW NN = NN - NW IF (NN.EQ.0) RETURN DFNU = FNU+(NN-1) IF (DFNU.GT.FNUL) GO TO 110 IF (AZ.GT.FNUL) GO TO 110 60 CONTINUE IF (AZ.GT.RL) GO TO 80 70 CONTINUE C----------------------------------------------------------------------- C MILLER ALGORITHM NORMALIZED BY THE SERIES C----------------------------------------------------------------------- CALL ZMLRI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, TOL) IF(NW.LT.0) GO TO 130 GO TO 120 80 CONTINUE C----------------------------------------------------------------------- C MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN C----------------------------------------------------------------------- C----------------------------------------------------------------------- C OVERFLOW TEST ON K FUNCTIONS USED IN WRONSKIAN C----------------------------------------------------------------------- CALL ZUOIK(ZR, ZI, FNU, KODE, 2, 2, CWR, CWI, NW, TOL, ELIM, * ALIM) IF (NW.GE.0) GO TO 100 NZ = NN DO 90 I=1,NN CYR(I) = ZEROR CYI(I) = ZEROI 90 CONTINUE RETURN 100 CONTINUE IF (NW.GT.0) GO TO 130 CALL ZWRSK(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, CWR, CWI, TOL, * ELIM, ALIM) IF (NW.LT.0) GO TO 130 GO TO 120 110 CONTINUE C----------------------------------------------------------------------- C INCREMENT FNU+NN-1 UP TO FNUL, COMPUTE AND RECUR BACKWARD C----------------------------------------------------------------------- NUI = FNUL-DFNU + 1 NUI = MAX(NUI,0) CALL ZBUNI(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, NUI, NLAST, FNUL, * TOL, ELIM, ALIM) IF (NW.LT.0) GO TO 130 NZ = NZ + NW IF (NLAST.EQ.0) GO TO 120 NN = NLAST GO TO 60 120 CONTINUE RETURN 130 CONTINUE NZ = -1 IF(NW.EQ.(-2)) NZ=-2 RETURN END SUBROUTINE ZBIRY (ZR, ZI, ID, KODE, BIR, BII, IERR) C***BEGIN PROLOGUE ZBIRY C***PURPOSE Compute the Airy function Bi(z) or its derivative dBi/dz C for complex argument z. A scaling option is available C to help avoid overflow. C***LIBRARY SLATEC C***CATEGORY C10D C***TYPE COMPLEX (CBIRY-C, ZBIRY-C) C***KEYWORDS AIRY FUNCTION, BESSEL FUNCTION OF ORDER ONE THIRD, C BESSEL FUNCTION OF ORDER TWO THIRDS C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ***A DOUBLE PRECISION ROUTINE*** C On KODE=1, ZBIRY computes the complex Airy function Bi(z) C or its derivative dBi/dz on ID=0 or ID=1 respectively. C On KODE=2, a scaling option exp(abs(Re(zeta)))*Bi(z) or C exp(abs(Re(zeta)))*dBi/dz is provided to remove the C exponential behavior in both the left and right half planes C where zeta=(2/3)*z**(3/2). C C The Airy functions Bi(z) and dBi/dz are analytic in the C whole z-plane, and the scaling option does not destroy this C property. C C Input C ZR - DOUBLE PRECISION real part of argument Z C ZI - DOUBLE PRECISION imag part of argument Z C ID - Order of derivative, ID=0 or ID=1 C KODE - A parameter to indicate the scaling option C KODE=1 returns C BI=Bi(z) on ID=0 C BI=dBi/dz on ID=1 C at z=Z C =2 returns C BI=exp(abs(Re(zeta)))*Bi(z) on ID=0 C BI=exp(abs(Re(zeta)))*dBi/dz on ID=1 C at z=Z where zeta=(2/3)*z**(3/2) C C Output C BIR - DOUBLE PRECISION real part of result C BII - DOUBLE PRECISION imag part of result C IERR - Error flag C IERR=0 Normal return - COMPUTATION COMPLETED C IERR=1 Input error - NO COMPUTATION C IERR=2 Overflow - NO COMPUTATION C (Re(Z) too large with KODE=1) C IERR=3 Precision warning - COMPUTATION COMPLETED C (Result has less than half precision) C IERR=4 Precision error - NO COMPUTATION C (Result has no precision) C IERR=5 Algorithmic error - NO COMPUTATION C (Termination condition not met) C C *Long Description: C C Bi(z) and dBi/dz are computed from I Bessel functions by C C Bi(z) = c*sqrt(z)*( I(-1/3,zeta) + I(1/3,zeta) ) C dBi/dz = c* z *( I(-2/3,zeta) + I(2/3,zeta) ) C c = 1/sqrt(3) C zeta = (2/3)*z**(3/2) C C when abs(z)>1 and from power series when abs(z)<=1. C C In most complex variable computation, one must evaluate ele- C mentary functions. When the magnitude of Z is large, losses C of significance by argument reduction occur. Consequently, if C the magnitude of ZETA=(2/3)*Z**(3/2) exceeds U1=SQRT(0.5/UR), C then losses exceeding half precision are likely and an error C flag IERR=3 is triggered where UR=MAX(D1MACH(4),1.0D-18) is C double precision unit roundoff limited to 18 digits precision. C Also, if the magnitude of ZETA is larger than U2=0.5/UR, then C all significance is lost and IERR=4. In order to use the INT C function, ZETA must be further restricted not to exceed C U3=I1MACH(9)=LARGEST INTEGER. Thus, the magnitude of ZETA C must be restricted by MIN(U2,U3). In IEEE arithmetic, U1,U2, C and U3 are approximately 2.0E+3, 4.2E+6, 2.1E+9 in single C precision and 4.7E+7, 2.3E+15, 2.1E+9 in double precision. C This makes U2 limiting is single precision and U3 limiting C in double precision. This means that the magnitude of Z C cannot exceed approximately 3.4E+4 in single precision and C 2.1E+6 in double precision. This also means that one can C expect to retain, in the worst cases on 32-bit machines, C no digits in single precision and only 6 digits in double C precision. C C The approximate relative error in the magnitude of a complex C Bessel function can be expressed as P*10**S where P=MAX(UNIT C ROUNDOFF,1.0E-18) is the nominal precision and 10**S repre- C sents the increase in error due to argument reduction in the C elementary functions. Here, S=MAX(1,ABS(LOG10(ABS(Z))), C ABS(LOG10(FNU))) approximately (i.e., S=MAX(1,ABS(EXPONENT OF C ABS(Z),ABS(EXPONENT OF FNU)) ). However, the phase angle may C have only absolute accuracy. This is most likely to occur C when one component (in magnitude) is larger than the other by C several orders of magnitude. If one component is 10**K larger C than the other, then one can expect only MAX(ABS(LOG10(P))-K, C 0) significant digits; or, stated another way, when K exceeds C the exponent of P, no significant digits remain in the smaller C component. However, the phase angle retains absolute accuracy C because, in complex arithmetic with precision P, the smaller C component will not (as a rule) decrease below P times the C magnitude of the larger component. In these extreme cases, C the principal phase angle is on the order of +P, -P, PI/2-P, C or -PI/2+P. C C***REFERENCES 1. M. Abramowitz and I. A. Stegun, Handbook of Mathe- C matical Functions, National Bureau of Standards C Applied Mathematics Series 55, U. S. Department C of Commerce, Tenth Printing (1972) or later. C 2. D. E. Amos, Computation of Bessel Functions of C Complex Argument and Large Order, Report SAND83-0643, C Sandia National Laboratories, Albuquerque, NM, May C 1983. C 3. D. E. Amos, A Subroutine Package for Bessel Functions C of a Complex Argument and Nonnegative Order, Report C SAND85-1018, Sandia National Laboratory, Albuquerque, C NM, May 1985. C 4. D. E. Amos, A portable package for Bessel functions C of a complex argument and nonnegative order, ACM C Transactions on Mathematical Software, 12 (September C 1986), pp. 265-273. C C***ROUTINES CALLED D1MACH, I1MACH, ZABS, ZBINU, ZDIV, ZSQRT C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 890801 REVISION DATE from Version 3.2 C 910415 Prologue converted to Version 4.0 format. (BAB) C 920128 Category corrected. (WRB) C 920811 Prologue revised. (DWL) C***END PROLOGUE ZBIRY C COMPLEX BI,CONE,CSQ,CY,S1,S2,TRM1,TRM2,Z,ZTA,Z3 DOUBLE PRECISION AA, AD, AK, ALIM, ATRM, AZ, AZ3, BB, BII, BIR, * BK, CC, CK, COEF, CONEI, CONER, CSQI, CSQR, CYI, CYR, C1, C2, * DIG, DK, D1, D2, EAA, ELIM, FID, FMR, FNU, FNUL, PI, RL, R1M5, * SFAC, STI, STR, S1I, S1R, S2I, S2R, TOL, TRM1I, TRM1R, TRM2I, * TRM2R, TTH, ZI, ZR, ZTAI, ZTAR, Z3I, Z3R, D1MACH, ZABS INTEGER ID, IERR, K, KODE, K1, K2, NZ, I1MACH DIMENSION CYR(2), CYI(2) EXTERNAL ZABS DATA TTH, C1, C2, COEF, PI /6.66666666666666667D-01, * 6.14926627446000736D-01,4.48288357353826359D-01, * 5.77350269189625765D-01,3.14159265358979324D+00/ DATA CONER, CONEI /1.0D0,0.0D0/ C***FIRST EXECUTABLE STATEMENT ZBIRY IERR = 0 NZ=0 IF (ID.LT.0 .OR. ID.GT.1) IERR=1 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 IF (IERR.NE.0) RETURN AZ = ZABS(ZR,ZI) TOL = MAX(D1MACH(4),1.0D-18) FID = ID IF (AZ.GT.1.0E0) GO TO 70 C----------------------------------------------------------------------- C POWER SERIES FOR ABS(Z).LE.1. C----------------------------------------------------------------------- S1R = CONER S1I = CONEI S2R = CONER S2I = CONEI IF (AZ.LT.TOL) GO TO 130 AA = AZ*AZ IF (AA.LT.TOL/AZ) GO TO 40 TRM1R = CONER TRM1I = CONEI TRM2R = CONER TRM2I = CONEI ATRM = 1.0D0 STR = ZR*ZR - ZI*ZI STI = ZR*ZI + ZI*ZR Z3R = STR*ZR - STI*ZI Z3I = STR*ZI + STI*ZR AZ3 = AZ*AA AK = 2.0D0 + FID BK = 3.0D0 - FID - FID CK = 4.0D0 - FID DK = 3.0D0 + FID + FID D1 = AK*DK D2 = BK*CK AD = MIN(D1,D2) AK = 24.0D0 + 9.0D0*FID BK = 30.0D0 - 9.0D0*FID DO 30 K=1,25 STR = (TRM1R*Z3R-TRM1I*Z3I)/D1 TRM1I = (TRM1R*Z3I+TRM1I*Z3R)/D1 TRM1R = STR S1R = S1R + TRM1R S1I = S1I + TRM1I STR = (TRM2R*Z3R-TRM2I*Z3I)/D2 TRM2I = (TRM2R*Z3I+TRM2I*Z3R)/D2 TRM2R = STR S2R = S2R + TRM2R S2I = S2I + TRM2I ATRM = ATRM*AZ3/AD D1 = D1 + AK D2 = D2 + BK AD = MIN(D1,D2) IF (ATRM.LT.TOL*AD) GO TO 40 AK = AK + 18.0D0 BK = BK + 18.0D0 30 CONTINUE 40 CONTINUE IF (ID.EQ.1) GO TO 50 BIR = C1*S1R + C2*(ZR*S2R-ZI*S2I) BII = C1*S1I + C2*(ZR*S2I+ZI*S2R) IF (KODE.EQ.1) RETURN CALL ZSQRT(ZR, ZI, STR, STI) ZTAR = TTH*(ZR*STR-ZI*STI) ZTAI = TTH*(ZR*STI+ZI*STR) AA = ZTAR AA = -ABS(AA) EAA = EXP(AA) BIR = BIR*EAA BII = BII*EAA RETURN 50 CONTINUE BIR = S2R*C2 BII = S2I*C2 IF (AZ.LE.TOL) GO TO 60 CC = C1/(1.0D0+FID) STR = S1R*ZR - S1I*ZI STI = S1R*ZI + S1I*ZR BIR = BIR + CC*(STR*ZR-STI*ZI) BII = BII + CC*(STR*ZI+STI*ZR) 60 CONTINUE IF (KODE.EQ.1) RETURN CALL ZSQRT(ZR, ZI, STR, STI) ZTAR = TTH*(ZR*STR-ZI*STI) ZTAI = TTH*(ZR*STI+ZI*STR) AA = ZTAR AA = -ABS(AA) EAA = EXP(AA) BIR = BIR*EAA BII = BII*EAA RETURN C----------------------------------------------------------------------- C CASE FOR ABS(Z).GT.1.0 C----------------------------------------------------------------------- 70 CONTINUE FNU = (1.0D0+FID)/3.0D0 C----------------------------------------------------------------------- C SET PARAMETERS RELATED TO MACHINE CONSTANTS. C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU. C----------------------------------------------------------------------- K1 = I1MACH(15) K2 = I1MACH(16) R1M5 = D1MACH(5) K = MIN(ABS(K1),ABS(K2)) ELIM = 2.303D0*(K*R1M5-3.0D0) K1 = I1MACH(14) - 1 AA = R1M5*K1 DIG = MIN(AA,18.0D0) AA = AA*2.303D0 ALIM = ELIM + MAX(-AA,-41.45D0) RL = 1.2D0*DIG + 3.0D0 FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) C----------------------------------------------------------------------- C TEST FOR RANGE C----------------------------------------------------------------------- AA=0.5D0/TOL BB=I1MACH(9)*0.5D0 AA=MIN(AA,BB) AA=AA**TTH IF (AZ.GT.AA) GO TO 260 AA=SQRT(AA) IF (AZ.GT.AA) IERR=3 CALL ZSQRT(ZR, ZI, CSQR, CSQI) ZTAR = TTH*(ZR*CSQR-ZI*CSQI) ZTAI = TTH*(ZR*CSQI+ZI*CSQR) C----------------------------------------------------------------------- C RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL C----------------------------------------------------------------------- SFAC = 1.0D0 AK = ZTAI IF (ZR.GE.0.0D0) GO TO 80 BK = ZTAR CK = -ABS(BK) ZTAR = CK ZTAI = AK 80 CONTINUE IF (ZI.NE.0.0D0 .OR. ZR.GT.0.0D0) GO TO 90 ZTAR = 0.0D0 ZTAI = AK 90 CONTINUE AA = ZTAR IF (KODE.EQ.2) GO TO 100 C----------------------------------------------------------------------- C OVERFLOW TEST C----------------------------------------------------------------------- BB = ABS(AA) IF (BB.LT.ALIM) GO TO 100 BB = BB + 0.25D0*LOG(AZ) SFAC = TOL IF (BB.GT.ELIM) GO TO 190 100 CONTINUE FMR = 0.0D0 IF (AA.GE.0.0D0 .AND. ZR.GT.0.0D0) GO TO 110 FMR = PI IF (ZI.LT.0.0D0) FMR = -PI ZTAR = -ZTAR ZTAI = -ZTAI 110 CONTINUE C----------------------------------------------------------------------- C AA=FACTOR FOR ANALYTIC CONTINUATION OF I(FNU,ZTA) C KODE=2 RETURNS EXP(-ABS(XZTA))*I(FNU,ZTA) FROM CBESI C----------------------------------------------------------------------- CALL ZBINU(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, NZ, RL, FNUL, TOL, * ELIM, ALIM) IF (NZ.LT.0) GO TO 200 AA = FMR*FNU Z3R = SFAC STR = COS(AA) STI = SIN(AA) S1R = (STR*CYR(1)-STI*CYI(1))*Z3R S1I = (STR*CYI(1)+STI*CYR(1))*Z3R FNU = (2.0D0-FID)/3.0D0 CALL ZBINU(ZTAR, ZTAI, FNU, KODE, 2, CYR, CYI, NZ, RL, FNUL, TOL, * ELIM, ALIM) CYR(1) = CYR(1)*Z3R CYI(1) = CYI(1)*Z3R CYR(2) = CYR(2)*Z3R CYI(2) = CYI(2)*Z3R C----------------------------------------------------------------------- C BACKWARD RECUR ONE STEP FOR ORDERS -1/3 OR -2/3 C----------------------------------------------------------------------- CALL ZDIV(CYR(1), CYI(1), ZTAR, ZTAI, STR, STI) S2R = (FNU+FNU)*STR + CYR(2) S2I = (FNU+FNU)*STI + CYI(2) AA = FMR*(FNU-1.0D0) STR = COS(AA) STI = SIN(AA) S1R = COEF*(S1R+S2R*STR-S2I*STI) S1I = COEF*(S1I+S2R*STI+S2I*STR) IF (ID.EQ.1) GO TO 120 STR = CSQR*S1R - CSQI*S1I S1I = CSQR*S1I + CSQI*S1R S1R = STR BIR = S1R/SFAC BII = S1I/SFAC RETURN 120 CONTINUE STR = ZR*S1R - ZI*S1I S1I = ZR*S1I + ZI*S1R S1R = STR BIR = S1R/SFAC BII = S1I/SFAC RETURN 130 CONTINUE AA = C1*(1.0D0-FID) + FID*C2 BIR = AA BII = 0.0D0 RETURN 190 CONTINUE IERR=2 NZ=0 RETURN 200 CONTINUE IF(NZ.EQ.(-1)) GO TO 190 NZ=0 IERR=5 RETURN 260 CONTINUE IERR=4 NZ=0 RETURN END SUBROUTINE ZBKNU (ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, + ALIM) C***BEGIN PROLOGUE ZBKNU C***SUBSIDIARY C***PURPOSE Subsidiary to ZAIRY, ZBESH, ZBESI and ZBESK C***LIBRARY SLATEC C***TYPE ALL (CBKNU-A, ZBKNU-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE. C C***SEE ALSO ZAIRY, ZBESH, ZBESI, ZBESK C***ROUTINES CALLED D1MACH, DGAMLN, I1MACH, ZABS, ZDIV, ZEXP, ZKSCL, C ZLOG, ZMLT, ZSHCH, ZSQRT, ZUCHK C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZBKNU C DOUBLE PRECISION AA, AK, ALIM, ASCLE, A1, A2, BB, BK, BRY, CAZ, * CBI, CBR, CC, CCHI, CCHR, CKI, CKR, COEFI, COEFR, CONEI, CONER, * CRSCR, CSCLR, CSHI, CSHR, CSI, CSR, CSRR, CSSR, CTWOR, * CZEROI, CZEROR, CZI, CZR, DNU, DNU2, DPI, ELIM, ETEST, FC, FHS, * FI, FK, FKS, FMUI, FMUR, FNU, FPI, FR, G1, G2, HPI, PI, PR, PTI, * PTR, P1I, P1R, P2I, P2M, P2R, QI, QR, RAK, RCAZ, RTHPI, RZI, * RZR, R1, S, SMUI, SMUR, SPI, STI, STR, S1I, S1R, S2I, S2R, TM, * TOL, TTH, T1, T2, YI, YR, ZI, ZR, DGAMLN, D1MACH, ZABS, ELM, * CELMR, ZDR, ZDI, AS, ALAS, HELIM, CYR, CYI INTEGER I, IFLAG, INU, K, KFLAG, KK, KMAX, KODE, KODED, N, NZ, * IDUM, I1MACH, J, IC, INUB, NW DIMENSION YR(N), YI(N), CC(8), CSSR(3), CSRR(3), BRY(3), CYR(2), * CYI(2) EXTERNAL ZABS C COMPLEX Z,Y,A,B,RZ,SMU,FU,FMU,F,FLRZ,CZ,S1,S2,CSH,CCH C COMPLEX CK,P,Q,COEF,P1,P2,CBK,PT,CZERO,CONE,CTWO,ST,EZ,CS,DK C DATA KMAX / 30 / DATA CZEROR,CZEROI,CONER,CONEI,CTWOR,R1/ 1 0.0D0 , 0.0D0 , 1.0D0 , 0.0D0 , 2.0D0 , 2.0D0 / DATA DPI, RTHPI, SPI ,HPI, FPI, TTH / 1 3.14159265358979324D0, 1.25331413731550025D0, 2 1.90985931710274403D0, 1.57079632679489662D0, 3 1.89769999331517738D0, 6.66666666666666666D-01/ DATA CC(1), CC(2), CC(3), CC(4), CC(5), CC(6), CC(7), CC(8)/ 1 5.77215664901532861D-01, -4.20026350340952355D-02, 2 -4.21977345555443367D-02, 7.21894324666309954D-03, 3 -2.15241674114950973D-04, -2.01348547807882387D-05, 4 1.13302723198169588D-06, 6.11609510448141582D-09/ C***FIRST EXECUTABLE STATEMENT ZBKNU CAZ = ZABS(ZR,ZI) CSCLR = 1.0D0/TOL CRSCR = TOL CSSR(1) = CSCLR CSSR(2) = 1.0D0 CSSR(3) = CRSCR CSRR(1) = CRSCR CSRR(2) = 1.0D0 CSRR(3) = CSCLR BRY(1) = 1.0D+3*D1MACH(1)/TOL BRY(2) = 1.0D0/BRY(1) BRY(3) = D1MACH(2) NZ = 0 IFLAG = 0 KODED = KODE RCAZ = 1.0D0/CAZ STR = ZR*RCAZ STI = -ZI*RCAZ RZR = (STR+STR)*RCAZ RZI = (STI+STI)*RCAZ INU = FNU+0.5D0 DNU = FNU - INU IF (ABS(DNU).EQ.0.5D0) GO TO 110 DNU2 = 0.0D0 IF (ABS(DNU).GT.TOL) DNU2 = DNU*DNU IF (CAZ.GT.R1) GO TO 110 C----------------------------------------------------------------------- C SERIES FOR ABS(Z).LE.R1 C----------------------------------------------------------------------- FC = 1.0D0 CALL ZLOG(RZR, RZI, SMUR, SMUI, IDUM) FMUR = SMUR*DNU FMUI = SMUI*DNU CALL ZSHCH(FMUR, FMUI, CSHR, CSHI, CCHR, CCHI) IF (DNU.EQ.0.0D0) GO TO 10 FC = DNU*DPI FC = FC/SIN(FC) SMUR = CSHR/DNU SMUI = CSHI/DNU 10 CONTINUE A2 = 1.0D0 + DNU C----------------------------------------------------------------------- C GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU) C----------------------------------------------------------------------- T2 = EXP(-DGAMLN(A2,IDUM)) T1 = 1.0D0/(T2*FC) IF (ABS(DNU).GT.0.1D0) GO TO 40 C----------------------------------------------------------------------- C SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU) C----------------------------------------------------------------------- AK = 1.0D0 S = CC(1) DO 20 K=2,8 AK = AK*DNU2 TM = CC(K)*AK S = S + TM IF (ABS(TM).LT.TOL) GO TO 30 20 CONTINUE 30 G1 = -S GO TO 50 40 CONTINUE G1 = (T1-T2)/(DNU+DNU) 50 CONTINUE G2 = (T1+T2)*0.5D0 FR = FC*(CCHR*G1+SMUR*G2) FI = FC*(CCHI*G1+SMUI*G2) CALL ZEXP(FMUR, FMUI, STR, STI) PR = 0.5D0*STR/T2 PI = 0.5D0*STI/T2 CALL ZDIV(0.5D0, 0.0D0, STR, STI, PTR, PTI) QR = PTR/T1 QI = PTI/T1 S1R = FR S1I = FI S2R = PR S2I = PI AK = 1.0D0 A1 = 1.0D0 CKR = CONER CKI = CONEI BK = 1.0D0 - DNU2 IF (INU.GT.0 .OR. N.GT.1) GO TO 80 C----------------------------------------------------------------------- C GENERATE K(FNU,Z), 0.0D0 .LE. FNU .LT. 0.5D0 AND N=1 C----------------------------------------------------------------------- IF (CAZ.LT.TOL) GO TO 70 CALL ZMLT(ZR, ZI, ZR, ZI, CZR, CZI) CZR = 0.25D0*CZR CZI = 0.25D0*CZI T1 = 0.25D0*CAZ*CAZ 60 CONTINUE FR = (FR*AK+PR+QR)/BK FI = (FI*AK+PI+QI)/BK STR = 1.0D0/(AK-DNU) PR = PR*STR PI = PI*STR STR = 1.0D0/(AK+DNU) QR = QR*STR QI = QI*STR STR = CKR*CZR - CKI*CZI RAK = 1.0D0/AK CKI = (CKR*CZI+CKI*CZR)*RAK CKR = STR*RAK S1R = CKR*FR - CKI*FI + S1R S1I = CKR*FI + CKI*FR + S1I A1 = A1*T1*RAK BK = BK + AK + AK + 1.0D0 AK = AK + 1.0D0 IF (A1.GT.TOL) GO TO 60 70 CONTINUE YR(1) = S1R YI(1) = S1I IF (KODED.EQ.1) RETURN CALL ZEXP(ZR, ZI, STR, STI) CALL ZMLT(S1R, S1I, STR, STI, YR(1), YI(1)) RETURN C----------------------------------------------------------------------- C GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE C----------------------------------------------------------------------- 80 CONTINUE IF (CAZ.LT.TOL) GO TO 100 CALL ZMLT(ZR, ZI, ZR, ZI, CZR, CZI) CZR = 0.25D0*CZR CZI = 0.25D0*CZI T1 = 0.25D0*CAZ*CAZ 90 CONTINUE FR = (FR*AK+PR+QR)/BK FI = (FI*AK+PI+QI)/BK STR = 1.0D0/(AK-DNU) PR = PR*STR PI = PI*STR STR = 1.0D0/(AK+DNU) QR = QR*STR QI = QI*STR STR = CKR*CZR - CKI*CZI RAK = 1.0D0/AK CKI = (CKR*CZI+CKI*CZR)*RAK CKR = STR*RAK S1R = CKR*FR - CKI*FI + S1R S1I = CKR*FI + CKI*FR + S1I STR = PR - FR*AK STI = PI - FI*AK S2R = CKR*STR - CKI*STI + S2R S2I = CKR*STI + CKI*STR + S2I A1 = A1*T1*RAK BK = BK + AK + AK + 1.0D0 AK = AK + 1.0D0 IF (A1.GT.TOL) GO TO 90 100 CONTINUE KFLAG = 2 A1 = FNU + 1.0D0 AK = A1*ABS(SMUR) IF (AK.GT.ALIM) KFLAG = 3 STR = CSSR(KFLAG) P2R = S2R*STR P2I = S2I*STR CALL ZMLT(P2R, P2I, RZR, RZI, S2R, S2I) S1R = S1R*STR S1I = S1I*STR IF (KODED.EQ.1) GO TO 210 CALL ZEXP(ZR, ZI, FR, FI) CALL ZMLT(S1R, S1I, FR, FI, S1R, S1I) CALL ZMLT(S2R, S2I, FR, FI, S2R, S2I) GO TO 210 C----------------------------------------------------------------------- C IFLAG=0 MEANS NO UNDERFLOW OCCURRED C IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH C KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD C RECURSION C----------------------------------------------------------------------- 110 CONTINUE CALL ZSQRT(ZR, ZI, STR, STI) CALL ZDIV(RTHPI, CZEROI, STR, STI, COEFR, COEFI) KFLAG = 2 IF (KODED.EQ.2) GO TO 120 IF (ZR.GT.ALIM) GO TO 290 C BLANK LINE STR = EXP(-ZR)*CSSR(KFLAG) STI = -STR*SIN(ZI) STR = STR*COS(ZI) CALL ZMLT(COEFR, COEFI, STR, STI, COEFR, COEFI) 120 CONTINUE IF (ABS(DNU).EQ.0.5D0) GO TO 300 C----------------------------------------------------------------------- C MILLER ALGORITHM FOR ABS(Z).GT.R1 C----------------------------------------------------------------------- AK = COS(DPI*DNU) AK = ABS(AK) IF (AK.EQ.CZEROR) GO TO 300 FHS = ABS(0.25D0-DNU2) IF (FHS.EQ.CZEROR) GO TO 300 C----------------------------------------------------------------------- C COMPUTE R2=F(E). IF ABS(Z).GE.R2, USE FORWARD RECURRENCE TO C DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON C 12.LE.E.LE.60. E IS COMPUTED FROM 2**(-E)=B**(1-I1MACH(14))= C TOL WHERE B IS THE BASE OF THE ARITHMETIC. C----------------------------------------------------------------------- T1 = I1MACH(14)-1 T1 = T1*D1MACH(5)*3.321928094D0 T1 = MAX(T1,12.0D0) T1 = MIN(T1,60.0D0) T2 = TTH*T1 - 6.0D0 IF (ZR.NE.0.0D0) GO TO 130 T1 = HPI GO TO 140 130 CONTINUE T1 = DATAN(ZI/ZR) T1 = ABS(T1) 140 CONTINUE IF (T2.GT.CAZ) GO TO 170 C----------------------------------------------------------------------- C FORWARD RECURRENCE LOOP WHEN ABS(Z).GE.R2 C----------------------------------------------------------------------- ETEST = AK/(DPI*CAZ*TOL) FK = CONER IF (ETEST.LT.CONER) GO TO 180 FKS = CTWOR CKR = CAZ + CAZ + CTWOR P1R = CZEROR P2R = CONER DO 150 I=1,KMAX AK = FHS/FKS CBR = CKR/(FK+CONER) PTR = P2R P2R = CBR*P2R - P1R*AK P1R = PTR CKR = CKR + CTWOR FKS = FKS + FK + FK + CTWOR FHS = FHS + FK + FK FK = FK + CONER STR = ABS(P2R)*FK IF (ETEST.LT.STR) GO TO 160 150 CONTINUE GO TO 310 160 CONTINUE FK = FK + SPI*T1*SQRT(T2/CAZ) FHS = ABS(0.25D0-DNU2) GO TO 180 170 CONTINUE C----------------------------------------------------------------------- C COMPUTE BACKWARD INDEX K FOR ABS(Z).LT.R2 C----------------------------------------------------------------------- A2 = SQRT(CAZ) AK = FPI*AK/(TOL*SQRT(A2)) AA = 3.0D0*T1/(1.0D0+CAZ) BB = 14.7D0*T1/(28.0D0+CAZ) AK = (LOG(AK)+CAZ*COS(AA)/(1.0D0+0.008D0*CAZ))/COS(BB) FK = 0.12125D0*AK*AK/CAZ + 1.5D0 180 CONTINUE C----------------------------------------------------------------------- C BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM C----------------------------------------------------------------------- K = FK FK = K FKS = FK*FK P1R = CZEROR P1I = CZEROI P2R = TOL P2I = CZEROI CSR = P2R CSI = P2I DO 190 I=1,K A1 = FKS - FK AK = (FKS+FK)/(A1+FHS) RAK = 2.0D0/(FK+CONER) CBR = (FK+ZR)*RAK CBI = ZI*RAK PTR = P2R PTI = P2I P2R = (PTR*CBR-PTI*CBI-P1R)*AK P2I = (PTI*CBR+PTR*CBI-P1I)*AK P1R = PTR P1I = PTI CSR = CSR + P2R CSI = CSI + P2I FKS = A1 - FK + CONER FK = FK - CONER 190 CONTINUE C----------------------------------------------------------------------- C COMPUTE (P2/CS)=(P2/ABS(CS))*(CONJG(CS)/ABS(CS)) FOR BETTER C SCALING C----------------------------------------------------------------------- TM = ZABS(CSR,CSI) PTR = 1.0D0/TM S1R = P2R*PTR S1I = P2I*PTR CSR = CSR*PTR CSI = -CSI*PTR CALL ZMLT(COEFR, COEFI, S1R, S1I, STR, STI) CALL ZMLT(STR, STI, CSR, CSI, S1R, S1I) IF (INU.GT.0 .OR. N.GT.1) GO TO 200 ZDR = ZR ZDI = ZI IF(IFLAG.EQ.1) GO TO 270 GO TO 240 200 CONTINUE C----------------------------------------------------------------------- C COMPUTE P1/P2=(P1/ABS(P2)*CONJG(P2)/ABS(P2) FOR SCALING C----------------------------------------------------------------------- TM = ZABS(P2R,P2I) PTR = 1.0D0/TM P1R = P1R*PTR P1I = P1I*PTR P2R = P2R*PTR P2I = -P2I*PTR CALL ZMLT(P1R, P1I, P2R, P2I, PTR, PTI) STR = DNU + 0.5D0 - PTR STI = -PTI CALL ZDIV(STR, STI, ZR, ZI, STR, STI) STR = STR + 1.0D0 CALL ZMLT(STR, STI, S1R, S1I, S2R, S2I) C----------------------------------------------------------------------- C FORWARD RECURSION ON THE THREE TERM RECURSION WITH RELATION WITH C SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3 C----------------------------------------------------------------------- 210 CONTINUE STR = DNU + 1.0D0 CKR = STR*RZR CKI = STR*RZI IF (N.EQ.1) INU = INU - 1 IF (INU.GT.0) GO TO 220 IF (N.GT.1) GO TO 215 S1R = S2R S1I = S2I 215 CONTINUE ZDR = ZR ZDI = ZI IF(IFLAG.EQ.1) GO TO 270 GO TO 240 220 CONTINUE INUB = 1 IF(IFLAG.EQ.1) GO TO 261 225 CONTINUE P1R = CSRR(KFLAG) ASCLE = BRY(KFLAG) DO 230 I=INUB,INU STR = S2R STI = S2I S2R = CKR*STR - CKI*STI + S1R S2I = CKR*STI + CKI*STR + S1I S1R = STR S1I = STI CKR = CKR + RZR CKI = CKI + RZI IF (KFLAG.GE.3) GO TO 230 P2R = S2R*P1R P2I = S2I*P1R STR = ABS(P2R) STI = ABS(P2I) P2M = MAX(STR,STI) IF (P2M.LE.ASCLE) GO TO 230 KFLAG = KFLAG + 1 ASCLE = BRY(KFLAG) S1R = S1R*P1R S1I = S1I*P1R S2R = P2R S2I = P2I STR = CSSR(KFLAG) S1R = S1R*STR S1I = S1I*STR S2R = S2R*STR S2I = S2I*STR P1R = CSRR(KFLAG) 230 CONTINUE IF (N.NE.1) GO TO 240 S1R = S2R S1I = S2I 240 CONTINUE STR = CSRR(KFLAG) YR(1) = S1R*STR YI(1) = S1I*STR IF (N.EQ.1) RETURN YR(2) = S2R*STR YI(2) = S2I*STR IF (N.EQ.2) RETURN KK = 2 250 CONTINUE KK = KK + 1 IF (KK.GT.N) RETURN P1R = CSRR(KFLAG) ASCLE = BRY(KFLAG) DO 260 I=KK,N P2R = S2R P2I = S2I S2R = CKR*P2R - CKI*P2I + S1R S2I = CKI*P2R + CKR*P2I + S1I S1R = P2R S1I = P2I CKR = CKR + RZR CKI = CKI + RZI P2R = S2R*P1R P2I = S2I*P1R YR(I) = P2R YI(I) = P2I IF (KFLAG.GE.3) GO TO 260 STR = ABS(P2R) STI = ABS(P2I) P2M = MAX(STR,STI) IF (P2M.LE.ASCLE) GO TO 260 KFLAG = KFLAG + 1 ASCLE = BRY(KFLAG) S1R = S1R*P1R S1I = S1I*P1R S2R = P2R S2I = P2I STR = CSSR(KFLAG) S1R = S1R*STR S1I = S1I*STR S2R = S2R*STR S2I = S2I*STR P1R = CSRR(KFLAG) 260 CONTINUE RETURN C----------------------------------------------------------------------- C IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW C----------------------------------------------------------------------- 261 CONTINUE HELIM = 0.5D0*ELIM ELM = EXP(-ELIM) CELMR = ELM ASCLE = BRY(1) ZDR = ZR ZDI = ZI IC = -1 J = 2 DO 262 I=1,INU STR = S2R STI = S2I S2R = STR*CKR-STI*CKI+S1R S2I = STI*CKR+STR*CKI+S1I S1R = STR S1I = STI CKR = CKR+RZR CKI = CKI+RZI AS = ZABS(S2R,S2I) ALAS = LOG(AS) P2R = -ZDR+ALAS IF(P2R.LT.(-ELIM)) GO TO 263 CALL ZLOG(S2R,S2I,STR,STI,IDUM) P2R = -ZDR+STR P2I = -ZDI+STI P2M = EXP(P2R)/TOL P1R = P2M*COS(P2I) P1I = P2M*SIN(P2I) CALL ZUCHK(P1R,P1I,NW,ASCLE,TOL) IF(NW.NE.0) GO TO 263 J = 3 - J CYR(J) = P1R CYI(J) = P1I IF(IC.EQ.(I-1)) GO TO 264 IC = I GO TO 262 263 CONTINUE IF(ALAS.LT.HELIM) GO TO 262 ZDR = ZDR-ELIM S1R = S1R*CELMR S1I = S1I*CELMR S2R = S2R*CELMR S2I = S2I*CELMR 262 CONTINUE IF(N.NE.1) GO TO 270 S1R = S2R S1I = S2I GO TO 270 264 CONTINUE KFLAG = 1 INUB = I+1 S2R = CYR(J) S2I = CYI(J) J = 3 - J S1R = CYR(J) S1I = CYI(J) IF(INUB.LE.INU) GO TO 225 IF(N.NE.1) GO TO 240 S1R = S2R S1I = S2I GO TO 240 270 CONTINUE YR(1) = S1R YI(1) = S1I IF(N.EQ.1) GO TO 280 YR(2) = S2R YI(2) = S2I 280 CONTINUE ASCLE = BRY(1) CALL ZKSCL(ZDR,ZDI,FNU,N,YR,YI,NZ,RZR,RZI,ASCLE,TOL,ELIM) INU = N - NZ IF (INU.LE.0) RETURN KK = NZ + 1 S1R = YR(KK) S1I = YI(KK) YR(KK) = S1R*CSRR(1) YI(KK) = S1I*CSRR(1) IF (INU.EQ.1) RETURN KK = NZ + 2 S2R = YR(KK) S2I = YI(KK) YR(KK) = S2R*CSRR(1) YI(KK) = S2I*CSRR(1) IF (INU.EQ.2) RETURN T2 = FNU + (KK-1) CKR = T2*RZR CKI = T2*RZI KFLAG = 1 GO TO 250 290 CONTINUE C----------------------------------------------------------------------- C SCALE BY EXP(Z), IFLAG = 1 CASES C----------------------------------------------------------------------- KODED = 2 IFLAG = 1 KFLAG = 2 GO TO 120 C----------------------------------------------------------------------- C FNU=HALF ODD INTEGER CASE, DNU=-0.5 C----------------------------------------------------------------------- 300 CONTINUE S1R = COEFR S1I = COEFI S2R = COEFR S2I = COEFI GO TO 210 C C 310 CONTINUE NZ=-2 RETURN END SUBROUTINE ZBUNI (ZR, ZI, FNU, KODE, N, YR, YI, NZ, NUI, NLAST, + FNUL, TOL, ELIM, ALIM) C***BEGIN PROLOGUE ZBUNI C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESI and ZBESK C***LIBRARY SLATEC C***TYPE ALL (CBUNI-A, ZBUNI-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZBUNI COMPUTES THE I BESSEL FUNCTION FOR LARGE ABS(Z).GT. C FNUL AND FNU+N-1.LT.FNUL. THE ORDER IS INCREASED FROM C FNU+N-1 GREATER THAN FNUL BY ADDING NUI AND COMPUTING C ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR I(FNU,Z) C ON IFORM=1 AND THE EXPANSION FOR J(FNU,Z) ON IFORM=2 C C***SEE ALSO ZBESI, ZBESK C***ROUTINES CALLED D1MACH, ZABS, ZUNI1, ZUNI2 C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZBUNI C COMPLEX CSCL,CSCR,CY,RZ,ST,S1,S2,Y,Z DOUBLE PRECISION ALIM, AX, AY, CSCLR, CSCRR, CYI, CYR, DFNU, * ELIM, FNU, FNUI, FNUL, GNU, RAZ, RZI, RZR, STI, STR, S1I, S1R, * S2I, S2R, TOL, YI, YR, ZI, ZR, ZABS, ASCLE, BRY, C1R, C1I, C1M, * D1MACH INTEGER I, IFLAG, IFORM, K, KODE, N, NL, NLAST, NUI, NW, NZ DIMENSION YR(N), YI(N), CYR(2), CYI(2), BRY(3) EXTERNAL ZABS C***FIRST EXECUTABLE STATEMENT ZBUNI NZ = 0 AX = ABS(ZR)*1.7321D0 AY = ABS(ZI) IFORM = 1 IF (AY.GT.AX) IFORM = 2 IF (NUI.EQ.0) GO TO 60 FNUI = NUI DFNU = FNU + (N-1) GNU = DFNU + FNUI IF (IFORM.EQ.2) GO TO 10 C----------------------------------------------------------------------- C ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN C -PI/3.LE.ARG(Z).LE.PI/3 C----------------------------------------------------------------------- CALL ZUNI1(ZR, ZI, GNU, KODE, 2, CYR, CYI, NW, NLAST, FNUL, TOL, * ELIM, ALIM) GO TO 20 10 CONTINUE C----------------------------------------------------------------------- C ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU C APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I C AND HPI=PI/2 C----------------------------------------------------------------------- CALL ZUNI2(ZR, ZI, GNU, KODE, 2, CYR, CYI, NW, NLAST, FNUL, TOL, * ELIM, ALIM) 20 CONTINUE IF (NW.LT.0) GO TO 50 IF (NW.NE.0) GO TO 90 STR = ZABS(CYR(1),CYI(1)) C---------------------------------------------------------------------- C SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED C---------------------------------------------------------------------- BRY(1)=1.0D+3*D1MACH(1)/TOL BRY(2) = 1.0D0/BRY(1) BRY(3) = BRY(2) IFLAG = 2 ASCLE = BRY(2) CSCLR = 1.0D0 IF (STR.GT.BRY(1)) GO TO 21 IFLAG = 1 ASCLE = BRY(1) CSCLR = 1.0D0/TOL GO TO 25 21 CONTINUE IF (STR.LT.BRY(2)) GO TO 25 IFLAG = 3 ASCLE=BRY(3) CSCLR = TOL 25 CONTINUE CSCRR = 1.0D0/CSCLR S1R = CYR(2)*CSCLR S1I = CYI(2)*CSCLR S2R = CYR(1)*CSCLR S2I = CYI(1)*CSCLR RAZ = 1.0D0/ZABS(ZR,ZI) STR = ZR*RAZ STI = -ZI*RAZ RZR = (STR+STR)*RAZ RZI = (STI+STI)*RAZ DO 30 I=1,NUI STR = S2R STI = S2I S2R = (DFNU+FNUI)*(RZR*STR-RZI*STI) + S1R S2I = (DFNU+FNUI)*(RZR*STI+RZI*STR) + S1I S1R = STR S1I = STI FNUI = FNUI - 1.0D0 IF (IFLAG.GE.3) GO TO 30 STR = S2R*CSCRR STI = S2I*CSCRR C1R = ABS(STR) C1I = ABS(STI) C1M = MAX(C1R,C1I) IF (C1M.LE.ASCLE) GO TO 30 IFLAG = IFLAG+1 ASCLE = BRY(IFLAG) S1R = S1R*CSCRR S1I = S1I*CSCRR S2R = STR S2I = STI CSCLR = CSCLR*TOL CSCRR = 1.0D0/CSCLR S1R = S1R*CSCLR S1I = S1I*CSCLR S2R = S2R*CSCLR S2I = S2I*CSCLR 30 CONTINUE YR(N) = S2R*CSCRR YI(N) = S2I*CSCRR IF (N.EQ.1) RETURN NL = N - 1 FNUI = NL K = NL DO 40 I=1,NL STR = S2R STI = S2I S2R = (FNU+FNUI)*(RZR*STR-RZI*STI) + S1R S2I = (FNU+FNUI)*(RZR*STI+RZI*STR) + S1I S1R = STR S1I = STI STR = S2R*CSCRR STI = S2I*CSCRR YR(K) = STR YI(K) = STI FNUI = FNUI - 1.0D0 K = K - 1 IF (IFLAG.GE.3) GO TO 40 C1R = ABS(STR) C1I = ABS(STI) C1M = MAX(C1R,C1I) IF (C1M.LE.ASCLE) GO TO 40 IFLAG = IFLAG+1 ASCLE = BRY(IFLAG) S1R = S1R*CSCRR S1I = S1I*CSCRR S2R = STR S2I = STI CSCLR = CSCLR*TOL CSCRR = 1.0D0/CSCLR S1R = S1R*CSCLR S1I = S1I*CSCLR S2R = S2R*CSCLR S2I = S2I*CSCLR 40 CONTINUE RETURN 50 CONTINUE NZ = -1 IF(NW.EQ.(-2)) NZ=-2 RETURN 60 CONTINUE IF (IFORM.EQ.2) GO TO 70 C----------------------------------------------------------------------- C ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN C -PI/3.LE.ARG(Z).LE.PI/3 C----------------------------------------------------------------------- CALL ZUNI1(ZR, ZI, FNU, KODE, N, YR, YI, NW, NLAST, FNUL, TOL, * ELIM, ALIM) GO TO 80 70 CONTINUE C----------------------------------------------------------------------- C ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU C APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I C AND HPI=PI/2 C----------------------------------------------------------------------- CALL ZUNI2(ZR, ZI, FNU, KODE, N, YR, YI, NW, NLAST, FNUL, TOL, * ELIM, ALIM) 80 CONTINUE IF (NW.LT.0) GO TO 50 NZ = NW RETURN 90 CONTINUE NLAST = N RETURN END SUBROUTINE ZDIV (AR, AI, BR, BI, CR, CI) C***BEGIN PROLOGUE ZDIV C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESH, ZBESI, ZBESJ, ZBESK, ZBESY, ZAIRY and C ZBIRY C***LIBRARY SLATEC C***TYPE ALL (ZDIV-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C DOUBLE PRECISION COMPLEX DIVIDE C=A/B. C C***SEE ALSO ZAIRY, ZBESH, ZBESI, ZBESJ, ZBESK, ZBESY, ZBIRY C***ROUTINES CALLED ZABS C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZDIV DOUBLE PRECISION AR, AI, BR, BI, CR, CI, BM, CA, CB, CC, CD DOUBLE PRECISION ZABS EXTERNAL ZABS C***FIRST EXECUTABLE STATEMENT ZDIV BM = 1.0D0/ZABS(BR,BI) CC = BR*BM CD = BI*BM CA = (AR*CC+AI*CD)*BM CB = (AI*CC-AR*CD)*BM CR = CA CI = CB RETURN END SUBROUTINE ZEXP (AR, AI, BR, BI) C***BEGIN PROLOGUE ZEXP C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESH, ZBESI, ZBESJ, ZBESK, ZBESY, ZAIRY and C ZBIRY C***LIBRARY SLATEC C***TYPE ALL (ZEXP-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C DOUBLE PRECISION COMPLEX EXPONENTIAL FUNCTION B=EXP(A) C C***SEE ALSO ZAIRY, ZBESH, ZBESI, ZBESJ, ZBESK, ZBESY, ZBIRY C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZEXP DOUBLE PRECISION AR, AI, BR, BI, ZM, CA, CB C***FIRST EXECUTABLE STATEMENT ZEXP ZM = EXP(AR) CA = ZM*COS(AI) CB = ZM*SIN(AI) BR = CA BI = CB RETURN END SUBROUTINE ZKSCL (ZRR, ZRI, FNU, N, YR, YI, NZ, RZR, RZI, ASCLE, + TOL, ELIM) C***BEGIN PROLOGUE ZKSCL C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESK C***LIBRARY SLATEC C***TYPE ALL (CKSCL-A, ZKSCL-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE C ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN C RETURN WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL. C C***SEE ALSO ZBESK C***ROUTINES CALLED ZABS, ZLOG, ZUCHK C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZKSCL C COMPLEX CK,CS,CY,CZERO,RZ,S1,S2,Y,ZR,ZD,CELM DOUBLE PRECISION ACS, AS, ASCLE, CKI, CKR, CSI, CSR, CYI, * CYR, ELIM, FN, FNU, RZI, RZR, STR, S1I, S1R, S2I, * S2R, TOL, YI, YR, ZEROI, ZEROR, ZRI, ZRR, ZABS, * ZDR, ZDI, CELMR, ELM, HELIM, ALAS INTEGER I, IC, IDUM, KK, N, NN, NW, NZ DIMENSION YR(N), YI(N), CYR(2), CYI(2) EXTERNAL ZABS DATA ZEROR,ZEROI / 0.0D0 , 0.0D0 / C***FIRST EXECUTABLE STATEMENT ZKSCL NZ = 0 IC = 0 NN = MIN(2,N) DO 10 I=1,NN S1R = YR(I) S1I = YI(I) CYR(I) = S1R CYI(I) = S1I AS = ZABS(S1R,S1I) ACS = -ZRR + LOG(AS) NZ = NZ + 1 YR(I) = ZEROR YI(I) = ZEROI IF (ACS.LT.(-ELIM)) GO TO 10 CALL ZLOG(S1R, S1I, CSR, CSI, IDUM) CSR = CSR - ZRR CSI = CSI - ZRI STR = EXP(CSR)/TOL CSR = STR*COS(CSI) CSI = STR*SIN(CSI) CALL ZUCHK(CSR, CSI, NW, ASCLE, TOL) IF (NW.NE.0) GO TO 10 YR(I) = CSR YI(I) = CSI IC = I NZ = NZ - 1 10 CONTINUE IF (N.EQ.1) RETURN IF (IC.GT.1) GO TO 20 YR(1) = ZEROR YI(1) = ZEROI NZ = 2 20 CONTINUE IF (N.EQ.2) RETURN IF (NZ.EQ.0) RETURN FN = FNU + 1.0D0 CKR = FN*RZR CKI = FN*RZI S1R = CYR(1) S1I = CYI(1) S2R = CYR(2) S2I = CYI(2) HELIM = 0.5D0*ELIM ELM = EXP(-ELIM) CELMR = ELM ZDR = ZRR ZDI = ZRI C C FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF C S2 GETS LARGER THAN EXP(ELIM/2) C DO 30 I=3,N KK = I CSR = S2R CSI = S2I S2R = CKR*CSR - CKI*CSI + S1R S2I = CKI*CSR + CKR*CSI + S1I S1R = CSR S1I = CSI CKR = CKR + RZR CKI = CKI + RZI AS = ZABS(S2R,S2I) ALAS = LOG(AS) ACS = -ZDR + ALAS NZ = NZ + 1 YR(I) = ZEROR YI(I) = ZEROI IF (ACS.LT.(-ELIM)) GO TO 25 CALL ZLOG(S2R, S2I, CSR, CSI, IDUM) CSR = CSR - ZDR CSI = CSI - ZDI STR = EXP(CSR)/TOL CSR = STR*COS(CSI) CSI = STR*SIN(CSI) CALL ZUCHK(CSR, CSI, NW, ASCLE, TOL) IF (NW.NE.0) GO TO 25 YR(I) = CSR YI(I) = CSI NZ = NZ - 1 IF (IC.EQ.KK-1) GO TO 40 IC = KK GO TO 30 25 CONTINUE IF(ALAS.LT.HELIM) GO TO 30 ZDR = ZDR - ELIM S1R = S1R*CELMR S1I = S1I*CELMR S2R = S2R*CELMR S2I = S2I*CELMR 30 CONTINUE NZ = N IF(IC.EQ.N) NZ=N-1 GO TO 45 40 CONTINUE NZ = KK - 2 45 CONTINUE DO 50 I=1,NZ YR(I) = ZEROR YI(I) = ZEROI 50 CONTINUE RETURN END SUBROUTINE ZLOG (AR, AI, BR, BI, IERR) C***BEGIN PROLOGUE ZLOG C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESH, ZBESI, ZBESJ, ZBESK, ZBESY, ZAIRY and C ZBIRY C***LIBRARY SLATEC C***TYPE ALL (ZLOG-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C DOUBLE PRECISION COMPLEX LOGARITHM B=CLOG(A) C IERR=0,NORMAL RETURN IERR=1, Z=CMPLX(0.0,0.0) C***SEE ALSO ZAIRY, ZBESH, ZBESI, ZBESJ, ZBESK, ZBESY, ZBIRY C***ROUTINES CALLED ZABS C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZLOG DOUBLE PRECISION AR, AI, BR, BI, ZM, DTHETA, DPI, DHPI DOUBLE PRECISION ZABS INTEGER IERR EXTERNAL ZABS DATA DPI , DHPI / 3.141592653589793238462643383D+0, 1 1.570796326794896619231321696D+0/ C***FIRST EXECUTABLE STATEMENT ZLOG IERR=0 IF (AR.EQ.0.0D+0) GO TO 10 IF (AI.EQ.0.0D+0) GO TO 20 DTHETA = DATAN(AI/AR) IF (DTHETA.LE.0.0D+0) GO TO 40 IF (AR.LT.0.0D+0) DTHETA = DTHETA - DPI GO TO 50 10 IF (AI.EQ.0.0D+0) GO TO 60 BI = DHPI BR = LOG(ABS(AI)) IF (AI.LT.0.0D+0) BI = -BI RETURN 20 IF (AR.GT.0.0D+0) GO TO 30 BR = LOG(ABS(AR)) BI = DPI RETURN 30 BR = LOG(AR) BI = 0.0D+0 RETURN 40 IF (AR.LT.0.0D+0) DTHETA = DTHETA + DPI 50 ZM = ZABS(AR,AI) BR = LOG(ZM) BI = DTHETA RETURN 60 CONTINUE IERR=1 RETURN END SUBROUTINE ZMLRI (ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL) C***BEGIN PROLOGUE ZMLRI C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESI and ZBESK C***LIBRARY SLATEC C***TYPE ALL (CMLRI-A, ZMLRI-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE C MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES. C C***SEE ALSO ZBESI, ZBESK C***ROUTINES CALLED D1MACH, DGAMLN, ZABS, ZEXP, ZLOG, ZMLT C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZMLRI C COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z DOUBLE PRECISION ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI, * CNORMR, CONEI, CONER, FKAP, FKK, FLAM, FNF, FNU, PTI, PTR, P1I, * P1R, P2I, P2R, RAZ, RHO, RHO2, RZI, RZR, SCLE, STI, STR, SUMI, * SUMR, TFNF, TOL, TST, YI, YR, ZEROI, ZEROR, ZI, ZR, DGAMLN, * D1MACH, ZABS INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N, NZ DIMENSION YR(N), YI(N) EXTERNAL ZABS DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 / C***FIRST EXECUTABLE STATEMENT ZMLRI SCLE = D1MACH(1)/TOL NZ=0 AZ = ZABS(ZR,ZI) IAZ = AZ IFNU = FNU INU = IFNU + N - 1 AT = IAZ + 1.0D0 RAZ = 1.0D0/AZ STR = ZR*RAZ STI = -ZI*RAZ CKR = STR*AT*RAZ CKI = STI*AT*RAZ RZR = (STR+STR)*RAZ RZI = (STI+STI)*RAZ P1R = ZEROR P1I = ZEROI P2R = CONER P2I = CONEI ACK = (AT+1.0D0)*RAZ RHO = ACK + SQRT(ACK*ACK-1.0D0) RHO2 = RHO*RHO TST = (RHO2+RHO2)/((RHO2-1.0D0)*(RHO-1.0D0)) TST = TST/TOL C----------------------------------------------------------------------- C COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES C----------------------------------------------------------------------- AK = AT DO 10 I=1,80 PTR = P2R PTI = P2I P2R = P1R - (CKR*PTR-CKI*PTI) P2I = P1I - (CKI*PTR+CKR*PTI) P1R = PTR P1I = PTI CKR = CKR + RZR CKI = CKI + RZI AP = ZABS(P2R,P2I) IF (AP.GT.TST*AK*AK) GO TO 20 AK = AK + 1.0D0 10 CONTINUE GO TO 110 20 CONTINUE I = I + 1 K = 0 IF (INU.LT.IAZ) GO TO 40 C----------------------------------------------------------------------- C COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS C----------------------------------------------------------------------- P1R = ZEROR P1I = ZEROI P2R = CONER P2I = CONEI AT = INU + 1.0D0 STR = ZR*RAZ STI = -ZI*RAZ CKR = STR*AT*RAZ CKI = STI*AT*RAZ ACK = AT*RAZ TST = SQRT(ACK/TOL) ITIME = 1 DO 30 K=1,80 PTR = P2R PTI = P2I P2R = P1R - (CKR*PTR-CKI*PTI) P2I = P1I - (CKR*PTI+CKI*PTR) P1R = PTR P1I = PTI CKR = CKR + RZR CKI = CKI + RZI AP = ZABS(P2R,P2I) IF (AP.LT.TST) GO TO 30 IF (ITIME.EQ.2) GO TO 40 ACK = ZABS(CKR,CKI) FLAM = ACK + SQRT(ACK*ACK-1.0D0) FKAP = AP/ZABS(P1R,P1I) RHO = MIN(FLAM,FKAP) TST = TST*SQRT(RHO/(RHO*RHO-1.0D0)) ITIME = 2 30 CONTINUE GO TO 110 40 CONTINUE C----------------------------------------------------------------------- C BACKWARD RECURRENCE AND SUM NORMALIZING RELATION C----------------------------------------------------------------------- K = K + 1 KK = MAX(I+IAZ,K+INU) FKK = KK P1R = ZEROR P1I = ZEROI C----------------------------------------------------------------------- C SCALE P2 AND SUM BY SCLE C----------------------------------------------------------------------- P2R = SCLE P2I = ZEROI FNF = FNU - IFNU TFNF = FNF + FNF BK = DGAMLN(FKK+TFNF+1.0D0,IDUM) - DGAMLN(FKK+1.0D0,IDUM) - * DGAMLN(TFNF+1.0D0,IDUM) BK = EXP(BK) SUMR = ZEROR SUMI = ZEROI KM = KK - INU DO 50 I=1,KM PTR = P2R PTI = P2I P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI) P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI) P1R = PTR P1I = PTI AK = 1.0D0 - TFNF/(FKK+TFNF) ACK = BK*AK SUMR = SUMR + (ACK+BK)*P1R SUMI = SUMI + (ACK+BK)*P1I BK = ACK FKK = FKK - 1.0D0 50 CONTINUE YR(N) = P2R YI(N) = P2I IF (N.EQ.1) GO TO 70 DO 60 I=2,N PTR = P2R PTI = P2I P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI) P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI) P1R = PTR P1I = PTI AK = 1.0D0 - TFNF/(FKK+TFNF) ACK = BK*AK SUMR = SUMR + (ACK+BK)*P1R SUMI = SUMI + (ACK+BK)*P1I BK = ACK FKK = FKK - 1.0D0 M = N - I + 1 YR(M) = P2R YI(M) = P2I 60 CONTINUE 70 CONTINUE IF (IFNU.LE.0) GO TO 90 DO 80 I=1,IFNU PTR = P2R PTI = P2I P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI) P2I = P1I + (FKK+FNF)*(RZR*PTI+RZI*PTR) P1R = PTR P1I = PTI AK = 1.0D0 - TFNF/(FKK+TFNF) ACK = BK*AK SUMR = SUMR + (ACK+BK)*P1R SUMI = SUMI + (ACK+BK)*P1I BK = ACK FKK = FKK - 1.0D0 80 CONTINUE 90 CONTINUE PTR = ZR PTI = ZI IF (KODE.EQ.2) PTR = ZEROR CALL ZLOG(RZR, RZI, STR, STI, IDUM) P1R = -FNF*STR + PTR P1I = -FNF*STI + PTI AP = DGAMLN(1.0D0+FNF,IDUM) PTR = P1R - AP PTI = P1I C----------------------------------------------------------------------- C THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW C IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES C----------------------------------------------------------------------- P2R = P2R + SUMR P2I = P2I + SUMI AP = ZABS(P2R,P2I) P1R = 1.0D0/AP CALL ZEXP(PTR, PTI, STR, STI) CKR = STR*P1R CKI = STI*P1R PTR = P2R*P1R PTI = -P2I*P1R CALL ZMLT(CKR, CKI, PTR, PTI, CNORMR, CNORMI) DO 100 I=1,N STR = YR(I)*CNORMR - YI(I)*CNORMI YI(I) = YR(I)*CNORMI + YI(I)*CNORMR YR(I) = STR 100 CONTINUE RETURN 110 CONTINUE NZ=-2 RETURN END SUBROUTINE ZMLT (AR, AI, BR, BI, CR, CI) C***BEGIN PROLOGUE ZMLT C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESH, ZBESI, ZBESJ, ZBESK, ZBESY, ZAIRY and C ZBIRY C***LIBRARY SLATEC C***TYPE ALL (ZMLT-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C DOUBLE PRECISION COMPLEX MULTIPLY, C=A*B. C C***SEE ALSO ZAIRY, ZBESH, ZBESI, ZBESJ, ZBESK, ZBESY, ZBIRY C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZMLT DOUBLE PRECISION AR, AI, BR, BI, CR, CI, CA, CB C***FIRST EXECUTABLE STATEMENT ZMLT CA = AR*BR - AI*BI CB = AR*BI + AI*BR CR = CA CI = CB RETURN END SUBROUTINE ZRATI (ZR, ZI, FNU, N, CYR, CYI, TOL) C***BEGIN PROLOGUE ZRATI C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESH, ZBESI and ZBESK C***LIBRARY SLATEC C***TYPE ALL (CRATI-A, ZRATI-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD C RECURRENCE. THE STARTING INDEX IS DETERMINED BY FORWARD C RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B, C MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973, C BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER, C BY D. J. SOOKNE. C C***SEE ALSO ZBESH, ZBESI, ZBESK C***ROUTINES CALLED ZABS, ZDIV C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZRATI DOUBLE PRECISION AK, AMAGZ, AP1, AP2, ARG, AZ, CDFNUI, CDFNUR, * CONEI, CONER, CYI, CYR, CZEROI, CZEROR, DFNU, FDNU, FLAM, FNU, * FNUP, PTI, PTR, P1I, P1R, P2I, P2R, RAK, RAP1, RHO, RT2, RZI, * RZR, TEST, TEST1, TOL, TTI, TTR, T1I, T1R, ZI, ZR, ZABS INTEGER I, ID, IDNU, INU, ITIME, K, KK, MAGZ, N DIMENSION CYR(N), CYI(N) EXTERNAL ZABS DATA CZEROR,CZEROI,CONER,CONEI,RT2/ 1 0.0D0, 0.0D0, 1.0D0, 0.0D0, 1.41421356237309505D0 / C***FIRST EXECUTABLE STATEMENT ZRATI AZ = ZABS(ZR,ZI) INU = FNU IDNU = INU + N - 1 MAGZ = AZ AMAGZ = MAGZ+1 FDNU = IDNU FNUP = MAX(AMAGZ,FDNU) ID = IDNU - MAGZ - 1 ITIME = 1 K = 1 PTR = 1.0D0/AZ RZR = PTR*(ZR+ZR)*PTR RZI = -PTR*(ZI+ZI)*PTR T1R = RZR*FNUP T1I = RZI*FNUP P2R = -T1R P2I = -T1I P1R = CONER P1I = CONEI T1R = T1R + RZR T1I = T1I + RZI IF (ID.GT.0) ID = 0 AP2 = ZABS(P2R,P2I) AP1 = ZABS(P1R,P1I) C----------------------------------------------------------------------- C THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU C GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT C P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR C PREMATURELY. C----------------------------------------------------------------------- ARG = (AP2+AP2)/(AP1*TOL) TEST1 = SQRT(ARG) TEST = TEST1 RAP1 = 1.0D0/AP1 P1R = P1R*RAP1 P1I = P1I*RAP1 P2R = P2R*RAP1 P2I = P2I*RAP1 AP2 = AP2*RAP1 10 CONTINUE K = K + 1 AP1 = AP2 PTR = P2R PTI = P2I P2R = P1R - (T1R*PTR-T1I*PTI) P2I = P1I - (T1R*PTI+T1I*PTR) P1R = PTR P1I = PTI T1R = T1R + RZR T1I = T1I + RZI AP2 = ZABS(P2R,P2I) IF (AP1.LE.TEST) GO TO 10 IF (ITIME.EQ.2) GO TO 20 AK = ZABS(T1R,T1I)*0.5D0 FLAM = AK + SQRT(AK*AK-1.0D0) RHO = MIN(AP2/AP1,FLAM) TEST = TEST1*SQRT(RHO/(RHO*RHO-1.0D0)) ITIME = 2 GO TO 10 20 CONTINUE KK = K + 1 - ID AK = KK T1R = AK T1I = CZEROI DFNU = FNU + (N-1) P1R = 1.0D0/AP2 P1I = CZEROI P2R = CZEROR P2I = CZEROI DO 30 I=1,KK PTR = P1R PTI = P1I RAP1 = DFNU + T1R TTR = RZR*RAP1 TTI = RZI*RAP1 P1R = (PTR*TTR-PTI*TTI) + P2R P1I = (PTR*TTI+PTI*TTR) + P2I P2R = PTR P2I = PTI T1R = T1R - CONER 30 CONTINUE IF (P1R.NE.CZEROR .OR. P1I.NE.CZEROI) GO TO 40 P1R = TOL P1I = TOL 40 CONTINUE CALL ZDIV(P2R, P2I, P1R, P1I, CYR(N), CYI(N)) IF (N.EQ.1) RETURN K = N - 1 AK = K T1R = AK T1I = CZEROI CDFNUR = FNU*RZR CDFNUI = FNU*RZI DO 60 I=2,N PTR = CDFNUR + (T1R*RZR-T1I*RZI) + CYR(K+1) PTI = CDFNUI + (T1R*RZI+T1I*RZR) + CYI(K+1) AK = ZABS(PTR,PTI) IF (AK.NE.CZEROR) GO TO 50 PTR = TOL PTI = TOL AK = TOL*RT2 50 CONTINUE RAK = CONER/AK CYR(K) = RAK*PTR*RAK CYI(K) = -RAK*PTI*RAK T1R = T1R - CONER K = K - 1 60 CONTINUE RETURN END SUBROUTINE ZS1S2 (ZRR, ZRI, S1R, S1I, S2R, S2I, NZ, ASCLE, ALIM, + IUF) C***BEGIN PROLOGUE ZS1S2 C***SUBSIDIARY C***PURPOSE Subsidiary to ZAIRY and ZBESK C***LIBRARY SLATEC C***TYPE ALL (CS1S2-A, ZS1S2-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZS1S2 TESTS FOR A POSSIBLE UNDERFLOW RESULTING FROM THE C ADDITION OF THE I AND K FUNCTIONS IN THE ANALYTIC CON- C TINUATION FORMULA WHERE S1=K FUNCTION AND S2=I FUNCTION. C ON KODE=1 THE I AND K FUNCTIONS ARE DIFFERENT ORDERS OF C MAGNITUDE, BUT FOR KODE=2 THEY CAN BE OF THE SAME ORDER C OF MAGNITUDE AND THE MAXIMUM MUST BE AT LEAST ONE C PRECISION ABOVE THE UNDERFLOW LIMIT. C C***SEE ALSO ZAIRY, ZBESK C***ROUTINES CALLED ZABS, ZEXP, ZLOG C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZS1S2 C COMPLEX CZERO,C1,S1,S1D,S2,ZR DOUBLE PRECISION AA, ALIM, ALN, ASCLE, AS1, AS2, C1I, C1R, S1DI, * S1DR, S1I, S1R, S2I, S2R, ZEROI, ZEROR, ZRI, ZRR, ZABS INTEGER IUF, IDUM, NZ EXTERNAL ZABS DATA ZEROR,ZEROI / 0.0D0 , 0.0D0 / C***FIRST EXECUTABLE STATEMENT ZS1S2 NZ = 0 AS1 = ZABS(S1R,S1I) AS2 = ZABS(S2R,S2I) IF (S1R.EQ.0.0D0 .AND. S1I.EQ.0.0D0) GO TO 10 IF (AS1.EQ.0.0D0) GO TO 10 ALN = -ZRR - ZRR + LOG(AS1) S1DR = S1R S1DI = S1I S1R = ZEROR S1I = ZEROI AS1 = ZEROR IF (ALN.LT.(-ALIM)) GO TO 10 CALL ZLOG(S1DR, S1DI, C1R, C1I, IDUM) C1R = C1R - ZRR - ZRR C1I = C1I - ZRI - ZRI CALL ZEXP(C1R, C1I, S1R, S1I) AS1 = ZABS(S1R,S1I) IUF = IUF + 1 10 CONTINUE AA = MAX(AS1,AS2) IF (AA.GT.ASCLE) RETURN S1R = ZEROR S1I = ZEROI S2R = ZEROR S2I = ZEROI NZ = 1 IUF = 0 RETURN END SUBROUTINE ZSERI (ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, + ALIM) C***BEGIN PROLOGUE ZSERI C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESI and ZBESK C***LIBRARY SLATEC C***TYPE ALL (CSERI-A, ZSERI-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY C MEANS OF THE POWER SERIES FOR LARGE ABS(Z) IN THE C REGION ABS(Z).LE.2*SQRT(FNU+1). NZ=0 IS A NORMAL RETURN. C NZ.GT.0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO C DUE TO UNDERFLOW. NZ.LT.0 MEANS UNDERFLOW OCCURRED, BUT THE C CONDITION ABS(Z).LE.2*SQRT(FNU+1) WAS VIOLATED AND THE C COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ). C C***SEE ALSO ZBESI, ZBESK C***ROUTINES CALLED D1MACH, DGAMLN, ZABS, ZDIV, ZLOG, ZMLT, ZUCHK C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZSERI C COMPLEX AK1,CK,COEF,CONE,CRSC,CSCL,CZ,CZERO,HZ,RZ,S1,S2,Y,Z DOUBLE PRECISION AA, ACZ, AK, AK1I, AK1R, ALIM, ARM, ASCLE, ATOL, * AZ, CKI, CKR, COEFI, COEFR, CONEI, CONER, CRSCR, CZI, CZR, DFNU, * ELIM, FNU, FNUP, HZI, HZR, RAZ, RS, RTR1, RZI, RZR, S, SS, STI, * STR, S1I, S1R, S2I, S2R, TOL, YI, YR, WI, WR, ZEROI, ZEROR, ZI, * ZR, DGAMLN, D1MACH, ZABS INTEGER I, IB, IDUM, IFLAG, IL, K, KODE, L, M, N, NN, NZ, NW DIMENSION YR(N), YI(N), WR(2), WI(2) EXTERNAL ZABS DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 / C***FIRST EXECUTABLE STATEMENT ZSERI NZ = 0 AZ = ZABS(ZR,ZI) IF (AZ.EQ.0.0D0) GO TO 160 ARM = 1.0D+3*D1MACH(1) RTR1 = SQRT(ARM) CRSCR = 1.0D0 IFLAG = 0 IF (AZ.LT.ARM) GO TO 150 HZR = 0.5D0*ZR HZI = 0.5D0*ZI CZR = ZEROR CZI = ZEROI IF (AZ.LE.RTR1) GO TO 10 CALL ZMLT(HZR, HZI, HZR, HZI, CZR, CZI) 10 CONTINUE ACZ = ZABS(CZR,CZI) NN = N CALL ZLOG(HZR, HZI, CKR, CKI, IDUM) 20 CONTINUE DFNU = FNU + (NN-1) FNUP = DFNU + 1.0D0 C----------------------------------------------------------------------- C UNDERFLOW TEST C----------------------------------------------------------------------- AK1R = CKR*DFNU AK1I = CKI*DFNU AK = DGAMLN(FNUP,IDUM) AK1R = AK1R - AK IF (KODE.EQ.2) AK1R = AK1R - ZR IF (AK1R.GT.(-ELIM)) GO TO 40 30 CONTINUE NZ = NZ + 1 YR(NN) = ZEROR YI(NN) = ZEROI IF (ACZ.GT.DFNU) GO TO 190 NN = NN - 1 IF (NN.EQ.0) RETURN GO TO 20 40 CONTINUE IF (AK1R.GT.(-ALIM)) GO TO 50 IFLAG = 1 SS = 1.0D0/TOL CRSCR = TOL ASCLE = ARM*SS 50 CONTINUE AA = EXP(AK1R) IF (IFLAG.EQ.1) AA = AA*SS COEFR = AA*COS(AK1I) COEFI = AA*SIN(AK1I) ATOL = TOL*ACZ/FNUP IL = MIN(2,NN) DO 90 I=1,IL DFNU = FNU + (NN-I) FNUP = DFNU + 1.0D0 S1R = CONER S1I = CONEI IF (ACZ.LT.TOL*FNUP) GO TO 70 AK1R = CONER AK1I = CONEI AK = FNUP + 2.0D0 S = FNUP AA = 2.0D0 60 CONTINUE RS = 1.0D0/S STR = AK1R*CZR - AK1I*CZI STI = AK1R*CZI + AK1I*CZR AK1R = STR*RS AK1I = STI*RS S1R = S1R + AK1R S1I = S1I + AK1I S = S + AK AK = AK + 2.0D0 AA = AA*ACZ*RS IF (AA.GT.ATOL) GO TO 60 70 CONTINUE S2R = S1R*COEFR - S1I*COEFI S2I = S1R*COEFI + S1I*COEFR WR(I) = S2R WI(I) = S2I IF (IFLAG.EQ.0) GO TO 80 CALL ZUCHK(S2R, S2I, NW, ASCLE, TOL) IF (NW.NE.0) GO TO 30 80 CONTINUE M = NN - I + 1 YR(M) = S2R*CRSCR YI(M) = S2I*CRSCR IF (I.EQ.IL) GO TO 90 CALL ZDIV(COEFR, COEFI, HZR, HZI, STR, STI) COEFR = STR*DFNU COEFI = STI*DFNU 90 CONTINUE IF (NN.LE.2) RETURN K = NN - 2 AK = K RAZ = 1.0D0/AZ STR = ZR*RAZ STI = -ZI*RAZ RZR = (STR+STR)*RAZ RZI = (STI+STI)*RAZ IF (IFLAG.EQ.1) GO TO 120 IB = 3 100 CONTINUE DO 110 I=IB,NN YR(K) = (AK+FNU)*(RZR*YR(K+1)-RZI*YI(K+1)) + YR(K+2) YI(K) = (AK+FNU)*(RZR*YI(K+1)+RZI*YR(K+1)) + YI(K+2) AK = AK - 1.0D0 K = K - 1 110 CONTINUE RETURN C----------------------------------------------------------------------- C RECUR BACKWARD WITH SCALED VALUES C----------------------------------------------------------------------- 120 CONTINUE C----------------------------------------------------------------------- C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE C UNDERFLOW LIMIT = ASCLE = D1MACH(1)*SS*1.0D+3 C----------------------------------------------------------------------- S1R = WR(1) S1I = WI(1) S2R = WR(2) S2I = WI(2) DO 130 L=3,NN CKR = S2R CKI = S2I S2R = S1R + (AK+FNU)*(RZR*CKR-RZI*CKI) S2I = S1I + (AK+FNU)*(RZR*CKI+RZI*CKR) S1R = CKR S1I = CKI CKR = S2R*CRSCR CKI = S2I*CRSCR YR(K) = CKR YI(K) = CKI AK = AK - 1.0D0 K = K - 1 IF (ZABS(CKR,CKI).GT.ASCLE) GO TO 140 130 CONTINUE RETURN 140 CONTINUE IB = L + 1 IF (IB.GT.NN) RETURN GO TO 100 150 CONTINUE NZ = N IF (FNU.EQ.0.0D0) NZ = NZ - 1 160 CONTINUE YR(1) = ZEROR YI(1) = ZEROI IF (FNU.NE.0.0D0) GO TO 170 YR(1) = CONER YI(1) = CONEI 170 CONTINUE IF (N.EQ.1) RETURN DO 180 I=2,N YR(I) = ZEROR YI(I) = ZEROI 180 CONTINUE RETURN C----------------------------------------------------------------------- C RETURN WITH NZ.LT.0 IF ABS(Z*Z/4).GT.FNU+N-NZ-1 COMPLETE C THE CALCULATION IN CBINU WITH N=N-ABS(NZ) C----------------------------------------------------------------------- 190 CONTINUE NZ = -NZ RETURN END SUBROUTINE ZSHCH (ZR, ZI, CSHR, CSHI, CCHR, CCHI) C***BEGIN PROLOGUE ZSHCH C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESH and ZBESK C***LIBRARY SLATEC C***TYPE ALL (CSHCH-A, ZSHCH-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZSHCH COMPUTES THE COMPLEX HYPERBOLIC FUNCTIONS CSH=SINH(X+I*Y) C AND CCH=COSH(X+I*Y), WHERE I**2=-1. C C***SEE ALSO ZBESH, ZBESK C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZSHCH C DOUBLE PRECISION CCHI, CCHR, CH, CN, CSHI, CSHR, SH, SN, ZI, ZR C***FIRST EXECUTABLE STATEMENT ZSHCH SH = SINH(ZR) CH = COSH(ZR) SN = SIN(ZI) CN = COS(ZI) CSHR = SH*CN CSHI = CH*SN CCHR = CH*CN CCHI = SH*SN RETURN END SUBROUTINE ZSQRT (AR, AI, BR, BI) C***BEGIN PROLOGUE ZSQRT C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESH, ZBESI, ZBESJ, ZBESK, ZBESY, ZAIRY and C ZBIRY C***LIBRARY SLATEC C***TYPE ALL (ZSQRT-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C DOUBLE PRECISION COMPLEX SQUARE ROOT, B=CSQRT(A) C C***SEE ALSO ZAIRY, ZBESH, ZBESI, ZBESJ, ZBESK, ZBESY, ZBIRY C***ROUTINES CALLED ZABS C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZSQRT DOUBLE PRECISION AR, AI, BR, BI, ZM, DTHETA, DPI, DRT DOUBLE PRECISION ZABS EXTERNAL ZABS DATA DRT , DPI / 7.071067811865475244008443621D-1, 1 3.141592653589793238462643383D+0/ C***FIRST EXECUTABLE STATEMENT ZSQRT ZM = ZABS(AR,AI) ZM = SQRT(ZM) IF (AR.EQ.0.0D+0) GO TO 10 IF (AI.EQ.0.0D+0) GO TO 20 DTHETA = DATAN(AI/AR) IF (DTHETA.LE.0.0D+0) GO TO 40 IF (AR.LT.0.0D+0) DTHETA = DTHETA - DPI GO TO 50 10 IF (AI.GT.0.0D+0) GO TO 60 IF (AI.LT.0.0D+0) GO TO 70 BR = 0.0D+0 BI = 0.0D+0 RETURN 20 IF (AR.GT.0.0D+0) GO TO 30 BR = 0.0D+0 BI = SQRT(ABS(AR)) RETURN 30 BR = SQRT(AR) BI = 0.0D+0 RETURN 40 IF (AR.LT.0.0D+0) DTHETA = DTHETA + DPI 50 DTHETA = DTHETA*0.5D+0 BR = ZM*COS(DTHETA) BI = ZM*SIN(DTHETA) RETURN 60 BR = ZM*DRT BI = ZM*DRT RETURN 70 BR = ZM*DRT BI = -ZM*DRT RETURN END SUBROUTINE ZUCHK (YR, YI, NZ, ASCLE, TOL) C***BEGIN PROLOGUE ZUCHK C***SUBSIDIARY C***PURPOSE Subsidiary to SERI, ZUOIK, ZUNK1, ZUNK2, ZUNI1, ZUNI2 and C ZKSCL C***LIBRARY SLATEC C***TYPE ALL (CUCHK-A, ZUCHK-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C Y ENTERS AS A SCALED QUANTITY WHOSE MAGNITUDE IS GREATER THAN C EXP(-ALIM)=ASCLE=1.0E+3*D1MACH(1)/TOL. THE TEST IS MADE TO SEE C IF THE MAGNITUDE OF THE REAL OR IMAGINARY PART WOULD UNDERFLOW C WHEN Y IS SCALED (BY TOL) TO ITS PROPER VALUE. Y IS ACCEPTED C IF THE UNDERFLOW IS AT LEAST ONE PRECISION BELOW THE MAGNITUDE C OF THE LARGEST COMPONENT; OTHERWISE THE PHASE ANGLE DOES NOT HAVE C ABSOLUTE ACCURACY AND AN UNDERFLOW IS ASSUMED. C C***SEE ALSO SERI, ZKSCL, ZUNI1, ZUNI2, ZUNK1, ZUNK2, ZUOIK C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C ?????? DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZUCHK C C COMPLEX Y DOUBLE PRECISION ASCLE, SS, ST, TOL, WR, WI, YR, YI INTEGER NZ C***FIRST EXECUTABLE STATEMENT ZUCHK NZ = 0 WR = ABS(YR) WI = ABS(YI) ST = MIN(WR,WI) IF (ST.GT.ASCLE) RETURN SS = MAX(WR,WI) ST = ST/TOL IF (SS.LT.ST) NZ = 1 RETURN END SUBROUTINE ZUNHJ (ZR, ZI, FNU, IPMTR, TOL, PHIR, PHII, ARGR, ARGI, + ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) C***BEGIN PROLOGUE ZUNHJ C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESI and ZBESK C***LIBRARY SLATEC C***TYPE ALL (CUNHJ-A, ZUNHJ-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C REFERENCES C HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ AND I.A. C STEGUN, AMS55, NATIONAL BUREAU OF STANDARDS, 1965, CHAPTER 9. C C ASYMPTOTICS AND SPECIAL FUNCTIONS BY F.W.J. OLVER, ACADEMIC C PRESS, N.Y., 1974, PAGE 420 C C ABSTRACT C ZUNHJ COMPUTES PARAMETERS FOR BESSEL FUNCTIONS C(FNU,Z) = C J(FNU,Z), Y(FNU,Z) OR H(I,FNU,Z) I=1,2 FOR LARGE ORDERS FNU C BY MEANS OF THE UNIFORM ASYMPTOTIC EXPANSION C C C(FNU,Z)=C1*PHI*( ASUM*AIRY(ARG) + C2*BSUM*DAIRY(ARG) ) C C FOR PROPER CHOICES OF C1, C2, AIRY AND DAIRY WHERE AIRY IS C AN AIRY FUNCTION AND DAIRY IS ITS DERIVATIVE. C C (2/3)*FNU*ZETA**1.5 = ZETA1-ZETA2, C C ZETA1=0.5*FNU*CLOG((1+W)/(1-W)), ZETA2=FNU*W FOR SCALING C PURPOSES IN AIRY FUNCTIONS FROM CAIRY OR CBIRY. C C MCONJ=SIGN OF AIMAG(Z), BUT IS AMBIGUOUS WHEN Z IS REAL AND C MUST BE SPECIFIED. IPMTR=0 RETURNS ALL PARAMETERS. IPMTR= C 1 COMPUTES ALL EXCEPT ASUM AND BSUM. C C***SEE ALSO ZBESI, ZBESK C***ROUTINES CALLED D1MACH, ZABS, ZDIV, ZLOG, ZSQRT C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZUNHJ C COMPLEX ARG,ASUM,BSUM,CFNU,CONE,CR,CZERO,DR,P,PHI,PRZTH,PTFN, C *RFN13,RTZTA,RZTH,SUMA,SUMB,TFN,T2,UP,W,W2,Z,ZA,ZB,ZC,ZETA,ZETA1, C *ZETA2,ZTH DOUBLE PRECISION ALFA, ANG, AP, AR, ARGI, ARGR, ASUMI, ASUMR, * ATOL, AW2, AZTH, BETA, BR, BSUMI, BSUMR, BTOL, C, CONEI, CONER, * CRI, CRR, DRI, DRR, EX1, EX2, FNU, FN13, FN23, GAMA, GPI, HPI, * PHII, PHIR, PI, PP, PR, PRZTHI, PRZTHR, PTFNI, PTFNR, RAW, RAW2, * RAZTH, RFNU, RFNU2, RFN13, RTZTI, RTZTR, RZTHI, RZTHR, STI, STR, * SUMAI, SUMAR, SUMBI, SUMBR, TEST, TFNI, TFNR, THPI, TOL, TZAI, * TZAR, T2I, T2R, UPI, UPR, WI, WR, W2I, W2R, ZAI, ZAR, ZBI, ZBR, * ZCI, ZCR, ZEROI, ZEROR, ZETAI, ZETAR, ZETA1I, ZETA1R, ZETA2I, * ZETA2R, ZI, ZR, ZTHI, ZTHR, ZABS, AC, D1MACH INTEGER IAS, IBS, IPMTR, IS, J, JR, JU, K, KMAX, KP1, KS, L, LR, * LRP1, L1, L2, M, IDUM DIMENSION AR(14), BR(14), C(105), ALFA(180), BETA(210), GAMA(30), * AP(30), PR(30), PI(30), UPR(14), UPI(14), CRR(14), CRI(14), * DRR(14), DRI(14) EXTERNAL ZABS DATA AR(1), AR(2), AR(3), AR(4), AR(5), AR(6), AR(7), AR(8), 1 AR(9), AR(10), AR(11), AR(12), AR(13), AR(14)/ 2 1.00000000000000000D+00, 1.04166666666666667D-01, 3 8.35503472222222222D-02, 1.28226574556327160D-01, 4 2.91849026464140464D-01, 8.81627267443757652D-01, 5 3.32140828186276754D+00, 1.49957629868625547D+01, 6 7.89230130115865181D+01, 4.74451538868264323D+02, 7 3.20749009089066193D+03, 2.40865496408740049D+04, 8 1.98923119169509794D+05, 1.79190200777534383D+06/ DATA BR(1), BR(2), BR(3), BR(4), BR(5), BR(6), BR(7), BR(8), 1 BR(9), BR(10), BR(11), BR(12), BR(13), BR(14)/ 2 1.00000000000000000D+00, -1.45833333333333333D-01, 3 -9.87413194444444444D-02, -1.43312053915895062D-01, 4 -3.17227202678413548D-01, -9.42429147957120249D-01, 5 -3.51120304082635426D+00, -1.57272636203680451D+01, 6 -8.22814390971859444D+01, -4.92355370523670524D+02, 7 -3.31621856854797251D+03, -2.48276742452085896D+04, 8 -2.04526587315129788D+05, -1.83844491706820990D+06/ 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 1.00000000000000000D+00, -2.08333333333333333D-01, 4 1.25000000000000000D-01, 3.34201388888888889D-01, 5 -4.01041666666666667D-01, 7.03125000000000000D-02, 6 -1.02581259645061728D+00, 1.84646267361111111D+00, 7 -8.91210937500000000D-01, 7.32421875000000000D-02, 8 4.66958442342624743D+00, -1.12070026162229938D+01, 9 8.78912353515625000D+00, -2.36408691406250000D+00, A 1.12152099609375000D-01, -2.82120725582002449D+01, B 8.46362176746007346D+01, -9.18182415432400174D+01, C 4.25349987453884549D+01, -7.36879435947963170D+00, D 2.27108001708984375D-01, 2.12570130039217123D+02, E -7.65252468141181642D+02, 1.05999045252799988D+03/ 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 -6.99579627376132541D+02, 2.18190511744211590D+02, 4 -2.64914304869515555D+01, 5.72501420974731445D-01, 5 -1.91945766231840700D+03, 8.06172218173730938D+03, 6 -1.35865500064341374D+04, 1.16553933368645332D+04, 7 -5.30564697861340311D+03, 1.20090291321635246D+03, 8 -1.08090919788394656D+02, 1.72772750258445740D+00, 9 2.02042913309661486D+04, -9.69805983886375135D+04, A 1.92547001232531532D+05, -2.03400177280415534D+05, B 1.22200464983017460D+05, -4.11926549688975513D+04, C 7.10951430248936372D+03, -4.93915304773088012D+02, D 6.07404200127348304D+00, -2.42919187900551333D+05, E 1.31176361466297720D+06, -2.99801591853810675D+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), C(66), C(67), C(68), C(69), C(70), C(71), C(72)/ 3 3.76327129765640400D+06, -2.81356322658653411D+06, 4 1.26836527332162478D+06, -3.31645172484563578D+05, 5 4.52187689813627263D+04, -2.49983048181120962D+03, 6 2.43805296995560639D+01, 3.28446985307203782D+06, 7 -1.97068191184322269D+07, 5.09526024926646422D+07, 8 -7.41051482115326577D+07, 6.63445122747290267D+07, 9 -3.75671766607633513D+07, 1.32887671664218183D+07, A -2.78561812808645469D+06, 3.08186404612662398D+05, B -1.38860897537170405D+04, 1.10017140269246738D+02, C -4.93292536645099620D+07, 3.25573074185765749D+08, D -9.39462359681578403D+08, 1.55359689957058006D+09, E -1.62108055210833708D+09, 1.10684281682301447D+09/ DATA C(73), C(74), C(75), C(76), C(77), C(78), C(79), C(80), 1 C(81), C(82), C(83), C(84), C(85), C(86), C(87), C(88), 2 C(89), C(90), C(91), C(92), C(93), C(94), C(95), C(96)/ 3 -4.95889784275030309D+08, 1.42062907797533095D+08, 4 -2.44740627257387285D+07, 2.24376817792244943D+06, 5 -8.40054336030240853D+04, 5.51335896122020586D+02, 6 8.14789096118312115D+08, -5.86648149205184723D+09, 7 1.86882075092958249D+10, -3.46320433881587779D+10, 8 4.12801855797539740D+10, -3.30265997498007231D+10, 9 1.79542137311556001D+10, -6.56329379261928433D+09, A 1.55927986487925751D+09, -2.25105661889415278D+08, B 1.73951075539781645D+07, -5.49842327572288687D+05, C 3.03809051092238427D+03, -1.46792612476956167D+10, D 1.14498237732025810D+11, -3.99096175224466498D+11, E 8.19218669548577329D+11, -1.09837515608122331D+12/ DATA C(97), C(98), C(99), C(100), C(101), C(102), C(103), C(104), 1 C(105)/ 2 1.00815810686538209D+12, -6.45364869245376503D+11, 3 2.87900649906150589D+11, -8.78670721780232657D+10, 4 1.76347306068349694D+10, -2.16716498322379509D+09, 5 1.43157876718888981D+08, -3.87183344257261262D+06, 6 1.82577554742931747D+04/ DATA ALFA(1), ALFA(2), ALFA(3), ALFA(4), ALFA(5), ALFA(6), 1 ALFA(7), ALFA(8), ALFA(9), ALFA(10), ALFA(11), ALFA(12), 2 ALFA(13), ALFA(14), ALFA(15), ALFA(16), ALFA(17), ALFA(18), 3 ALFA(19), ALFA(20), ALFA(21), ALFA(22)/ 4 -4.44444444444444444D-03, -9.22077922077922078D-04, 5 -8.84892884892884893D-05, 1.65927687832449737D-04, 6 2.46691372741792910D-04, 2.65995589346254780D-04, 7 2.61824297061500945D-04, 2.48730437344655609D-04, 8 2.32721040083232098D-04, 2.16362485712365082D-04, 9 2.00738858762752355D-04, 1.86267636637545172D-04, A 1.73060775917876493D-04, 1.61091705929015752D-04, B 1.50274774160908134D-04, 1.40503497391269794D-04, C 1.31668816545922806D-04, 1.23667445598253261D-04, D 1.16405271474737902D-04, 1.09798298372713369D-04, E 1.03772410422992823D-04, 9.82626078369363448D-05/ DATA ALFA(23), ALFA(24), ALFA(25), ALFA(26), ALFA(27), ALFA(28), 1 ALFA(29), ALFA(30), ALFA(31), ALFA(32), ALFA(33), ALFA(34), 2 ALFA(35), ALFA(36), ALFA(37), ALFA(38), ALFA(39), ALFA(40), 3 ALFA(41), ALFA(42), ALFA(43), ALFA(44)/ 4 9.32120517249503256D-05, 8.85710852478711718D-05, 5 8.42963105715700223D-05, 8.03497548407791151D-05, 6 7.66981345359207388D-05, 7.33122157481777809D-05, 7 7.01662625163141333D-05, 6.72375633790160292D-05, 8 6.93735541354588974D-04, 2.32241745182921654D-04, 9 -1.41986273556691197D-05, -1.16444931672048640D-04, A -1.50803558053048762D-04, -1.55121924918096223D-04, B -1.46809756646465549D-04, -1.33815503867491367D-04, C -1.19744975684254051D-04, -1.06184319207974020D-04, D -9.37699549891194492D-05, -8.26923045588193274D-05, E -7.29374348155221211D-05, -6.44042357721016283D-05/ DATA ALFA(45), ALFA(46), ALFA(47), ALFA(48), ALFA(49), ALFA(50), 1 ALFA(51), ALFA(52), ALFA(53), ALFA(54), ALFA(55), ALFA(56), 2 ALFA(57), ALFA(58), ALFA(59), ALFA(60), ALFA(61), ALFA(62), 3 ALFA(63), ALFA(64), ALFA(65), ALFA(66)/ 4 -5.69611566009369048D-05, -5.04731044303561628D-05, 5 -4.48134868008882786D-05, -3.98688727717598864D-05, 6 -3.55400532972042498D-05, -3.17414256609022480D-05, 7 -2.83996793904174811D-05, -2.54522720634870566D-05, 8 -2.28459297164724555D-05, -2.05352753106480604D-05, 9 -1.84816217627666085D-05, -1.66519330021393806D-05, A -1.50179412980119482D-05, -1.35554031379040526D-05, B -1.22434746473858131D-05, -1.10641884811308169D-05, C -3.54211971457743841D-04, -1.56161263945159416D-04, D 3.04465503594936410D-05, 1.30198655773242693D-04, E 1.67471106699712269D-04, 1.70222587683592569D-04/ DATA ALFA(67), ALFA(68), ALFA(69), ALFA(70), ALFA(71), ALFA(72), 1 ALFA(73), ALFA(74), ALFA(75), ALFA(76), ALFA(77), ALFA(78), 2 ALFA(79), ALFA(80), ALFA(81), ALFA(82), ALFA(83), ALFA(84), 3 ALFA(85), ALFA(86), ALFA(87), ALFA(88)/ 4 1.56501427608594704D-04, 1.36339170977445120D-04, 5 1.14886692029825128D-04, 9.45869093034688111D-05, 6 7.64498419250898258D-05, 6.07570334965197354D-05, 7 4.74394299290508799D-05, 3.62757512005344297D-05, 8 2.69939714979224901D-05, 1.93210938247939253D-05, 9 1.30056674793963203D-05, 7.82620866744496661D-06, A 3.59257485819351583D-06, 1.44040049814251817D-07, B -2.65396769697939116D-06, -4.91346867098485910D-06, C -6.72739296091248287D-06, -8.17269379678657923D-06, D -9.31304715093561232D-06, -1.02011418798016441D-05, E -1.08805962510592880D-05, -1.13875481509603555D-05/ DATA ALFA(89), ALFA(90), ALFA(91), ALFA(92), ALFA(93), ALFA(94), 1 ALFA(95), ALFA(96), ALFA(97), ALFA(98), ALFA(99), ALFA(100), 2 ALFA(101), ALFA(102), ALFA(103), ALFA(104), ALFA(105), 3 ALFA(106), ALFA(107), ALFA(108), ALFA(109), ALFA(110)/ 4 -1.17519675674556414D-05, -1.19987364870944141D-05, 5 3.78194199201772914D-04, 2.02471952761816167D-04, 6 -6.37938506318862408D-05, -2.38598230603005903D-04, 7 -3.10916256027361568D-04, -3.13680115247576316D-04, 8 -2.78950273791323387D-04, -2.28564082619141374D-04, 9 -1.75245280340846749D-04, -1.25544063060690348D-04, A -8.22982872820208365D-05, -4.62860730588116458D-05, B -1.72334302366962267D-05, 5.60690482304602267D-06, C 2.31395443148286800D-05, 3.62642745856793957D-05, D 4.58006124490188752D-05, 5.24595294959114050D-05, E 5.68396208545815266D-05, 5.94349820393104052D-05/ DATA ALFA(111), ALFA(112), ALFA(113), ALFA(114), ALFA(115), 1 ALFA(116), ALFA(117), ALFA(118), ALFA(119), ALFA(120), 2 ALFA(121), ALFA(122), ALFA(123), ALFA(124), ALFA(125), 3 ALFA(126), ALFA(127), ALFA(128), ALFA(129), ALFA(130)/ 4 6.06478527578421742D-05, 6.08023907788436497D-05, 5 6.01577894539460388D-05, 5.89199657344698500D-05, 6 5.72515823777593053D-05, 5.52804375585852577D-05, 7 5.31063773802880170D-05, 5.08069302012325706D-05, 8 4.84418647620094842D-05, 4.60568581607475370D-05, 9 -6.91141397288294174D-04, -4.29976633058871912D-04, A 1.83067735980039018D-04, 6.60088147542014144D-04, B 8.75964969951185931D-04, 8.77335235958235514D-04, C 7.49369585378990637D-04, 5.63832329756980918D-04, D 3.68059319971443156D-04, 1.88464535514455599D-04/ DATA ALFA(131), ALFA(132), ALFA(133), ALFA(134), ALFA(135), 1 ALFA(136), ALFA(137), ALFA(138), ALFA(139), ALFA(140), 2 ALFA(141), ALFA(142), ALFA(143), ALFA(144), ALFA(145), 3 ALFA(146), ALFA(147), ALFA(148), ALFA(149), ALFA(150)/ 4 3.70663057664904149D-05, -8.28520220232137023D-05, 5 -1.72751952869172998D-04, -2.36314873605872983D-04, 6 -2.77966150694906658D-04, -3.02079514155456919D-04, 7 -3.12594712643820127D-04, -3.12872558758067163D-04, 8 -3.05678038466324377D-04, -2.93226470614557331D-04, 9 -2.77255655582934777D-04, -2.59103928467031709D-04, A -2.39784014396480342D-04, -2.20048260045422848D-04, B -2.00443911094971498D-04, -1.81358692210970687D-04, C -1.63057674478657464D-04, -1.45712672175205844D-04, D -1.29425421983924587D-04, -1.14245691942445952D-04/ DATA ALFA(151), ALFA(152), ALFA(153), ALFA(154), ALFA(155), 1 ALFA(156), ALFA(157), ALFA(158), ALFA(159), ALFA(160), 2 ALFA(161), ALFA(162), ALFA(163), ALFA(164), ALFA(165), 3 ALFA(166), ALFA(167), ALFA(168), ALFA(169), ALFA(170)/ 4 1.92821964248775885D-03, 1.35592576302022234D-03, 5 -7.17858090421302995D-04, -2.58084802575270346D-03, 6 -3.49271130826168475D-03, -3.46986299340960628D-03, 7 -2.82285233351310182D-03, -1.88103076404891354D-03, 8 -8.89531718383947600D-04, 3.87912102631035228D-06, 9 7.28688540119691412D-04, 1.26566373053457758D-03, A 1.62518158372674427D-03, 1.83203153216373172D-03, B 1.91588388990527909D-03, 1.90588846755546138D-03, C 1.82798982421825727D-03, 1.70389506421121530D-03, D 1.55097127171097686D-03, 1.38261421852276159D-03/ DATA ALFA(171), ALFA(172), ALFA(173), ALFA(174), ALFA(175), 1 ALFA(176), ALFA(177), ALFA(178), ALFA(179), ALFA(180)/ 2 1.20881424230064774D-03, 1.03676532638344962D-03, 3 8.71437918068619115D-04, 7.16080155297701002D-04, 4 5.72637002558129372D-04, 4.42089819465802277D-04, 5 3.24724948503090564D-04, 2.20342042730246599D-04, 6 1.28412898401353882D-04, 4.82005924552095464D-05/ DATA BETA(1), BETA(2), BETA(3), BETA(4), BETA(5), BETA(6), 1 BETA(7), BETA(8), BETA(9), BETA(10), BETA(11), BETA(12), 2 BETA(13), BETA(14), BETA(15), BETA(16), BETA(17), BETA(18), 3 BETA(19), BETA(20), BETA(21), BETA(22)/ 4 1.79988721413553309D-02, 5.59964911064388073D-03, 5 2.88501402231132779D-03, 1.80096606761053941D-03, 6 1.24753110589199202D-03, 9.22878876572938311D-04, 7 7.14430421727287357D-04, 5.71787281789704872D-04, 8 4.69431007606481533D-04, 3.93232835462916638D-04, 9 3.34818889318297664D-04, 2.88952148495751517D-04, A 2.52211615549573284D-04, 2.22280580798883327D-04, B 1.97541838033062524D-04, 1.76836855019718004D-04, C 1.59316899661821081D-04, 1.44347930197333986D-04, D 1.31448068119965379D-04, 1.20245444949302884D-04, E 1.10449144504599392D-04, 1.01828770740567258D-04/ DATA BETA(23), BETA(24), BETA(25), BETA(26), BETA(27), BETA(28), 1 BETA(29), BETA(30), BETA(31), BETA(32), BETA(33), BETA(34), 2 BETA(35), BETA(36), BETA(37), BETA(38), BETA(39), BETA(40), 3 BETA(41), BETA(42), BETA(43), BETA(44)/ 4 9.41998224204237509D-05, 8.74130545753834437D-05, 5 8.13466262162801467D-05, 7.59002269646219339D-05, 6 7.09906300634153481D-05, 6.65482874842468183D-05, 7 6.25146958969275078D-05, 5.88403394426251749D-05, 8 -1.49282953213429172D-03, -8.78204709546389328D-04, 9 -5.02916549572034614D-04, -2.94822138512746025D-04, A -1.75463996970782828D-04, -1.04008550460816434D-04, B -5.96141953046457895D-05, -3.12038929076098340D-05, C -1.26089735980230047D-05, -2.42892608575730389D-07, D 8.05996165414273571D-06, 1.36507009262147391D-05, E 1.73964125472926261D-05, 1.98672978842133780D-05/ DATA BETA(45), BETA(46), BETA(47), BETA(48), BETA(49), BETA(50), 1 BETA(51), BETA(52), BETA(53), BETA(54), BETA(55), BETA(56), 2 BETA(57), BETA(58), BETA(59), BETA(60), BETA(61), BETA(62), 3 BETA(63), BETA(64), BETA(65), BETA(66)/ 4 2.14463263790822639D-05, 2.23954659232456514D-05, 5 2.28967783814712629D-05, 2.30785389811177817D-05, 6 2.30321976080909144D-05, 2.28236073720348722D-05, 7 2.25005881105292418D-05, 2.20981015361991429D-05, 8 2.16418427448103905D-05, 2.11507649256220843D-05, 9 2.06388749782170737D-05, 2.01165241997081666D-05, A 1.95913450141179244D-05, 1.90689367910436740D-05, B 1.85533719641636667D-05, 1.80475722259674218D-05, C 5.52213076721292790D-04, 4.47932581552384646D-04, D 2.79520653992020589D-04, 1.52468156198446602D-04, E 6.93271105657043598D-05, 1.76258683069991397D-05/ DATA BETA(67), BETA(68), BETA(69), BETA(70), BETA(71), BETA(72), 1 BETA(73), BETA(74), BETA(75), BETA(76), BETA(77), BETA(78), 2 BETA(79), BETA(80), BETA(81), BETA(82), BETA(83), BETA(84), 3 BETA(85), BETA(86), BETA(87), BETA(88)/ 4 -1.35744996343269136D-05, -3.17972413350427135D-05, 5 -4.18861861696693365D-05, -4.69004889379141029D-05, 6 -4.87665447413787352D-05, -4.87010031186735069D-05, 7 -4.74755620890086638D-05, -4.55813058138628452D-05, 8 -4.33309644511266036D-05, -4.09230193157750364D-05, 9 -3.84822638603221274D-05, -3.60857167535410501D-05, A -3.37793306123367417D-05, -3.15888560772109621D-05, B -2.95269561750807315D-05, -2.75978914828335759D-05, C -2.58006174666883713D-05, -2.41308356761280200D-05, D -2.25823509518346033D-05, -2.11479656768912971D-05, E -1.98200638885294927D-05, -1.85909870801065077D-05/ DATA BETA(89), BETA(90), BETA(91), BETA(92), BETA(93), BETA(94), 1 BETA(95), BETA(96), BETA(97), BETA(98), BETA(99), BETA(100), 2 BETA(101), BETA(102), BETA(103), BETA(104), BETA(105), 3 BETA(106), BETA(107), BETA(108), BETA(109), BETA(110)/ 4 -1.74532699844210224D-05, -1.63997823854497997D-05, 5 -4.74617796559959808D-04, -4.77864567147321487D-04, 6 -3.20390228067037603D-04, -1.61105016119962282D-04, 7 -4.25778101285435204D-05, 3.44571294294967503D-05, 8 7.97092684075674924D-05, 1.03138236708272200D-04, 9 1.12466775262204158D-04, 1.13103642108481389D-04, A 1.08651634848774268D-04, 1.01437951597661973D-04, B 9.29298396593363896D-05, 8.40293133016089978D-05, C 7.52727991349134062D-05, 6.69632521975730872D-05, D 5.92564547323194704D-05, 5.22169308826975567D-05, E 4.58539485165360646D-05, 4.01445513891486808D-05/ DATA BETA(111), BETA(112), BETA(113), BETA(114), BETA(115), 1 BETA(116), BETA(117), BETA(118), BETA(119), BETA(120), 2 BETA(121), BETA(122), BETA(123), BETA(124), BETA(125), 3 BETA(126), BETA(127), BETA(128), BETA(129), BETA(130)/ 4 3.50481730031328081D-05, 3.05157995034346659D-05, 5 2.64956119950516039D-05, 2.29363633690998152D-05, 6 1.97893056664021636D-05, 1.70091984636412623D-05, 7 1.45547428261524004D-05, 1.23886640995878413D-05, 8 1.04775876076583236D-05, 8.79179954978479373D-06, 9 7.36465810572578444D-04, 8.72790805146193976D-04, A 6.22614862573135066D-04, 2.85998154194304147D-04, B 3.84737672879366102D-06, -1.87906003636971558D-04, C -2.97603646594554535D-04, -3.45998126832656348D-04, D -3.53382470916037712D-04, -3.35715635775048757D-04/ DATA BETA(131), BETA(132), BETA(133), BETA(134), BETA(135), 1 BETA(136), BETA(137), BETA(138), BETA(139), BETA(140), 2 BETA(141), BETA(142), BETA(143), BETA(144), BETA(145), 3 BETA(146), BETA(147), BETA(148), BETA(149), BETA(150)/ 4 -3.04321124789039809D-04, -2.66722723047612821D-04, 5 -2.27654214122819527D-04, -1.89922611854562356D-04, 6 -1.55058918599093870D-04, -1.23778240761873630D-04, 7 -9.62926147717644187D-05, -7.25178327714425337D-05, 8 -5.22070028895633801D-05, -3.50347750511900522D-05, 9 -2.06489761035551757D-05, -8.70106096849767054D-06, A 1.13698686675100290D-06, 9.16426474122778849D-06, B 1.56477785428872620D-05, 2.08223629482466847D-05, C 2.48923381004595156D-05, 2.80340509574146325D-05, D 3.03987774629861915D-05, 3.21156731406700616D-05/ DATA BETA(151), BETA(152), BETA(153), BETA(154), BETA(155), 1 BETA(156), BETA(157), BETA(158), BETA(159), BETA(160), 2 BETA(161), BETA(162), BETA(163), BETA(164), BETA(165), 3 BETA(166), BETA(167), BETA(168), BETA(169), BETA(170)/ 4 -1.80182191963885708D-03, -2.43402962938042533D-03, 5 -1.83422663549856802D-03, -7.62204596354009765D-04, 6 2.39079475256927218D-04, 9.49266117176881141D-04, 7 1.34467449701540359D-03, 1.48457495259449178D-03, 8 1.44732339830617591D-03, 1.30268261285657186D-03, 9 1.10351597375642682D-03, 8.86047440419791759D-04, A 6.73073208165665473D-04, 4.77603872856582378D-04, B 3.05991926358789362D-04, 1.60315694594721630D-04, C 4.00749555270613286D-05, -5.66607461635251611D-05, D -1.32506186772982638D-04, -1.90296187989614057D-04/ DATA BETA(171), BETA(172), BETA(173), BETA(174), BETA(175), 1 BETA(176), BETA(177), BETA(178), BETA(179), BETA(180), 2 BETA(181), BETA(182), BETA(183), BETA(184), BETA(185), 3 BETA(186), BETA(187), BETA(188), BETA(189), BETA(190)/ 4 -2.32811450376937408D-04, -2.62628811464668841D-04, 5 -2.82050469867598672D-04, -2.93081563192861167D-04, 6 -2.97435962176316616D-04, -2.96557334239348078D-04, 7 -2.91647363312090861D-04, -2.83696203837734166D-04, 8 -2.73512317095673346D-04, -2.61750155806768580D-04, 9 6.38585891212050914D-03, 9.62374215806377941D-03, A 7.61878061207001043D-03, 2.83219055545628054D-03, B -2.09841352012720090D-03, -5.73826764216626498D-03, C -7.70804244495414620D-03, -8.21011692264844401D-03, D -7.65824520346905413D-03, -6.47209729391045177D-03/ DATA BETA(191), BETA(192), BETA(193), BETA(194), BETA(195), 1 BETA(196), BETA(197), BETA(198), BETA(199), BETA(200), 2 BETA(201), BETA(202), BETA(203), BETA(204), BETA(205), 3 BETA(206), BETA(207), BETA(208), BETA(209), BETA(210)/ 4 -4.99132412004966473D-03, -3.45612289713133280D-03, 5 -2.01785580014170775D-03, -7.59430686781961401D-04, 6 2.84173631523859138D-04, 1.10891667586337403D-03, 7 1.72901493872728771D-03, 2.16812590802684701D-03, 8 2.45357710494539735D-03, 2.61281821058334862D-03, 9 2.67141039656276912D-03, 2.65203073395980430D-03, A 2.57411652877287315D-03, 2.45389126236094427D-03, B 2.30460058071795494D-03, 2.13684837686712662D-03, C 1.95896528478870911D-03, 1.77737008679454412D-03, D 1.59690280765839059D-03, 1.42111975664438546D-03/ DATA GAMA(1), GAMA(2), GAMA(3), GAMA(4), GAMA(5), GAMA(6), 1 GAMA(7), GAMA(8), GAMA(9), GAMA(10), GAMA(11), GAMA(12), 2 GAMA(13), GAMA(14), GAMA(15), GAMA(16), GAMA(17), GAMA(18), 3 GAMA(19), GAMA(20), GAMA(21), GAMA(22)/ 4 6.29960524947436582D-01, 2.51984209978974633D-01, 5 1.54790300415655846D-01, 1.10713062416159013D-01, 6 8.57309395527394825D-02, 6.97161316958684292D-02, 7 5.86085671893713576D-02, 5.04698873536310685D-02, 8 4.42600580689154809D-02, 3.93720661543509966D-02, 9 3.54283195924455368D-02, 3.21818857502098231D-02, A 2.94646240791157679D-02, 2.71581677112934479D-02, B 2.51768272973861779D-02, 2.34570755306078891D-02, C 2.19508390134907203D-02, 2.06210828235646240D-02, D 1.94388240897880846D-02, 1.83810633800683158D-02, E 1.74293213231963172D-02, 1.65685837786612353D-02/ DATA GAMA(23), GAMA(24), GAMA(25), GAMA(26), GAMA(27), GAMA(28), 1 GAMA(29), GAMA(30)/ 2 1.57865285987918445D-02, 1.50729501494095594D-02, 3 1.44193250839954639D-02, 1.38184805735341786D-02, 4 1.32643378994276568D-02, 1.27517121970498651D-02, 5 1.22761545318762767D-02, 1.18338262398482403D-02/ DATA EX1, EX2, HPI, GPI, THPI / 1 3.33333333333333333D-01, 6.66666666666666667D-01, 2 1.57079632679489662D+00, 3.14159265358979324D+00, 3 4.71238898038468986D+00/ DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 / C***FIRST EXECUTABLE STATEMENT ZUNHJ RFNU = 1.0D0/FNU C----------------------------------------------------------------------- C OVERFLOW TEST (Z/FNU TOO SMALL) C----------------------------------------------------------------------- TEST = D1MACH(1)*1.0D+3 AC = FNU*TEST IF (ABS(ZR).GT.AC .OR. ABS(ZI).GT.AC) GO TO 15 ZETA1R = 2.0D0*ABS(LOG(TEST))+FNU ZETA1I = 0.0D0 ZETA2R = FNU ZETA2I = 0.0D0 PHIR = 1.0D0 PHII = 0.0D0 ARGR = 1.0D0 ARGI = 0.0D0 RETURN 15 CONTINUE ZBR = ZR*RFNU ZBI = ZI*RFNU RFNU2 = RFNU*RFNU C----------------------------------------------------------------------- C COMPUTE IN THE FOURTH QUADRANT C----------------------------------------------------------------------- FN13 = FNU**EX1 FN23 = FN13*FN13 RFN13 = 1.0D0/FN13 W2R = CONER - ZBR*ZBR + ZBI*ZBI W2I = CONEI - ZBR*ZBI - ZBR*ZBI AW2 = ZABS(W2R,W2I) IF (AW2.GT.0.25D0) GO TO 130 C----------------------------------------------------------------------- C POWER SERIES FOR ABS(W2).LE.0.25D0 C----------------------------------------------------------------------- K = 1 PR(1) = CONER PI(1) = CONEI SUMAR = GAMA(1) SUMAI = ZEROI AP(1) = 1.0D0 IF (AW2.LT.TOL) GO TO 20 DO 10 K=2,30 PR(K) = PR(K-1)*W2R - PI(K-1)*W2I PI(K) = PR(K-1)*W2I + PI(K-1)*W2R SUMAR = SUMAR + PR(K)*GAMA(K) SUMAI = SUMAI + PI(K)*GAMA(K) AP(K) = AP(K-1)*AW2 IF (AP(K).LT.TOL) GO TO 20 10 CONTINUE K = 30 20 CONTINUE KMAX = K ZETAR = W2R*SUMAR - W2I*SUMAI ZETAI = W2R*SUMAI + W2I*SUMAR ARGR = ZETAR*FN23 ARGI = ZETAI*FN23 CALL ZSQRT(SUMAR, SUMAI, ZAR, ZAI) CALL ZSQRT(W2R, W2I, STR, STI) ZETA2R = STR*FNU ZETA2I = STI*FNU STR = CONER + EX2*(ZETAR*ZAR-ZETAI*ZAI) STI = CONEI + EX2*(ZETAR*ZAI+ZETAI*ZAR) ZETA1R = STR*ZETA2R - STI*ZETA2I ZETA1I = STR*ZETA2I + STI*ZETA2R ZAR = ZAR + ZAR ZAI = ZAI + ZAI CALL ZSQRT(ZAR, ZAI, STR, STI) PHIR = STR*RFN13 PHII = STI*RFN13 IF (IPMTR.EQ.1) GO TO 120 C----------------------------------------------------------------------- C SUM SERIES FOR ASUM AND BSUM C----------------------------------------------------------------------- SUMBR = ZEROR SUMBI = ZEROI DO 30 K=1,KMAX SUMBR = SUMBR + PR(K)*BETA(K) SUMBI = SUMBI + PI(K)*BETA(K) 30 CONTINUE ASUMR = ZEROR ASUMI = ZEROI BSUMR = SUMBR BSUMI = SUMBI L1 = 0 L2 = 30 BTOL = TOL*(ABS(BSUMR)+ABS(BSUMI)) ATOL = TOL PP = 1.0D0 IAS = 0 IBS = 0 IF (RFNU2.LT.TOL) GO TO 110 DO 100 IS=2,7 ATOL = ATOL/RFNU2 PP = PP*RFNU2 IF (IAS.EQ.1) GO TO 60 SUMAR = ZEROR SUMAI = ZEROI DO 40 K=1,KMAX M = L1 + K SUMAR = SUMAR + PR(K)*ALFA(M) SUMAI = SUMAI + PI(K)*ALFA(M) IF (AP(K).LT.ATOL) GO TO 50 40 CONTINUE 50 CONTINUE ASUMR = ASUMR + SUMAR*PP ASUMI = ASUMI + SUMAI*PP IF (PP.LT.TOL) IAS = 1 60 CONTINUE IF (IBS.EQ.1) GO TO 90 SUMBR = ZEROR SUMBI = ZEROI DO 70 K=1,KMAX M = L2 + K SUMBR = SUMBR + PR(K)*BETA(M) SUMBI = SUMBI + PI(K)*BETA(M) IF (AP(K).LT.ATOL) GO TO 80 70 CONTINUE 80 CONTINUE BSUMR = BSUMR + SUMBR*PP BSUMI = BSUMI + SUMBI*PP IF (PP.LT.BTOL) IBS = 1 90 CONTINUE IF (IAS.EQ.1 .AND. IBS.EQ.1) GO TO 110 L1 = L1 + 30 L2 = L2 + 30 100 CONTINUE 110 CONTINUE ASUMR = ASUMR + CONER PP = RFNU*RFN13 BSUMR = BSUMR*PP BSUMI = BSUMI*PP 120 CONTINUE RETURN C----------------------------------------------------------------------- C ABS(W2).GT.0.25D0 C----------------------------------------------------------------------- 130 CONTINUE CALL ZSQRT(W2R, W2I, WR, WI) IF (WR.LT.0.0D0) WR = 0.0D0 IF (WI.LT.0.0D0) WI = 0.0D0 STR = CONER + WR STI = WI CALL ZDIV(STR, STI, ZBR, ZBI, ZAR, ZAI) CALL ZLOG(ZAR, ZAI, ZCR, ZCI, IDUM) IF (ZCI.LT.0.0D0) ZCI = 0.0D0 IF (ZCI.GT.HPI) ZCI = HPI IF (ZCR.LT.0.0D0) ZCR = 0.0D0 ZTHR = (ZCR-WR)*1.5D0 ZTHI = (ZCI-WI)*1.5D0 ZETA1R = ZCR*FNU ZETA1I = ZCI*FNU ZETA2R = WR*FNU ZETA2I = WI*FNU AZTH = ZABS(ZTHR,ZTHI) ANG = THPI IF (ZTHR.GE.0.0D0 .AND. ZTHI.LT.0.0D0) GO TO 140 ANG = HPI IF (ZTHR.EQ.0.0D0) GO TO 140 ANG = DATAN(ZTHI/ZTHR) IF (ZTHR.LT.0.0D0) ANG = ANG + GPI 140 CONTINUE PP = AZTH**EX2 ANG = ANG*EX2 ZETAR = PP*COS(ANG) ZETAI = PP*SIN(ANG) IF (ZETAI.LT.0.0D0) ZETAI = 0.0D0 ARGR = ZETAR*FN23 ARGI = ZETAI*FN23 CALL ZDIV(ZTHR, ZTHI, ZETAR, ZETAI, RTZTR, RTZTI) CALL ZDIV(RTZTR, RTZTI, WR, WI, ZAR, ZAI) TZAR = ZAR + ZAR TZAI = ZAI + ZAI CALL ZSQRT(TZAR, TZAI, STR, STI) PHIR = STR*RFN13 PHII = STI*RFN13 IF (IPMTR.EQ.1) GO TO 120 RAW = 1.0D0/SQRT(AW2) STR = WR*RAW STI = -WI*RAW TFNR = STR*RFNU*RAW TFNI = STI*RFNU*RAW RAZTH = 1.0D0/AZTH STR = ZTHR*RAZTH STI = -ZTHI*RAZTH RZTHR = STR*RAZTH*RFNU RZTHI = STI*RAZTH*RFNU ZCR = RZTHR*AR(2) ZCI = RZTHI*AR(2) RAW2 = 1.0D0/AW2 STR = W2R*RAW2 STI = -W2I*RAW2 T2R = STR*RAW2 T2I = STI*RAW2 STR = T2R*C(2) + C(3) STI = T2I*C(2) UPR(2) = STR*TFNR - STI*TFNI UPI(2) = STR*TFNI + STI*TFNR BSUMR = UPR(2) + ZCR BSUMI = UPI(2) + ZCI ASUMR = ZEROR ASUMI = ZEROI IF (RFNU.LT.TOL) GO TO 220 PRZTHR = RZTHR PRZTHI = RZTHI PTFNR = TFNR PTFNI = TFNI UPR(1) = CONER UPI(1) = CONEI PP = 1.0D0 BTOL = TOL*(ABS(BSUMR)+ABS(BSUMI)) KS = 0 KP1 = 2 L = 3 IAS = 0 IBS = 0 DO 210 LR=2,12,2 LRP1 = LR + 1 C----------------------------------------------------------------------- C COMPUTE TWO ADDITIONAL CR, DR, AND UP FOR TWO MORE TERMS IN C NEXT SUMA AND SUMB C----------------------------------------------------------------------- DO 160 K=LR,LRP1 KS = KS + 1 KP1 = KP1 + 1 L = L + 1 ZAR = C(L) ZAI = ZEROI DO 150 J=2,KP1 L = L + 1 STR = ZAR*T2R - T2I*ZAI + C(L) ZAI = ZAR*T2I + ZAI*T2R ZAR = STR 150 CONTINUE STR = PTFNR*TFNR - PTFNI*TFNI PTFNI = PTFNR*TFNI + PTFNI*TFNR PTFNR = STR UPR(KP1) = PTFNR*ZAR - PTFNI*ZAI UPI(KP1) = PTFNI*ZAR + PTFNR*ZAI CRR(KS) = PRZTHR*BR(KS+1) CRI(KS) = PRZTHI*BR(KS+1) STR = PRZTHR*RZTHR - PRZTHI*RZTHI PRZTHI = PRZTHR*RZTHI + PRZTHI*RZTHR PRZTHR = STR DRR(KS) = PRZTHR*AR(KS+2) DRI(KS) = PRZTHI*AR(KS+2) 160 CONTINUE PP = PP*RFNU2 IF (IAS.EQ.1) GO TO 180 SUMAR = UPR(LRP1) SUMAI = UPI(LRP1) JU = LRP1 DO 170 JR=1,LR JU = JU - 1 SUMAR = SUMAR + CRR(JR)*UPR(JU) - CRI(JR)*UPI(JU) SUMAI = SUMAI + CRR(JR)*UPI(JU) + CRI(JR)*UPR(JU) 170 CONTINUE ASUMR = ASUMR + SUMAR ASUMI = ASUMI + SUMAI TEST = ABS(SUMAR) + ABS(SUMAI) IF (PP.LT.TOL .AND. TEST.LT.TOL) IAS = 1 180 CONTINUE IF (IBS.EQ.1) GO TO 200 SUMBR = UPR(LR+2) + UPR(LRP1)*ZCR - UPI(LRP1)*ZCI SUMBI = UPI(LR+2) + UPR(LRP1)*ZCI + UPI(LRP1)*ZCR JU = LRP1 DO 190 JR=1,LR JU = JU - 1 SUMBR = SUMBR + DRR(JR)*UPR(JU) - DRI(JR)*UPI(JU) SUMBI = SUMBI + DRR(JR)*UPI(JU) + DRI(JR)*UPR(JU) 190 CONTINUE BSUMR = BSUMR + SUMBR BSUMI = BSUMI + SUMBI TEST = ABS(SUMBR) + ABS(SUMBI) IF (PP.LT.BTOL .AND. TEST.LT.BTOL) IBS = 1 200 CONTINUE IF (IAS.EQ.1 .AND. IBS.EQ.1) GO TO 220 210 CONTINUE 220 CONTINUE ASUMR = ASUMR + CONER STR = -BSUMR*RFN13 STI = -BSUMI*RFN13 CALL ZDIV(STR, STI, RTZTR, RTZTI, BSUMR, BSUMI) GO TO 120 END SUBROUTINE ZUNI1 (ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL, + TOL, ELIM, ALIM) C***BEGIN PROLOGUE ZUNI1 C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESI and ZBESK C***LIBRARY SLATEC C***TYPE ALL (CUNI1-A, ZUNI1-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZUNI1 COMPUTES I(FNU,Z) BY MEANS OF THE UNIFORM ASYMPTOTIC C EXPANSION FOR I(FNU,Z) IN -PI/3.LE.ARG Z.LE.PI/3. C C FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC C EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET. C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL. C Y(I)=CZERO FOR I=NLAST+1,N C C***SEE ALSO ZBESI, ZBESK C***ROUTINES CALLED D1MACH, ZABS, ZUCHK, ZUNIK, ZUOIK C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZUNI1 C COMPLEX CFN,CONE,CRSC,CSCL,CSR,CSS,CWRK,CZERO,C1,C2,PHI,RZ,SUM,S1, C *S2,Y,Z,ZETA1,ZETA2 DOUBLE PRECISION ALIM, APHI, ASCLE, BRY, CONER, CRSC, * CSCL, CSRR, CSSR, CWRKI, CWRKR, C1R, C2I, C2M, C2R, ELIM, FN, * FNU, FNUL, PHII, PHIR, RAST, RS1, RZI, RZR, STI, STR, SUMI, * SUMR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZEROI, ZEROR, ZETA1I, * ZETA1R, ZETA2I, ZETA2R, ZI, ZR, CYR, CYI, D1MACH, ZABS INTEGER I, IFLAG, INIT, K, KODE, M, N, ND, NLAST, NN, NUF, NW, NZ DIMENSION BRY(3), YR(N), YI(N), CWRKR(16), CWRKI(16), CSSR(3), * CSRR(3), CYR(2), CYI(2) EXTERNAL ZABS DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 / C***FIRST EXECUTABLE STATEMENT ZUNI1 NZ = 0 ND = N NLAST = 0 C----------------------------------------------------------------------- C COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG- C NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE, C EXP(ALIM)=EXP(ELIM)*TOL C----------------------------------------------------------------------- CSCL = 1.0D0/TOL CRSC = TOL CSSR(1) = CSCL CSSR(2) = CONER CSSR(3) = CRSC CSRR(1) = CRSC CSRR(2) = CONER CSRR(3) = CSCL BRY(1) = 1.0D+3*D1MACH(1)/TOL C----------------------------------------------------------------------- C CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER C----------------------------------------------------------------------- FN = MAX(FNU,1.0D0) INIT = 0 CALL ZUNIK(ZR, ZI, FN, 1, 1, TOL, INIT, PHIR, PHII, ZETA1R, * ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI) IF (KODE.EQ.1) GO TO 10 STR = ZR + ZETA2R STI = ZI + ZETA2I RAST = FN/ZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = -ZETA1R + STR S1I = -ZETA1I + STI GO TO 20 10 CONTINUE S1R = -ZETA1R + ZETA2R S1I = -ZETA1I + ZETA2I 20 CONTINUE RS1 = S1R IF (ABS(RS1).GT.ELIM) GO TO 130 30 CONTINUE NN = MIN(2,ND) DO 80 I=1,NN FN = FNU + (ND-I) INIT = 0 CALL ZUNIK(ZR, ZI, FN, 1, 0, TOL, INIT, PHIR, PHII, ZETA1R, * ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI) IF (KODE.EQ.1) GO TO 40 STR = ZR + ZETA2R STI = ZI + ZETA2I RAST = FN/ZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = -ZETA1R + STR S1I = -ZETA1I + STI + ZI GO TO 50 40 CONTINUE S1R = -ZETA1R + ZETA2R S1I = -ZETA1I + ZETA2I 50 CONTINUE C----------------------------------------------------------------------- C TEST FOR UNDERFLOW AND OVERFLOW C----------------------------------------------------------------------- RS1 = S1R IF (ABS(RS1).GT.ELIM) GO TO 110 IF (I.EQ.1) IFLAG = 2 IF (ABS(RS1).LT.ALIM) GO TO 60 C----------------------------------------------------------------------- C REFINE TEST AND SCALE C----------------------------------------------------------------------- APHI = ZABS(PHIR,PHII) RS1 = RS1 + LOG(APHI) IF (ABS(RS1).GT.ELIM) GO TO 110 IF (I.EQ.1) IFLAG = 1 IF (RS1.LT.0.0D0) GO TO 60 IF (I.EQ.1) IFLAG = 3 60 CONTINUE C----------------------------------------------------------------------- C SCALE S1 IF ABS(S1).LT.ASCLE C----------------------------------------------------------------------- S2R = PHIR*SUMR - PHII*SUMI S2I = PHIR*SUMI + PHII*SUMR STR = EXP(S1R)*CSSR(IFLAG) S1R = STR*COS(S1I) S1I = STR*SIN(S1I) STR = S2R*S1R - S2I*S1I S2I = S2R*S1I + S2I*S1R S2R = STR IF (IFLAG.NE.1) GO TO 70 CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL) IF (NW.NE.0) GO TO 110 70 CONTINUE CYR(I) = S2R CYI(I) = S2I M = ND - I + 1 YR(M) = S2R*CSRR(IFLAG) YI(M) = S2I*CSRR(IFLAG) 80 CONTINUE IF (ND.LE.2) GO TO 100 RAST = 1.0D0/ZABS(ZR,ZI) STR = ZR*RAST STI = -ZI*RAST RZR = (STR+STR)*RAST RZI = (STI+STI)*RAST BRY(2) = 1.0D0/BRY(1) BRY(3) = D1MACH(2) S1R = CYR(1) S1I = CYI(1) S2R = CYR(2) S2I = CYI(2) C1R = CSRR(IFLAG) ASCLE = BRY(IFLAG) K = ND - 2 FN = K DO 90 I=3,ND C2R = S2R C2I = S2I S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I) S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R) S1R = C2R S1I = C2I C2R = S2R*C1R C2I = S2I*C1R YR(K) = C2R YI(K) = C2I K = K - 1 FN = FN - 1.0D0 IF (IFLAG.GE.3) GO TO 90 STR = ABS(C2R) STI = ABS(C2I) C2M = MAX(STR,STI) IF (C2M.LE.ASCLE) GO TO 90 IFLAG = IFLAG + 1 ASCLE = BRY(IFLAG) S1R = S1R*C1R S1I = S1I*C1R S2R = C2R S2I = C2I S1R = S1R*CSSR(IFLAG) S1I = S1I*CSSR(IFLAG) S2R = S2R*CSSR(IFLAG) S2I = S2I*CSSR(IFLAG) C1R = CSRR(IFLAG) 90 CONTINUE 100 CONTINUE RETURN C----------------------------------------------------------------------- C SET UNDERFLOW AND UPDATE PARAMETERS C----------------------------------------------------------------------- 110 CONTINUE IF (RS1.GT.0.0D0) GO TO 120 YR(ND) = ZEROR YI(ND) = ZEROI NZ = NZ + 1 ND = ND - 1 IF (ND.EQ.0) GO TO 100 CALL ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, NUF, TOL, ELIM, ALIM) IF (NUF.LT.0) GO TO 120 ND = ND - NUF NZ = NZ + NUF IF (ND.EQ.0) GO TO 100 FN = FNU + (ND-1) IF (FN.GE.FNUL) GO TO 30 NLAST = ND RETURN 120 CONTINUE NZ = -1 RETURN 130 CONTINUE IF (RS1.GT.0.0D0) GO TO 120 NZ = N DO 140 I=1,N YR(I) = ZEROR YI(I) = ZEROI 140 CONTINUE RETURN END SUBROUTINE ZUNI2 (ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL, + TOL, ELIM, ALIM) C***BEGIN PROLOGUE ZUNI2 C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESI and ZBESK C***LIBRARY SLATEC C***TYPE ALL (CUNI2-A, ZUNI2-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF C UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I C OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO. C C FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC C EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET. C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL. C Y(I)=CZERO FOR I=NLAST+1,N C C***SEE ALSO ZBESI, ZBESK C***ROUTINES CALLED D1MACH, ZABS, ZAIRY, ZUCHK, ZUNHJ, ZUOIK C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZUNI2 C COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS, C *CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGI, * ARGR, ASCLE, ASUMI, ASUMR, BRY, BSUMI, BSUMR, CIDI, CIPI, CIPR, * CONER, CRSC, CSCL, CSRR, CSSR, C1R, C2I, C2M, C2R, DAII, * DAIR, ELIM, FN, FNU, FNUL, HPI, PHII, PHIR, RAST, RAZ, RS1, RZI, * RZR, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZBI, ZBR, ZEROI, * ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI, ZNI, ZNR, ZR, CYR, * CYI, D1MACH, ZABS, CAR, SAR INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST, * NN, NUF, NW, NZ, IDUM DIMENSION BRY(3), YR(N), YI(N), CIPR(4), CIPI(4), CSSR(3), * CSRR(3), CYR(2), CYI(2) EXTERNAL ZABS DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 / DATA CIPR(1),CIPI(1),CIPR(2),CIPI(2),CIPR(3),CIPI(3),CIPR(4), * CIPI(4)/ 1.0D0,0.0D0, 0.0D0,1.0D0, -1.0D0,0.0D0, 0.0D0,-1.0D0/ DATA HPI, AIC / 1 1.57079632679489662D+00, 1.265512123484645396D+00/ C***FIRST EXECUTABLE STATEMENT ZUNI2 NZ = 0 ND = N NLAST = 0 C----------------------------------------------------------------------- C COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG- C NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE, C EXP(ALIM)=EXP(ELIM)*TOL C----------------------------------------------------------------------- CSCL = 1.0D0/TOL CRSC = TOL CSSR(1) = CSCL CSSR(2) = CONER CSSR(3) = CRSC CSRR(1) = CRSC CSRR(2) = CONER CSRR(3) = CSCL BRY(1) = 1.0D+3*D1MACH(1)/TOL C----------------------------------------------------------------------- C ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI C----------------------------------------------------------------------- ZNR = ZI ZNI = -ZR ZBR = ZR ZBI = ZI CIDI = -CONER INU = FNU ANG = HPI*(FNU-INU) C2R = COS(ANG) C2I = SIN(ANG) CAR = C2R SAR = C2I IN = INU + N - 1 IN = MOD(IN,4) + 1 STR = C2R*CIPR(IN) - C2I*CIPI(IN) C2I = C2R*CIPI(IN) + C2I*CIPR(IN) C2R = STR IF (ZI.GT.0.0D0) GO TO 10 ZNR = -ZNR ZBI = -ZBI CIDI = -CIDI C2I = -C2I 10 CONTINUE C----------------------------------------------------------------------- C CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER C----------------------------------------------------------------------- FN = MAX(FNU,1.0D0) CALL ZUNHJ(ZNR, ZNI, FN, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) IF (KODE.EQ.1) GO TO 20 STR = ZBR + ZETA2R STI = ZBI + ZETA2I RAST = FN/ZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = -ZETA1R + STR S1I = -ZETA1I + STI GO TO 30 20 CONTINUE S1R = -ZETA1R + ZETA2R S1I = -ZETA1I + ZETA2I 30 CONTINUE RS1 = S1R IF (ABS(RS1).GT.ELIM) GO TO 150 40 CONTINUE NN = MIN(2,ND) DO 90 I=1,NN FN = FNU + (ND-I) CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIR, PHII, ARGR, ARGI, * ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) IF (KODE.EQ.1) GO TO 50 STR = ZBR + ZETA2R STI = ZBI + ZETA2I RAST = FN/ZABS(STR,STI) STR = STR*RAST*RAST STI = -STI*RAST*RAST S1R = -ZETA1R + STR S1I = -ZETA1I + STI + ABS(ZI) GO TO 60 50 CONTINUE S1R = -ZETA1R + ZETA2R S1I = -ZETA1I + ZETA2I 60 CONTINUE C----------------------------------------------------------------------- C TEST FOR UNDERFLOW AND OVERFLOW C----------------------------------------------------------------------- RS1 = S1R IF (ABS(RS1).GT.ELIM) GO TO 120 IF (I.EQ.1) IFLAG = 2 IF (ABS(RS1).LT.ALIM) GO TO 70 C----------------------------------------------------------------------- C REFINE TEST AND SCALE C----------------------------------------------------------------------- C----------------------------------------------------------------------- APHI = ZABS(PHIR,PHII) AARG = ZABS(ARGR,ARGI) RS1 = RS1 + LOG(APHI) - 0.25D0*LOG(AARG) - AIC IF (ABS(RS1).GT.ELIM) GO TO 120 IF (I.EQ.1) IFLAG = 1 IF (RS1.LT.0.0D0) GO TO 70 IF (I.EQ.1) IFLAG = 3 70 CONTINUE C----------------------------------------------------------------------- C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR C EXPONENT EXTREMES C----------------------------------------------------------------------- CALL ZAIRY(ARGR, ARGI, 0, 2, AIR, AII, NAI, IDUM) CALL ZAIRY(ARGR, ARGI, 1, 2, DAIR, DAII, NDAI, IDUM) STR = DAIR*BSUMR - DAII*BSUMI STI = DAIR*BSUMI + DAII*BSUMR STR = STR + (AIR*ASUMR-AII*ASUMI) STI = STI + (AIR*ASUMI+AII*ASUMR) S2R = PHIR*STR - PHII*STI S2I = PHIR*STI + PHII*STR STR = EXP(S1R)*CSSR(IFLAG) S1R = STR*COS(S1I) S1I = STR*SIN(S1I) STR = S2R*S1R - S2I*S1I S2I = S2R*S1I + S2I*S1R S2R = STR IF (IFLAG.NE.1) GO TO 80 CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL) IF (NW.NE.0) GO TO 120 80 CONTINUE IF (ZI.LE.0.0D0) S2I = -S2I STR = S2R*C2R - S2I*C2I S2I = S2R*C2I + S2I*C2R S2R = STR CYR(I) = S2R CYI(I) = S2I J = ND - I + 1 YR(J) = S2R*CSRR(IFLAG) YI(J) = S2I*CSRR(IFLAG) STR = -C2I*CIDI C2I = C2R*CIDI C2R = STR 90 CONTINUE IF (ND.LE.2) GO TO 110 RAZ = 1.0D0/ZABS(ZR,ZI) STR = ZR*RAZ STI = -ZI*RAZ RZR = (STR+STR)*RAZ RZI = (STI+STI)*RAZ BRY(2) = 1.0D0/BRY(1) BRY(3) = D1MACH(2) S1R = CYR(1) S1I = CYI(1) S2R = CYR(2) S2I = CYI(2) C1R = CSRR(IFLAG) ASCLE = BRY(IFLAG) K = ND - 2 FN = K DO 100 I=3,ND C2R = S2R C2I = S2I S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I) S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R) S1R = C2R S1I = C2I C2R = S2R*C1R C2I = S2I*C1R YR(K) = C2R YI(K) = C2I K = K - 1 FN = FN - 1.0D0 IF (IFLAG.GE.3) GO TO 100 STR = ABS(C2R) STI = ABS(C2I) C2M = MAX(STR,STI) IF (C2M.LE.ASCLE) GO TO 100 IFLAG = IFLAG + 1 ASCLE = BRY(IFLAG) S1R = S1R*C1R S1I = S1I*C1R S2R = C2R S2I = C2I S1R = S1R*CSSR(IFLAG) S1I = S1I*CSSR(IFLAG) S2R = S2R*CSSR(IFLAG) S2I = S2I*CSSR(IFLAG) C1R = CSRR(IFLAG) 100 CONTINUE 110 CONTINUE RETURN 120 CONTINUE IF (RS1.GT.0.0D0) GO TO 140 C----------------------------------------------------------------------- C SET UNDERFLOW AND UPDATE PARAMETERS C----------------------------------------------------------------------- YR(ND) = ZEROR YI(ND) = ZEROI NZ = NZ + 1 ND = ND - 1 IF (ND.EQ.0) GO TO 110 CALL ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, NUF, TOL, ELIM, ALIM) IF (NUF.LT.0) GO TO 140 ND = ND - NUF NZ = NZ + NUF IF (ND.EQ.0) GO TO 110 FN = FNU + (ND-1) IF (FN.LT.FNUL) GO TO 130 C FN = CIDI C J = NUF + 1 C K = MOD(J,4) + 1 C S1R = CIPR(K) C S1I = CIPI(K) C IF (FN.LT.0.0D0) S1I = -S1I C STR = C2R*S1R - C2I*S1I C C2I = C2R*S1I + C2I*S1R C C2R = STR IN = INU + ND - 1 IN = MOD(IN,4) + 1 C2R = CAR*CIPR(IN) - SAR*CIPI(IN) C2I = CAR*CIPI(IN) + SAR*CIPR(IN) IF (ZI.LE.0.0D0) C2I = -C2I GO TO 40 130 CONTINUE NLAST = ND RETURN 140 CONTINUE NZ = -1 RETURN 150 CONTINUE IF (RS1.GT.0.0D0) GO TO 140 NZ = N DO 160 I=1,N YR(I) = ZEROR YI(I) = ZEROI 160 CONTINUE RETURN END SUBROUTINE ZUNIK (ZRR, ZRI, FNU, IKFLG, IPMTR, TOL, INIT, PHIR, + PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI) C***BEGIN PROLOGUE ZUNIK C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESI and ZBESK C***LIBRARY SLATEC C***TYPE ALL (CUNIK-A, ZUNIK-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZUNIK COMPUTES PARAMETERS FOR THE UNIFORM ASYMPTOTIC C EXPANSIONS OF THE I AND K FUNCTIONS ON IKFLG= 1 OR 2 C RESPECTIVELY BY C C W(FNU,ZR) = PHI*EXP(ZETA)*SUM C C WHERE ZETA=-ZETA1 + ZETA2 OR C ZETA1 - ZETA2 C C THE FIRST CALL MUST HAVE INIT=0. SUBSEQUENT CALLS WITH THE C SAME ZR AND FNU WILL RETURN THE I OR K FUNCTION ON IKFLG= C 1 OR 2 WITH NO CHANGE IN INIT. CWRK IS A COMPLEX WORK C ARRAY. IPMTR=0 COMPUTES ALL PARAMETERS. IPMTR=1 COMPUTES PHI, C ZETA1,ZETA2. C C***SEE ALSO ZBESI, ZBESK C***ROUTINES CALLED D1MACH, ZDIV, ZLOG, ZSQRT C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZUNIK C COMPLEX CFN,CON,CONE,CRFN,CWRK,CZERO,PHI,S,SR,SUM,T,T2,ZETA1, C *ZETA2,ZN,ZR DOUBLE PRECISION AC, C, CON, CONEI, CONER, CRFNI, CRFNR, CWRKI, * CWRKR, FNU, PHII, PHIR, RFN, SI, SR, SRI, SRR, STI, STR, SUMI, * SUMR, TEST, TI, TOL, TR, T2I, T2R, ZEROI, ZEROR, ZETA1I, ZETA1R, * ZETA2I, ZETA2R, ZNI, ZNR, ZRI, ZRR, D1MACH INTEGER I, IDUM, IKFLG, INIT, IPMTR, J, K, L DIMENSION C(120), CWRKR(16), CWRKI(16), CON(2) DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 / DATA CON(1), CON(2) / 1 3.98942280401432678D-01, 1.25331413731550025D+00 / 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 1.00000000000000000D+00, -2.08333333333333333D-01, 4 1.25000000000000000D-01, 3.34201388888888889D-01, 5 -4.01041666666666667D-01, 7.03125000000000000D-02, 6 -1.02581259645061728D+00, 1.84646267361111111D+00, 7 -8.91210937500000000D-01, 7.32421875000000000D-02, 8 4.66958442342624743D+00, -1.12070026162229938D+01, 9 8.78912353515625000D+00, -2.36408691406250000D+00, A 1.12152099609375000D-01, -2.82120725582002449D+01, B 8.46362176746007346D+01, -9.18182415432400174D+01, C 4.25349987453884549D+01, -7.36879435947963170D+00, D 2.27108001708984375D-01, 2.12570130039217123D+02, E -7.65252468141181642D+02, 1.05999045252799988D+03/ 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 -6.99579627376132541D+02, 2.18190511744211590D+02, 4 -2.64914304869515555D+01, 5.72501420974731445D-01, 5 -1.91945766231840700D+03, 8.06172218173730938D+03, 6 -1.35865500064341374D+04, 1.16553933368645332D+04, 7 -5.30564697861340311D+03, 1.20090291321635246D+03, 8 -1.08090919788394656D+02, 1.72772750258445740D+00, 9 2.02042913309661486D+04, -9.69805983886375135D+04, A 1.92547001232531532D+05, -2.03400177280415534D+05, B 1.22200464983017460D+05, -4.11926549688975513D+04, C 7.10951430248936372D+03, -4.93915304773088012D+02, D 6.07404200127348304D+00, -2.42919187900551333D+05, E 1.31176361466297720D+06, -2.99801591853810675D+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), C(66), C(67), C(68), C(69), C(70), C(71), C(72)/ 3 3.76327129765640400D+06, -2.81356322658653411D+06, 4 1.26836527332162478D+06, -3.31645172484563578D+05, 5 4.52187689813627263D+04, -2.49983048181120962D+03, 6 2.43805296995560639D+01, 3.28446985307203782D+06, 7 -1.97068191184322269D+07, 5.09526024926646422D+07, 8 -7.41051482115326577D+07, 6.63445122747290267D+07, 9 -3.75671766607633513D+07, 1.32887671664218183D+07, A -2.78561812808645469D+06, 3.08186404612662398D+05, B -1.38860897537170405D+04, 1.10017140269246738D+02, C -4.93292536645099620D+07, 3.25573074185765749D+08, D -9.39462359681578403D+08, 1.55359689957058006D+09, E -1.62108055210833708D+09, 1.10684281682301447D+09/ DATA C(73), C(74), C(75), C(76), C(77), C(78), C(79), C(80), 1 C(81), C(82), C(83), C(84), C(85), C(86), C(87), C(88), 2 C(89), C(90), C(91), C(92), C(93), C(94), C(95), C(96)/ 3 -4.95889784275030309D+08, 1.42062907797533095D+08, 4 -2.44740627257387285D+07, 2.24376817792244943D+06, 5 -8.40054336030240853D+04, 5.51335896122020586D+02, 6 8.14789096118312115D+08, -5.86648149205184723D+09, 7 1.86882075092958249D+10, -3.46320433881587779D+10, 8 4.12801855797539740D+10, -3.30265997498007231D+10, 9 1.79542137311556001D+10, -6.56329379261928433D+09, A 1.55927986487925751D+09, -2.25105661889415278D+08, B 1.73951075539781645D+07, -5.49842327572288687D+05, C 3.03809051092238427D+03, -1.46792612476956167D+10, D 1.14498237732025810D+11, -3.99096175224466498D+11, E 8.19218669548577329D+11, -1.09837515608122331D+12/ DATA C(97), C(98), C(99), C(100), C(101), C(102), C(103), C(104), 1 C(105), C(106), C(107), C(108), C(109), C(110), C(111), 2 C(112), C(113), C(114), C(115), C(116), C(117), C(118)/ 3 1.00815810686538209D+12, -6.45364869245376503D+11, 4 2.87900649906150589D+11, -8.78670721780232657D+10, 5 1.76347306068349694D+10, -2.16716498322379509D+09, 6 1.43157876718888981D+08, -3.87183344257261262D+06, 7 1.82577554742931747D+04, 2.86464035717679043D+11, 8 -2.40629790002850396D+12, 9.10934118523989896D+12, 9 -2.05168994109344374D+13, 3.05651255199353206D+13, A -3.16670885847851584D+13, 2.33483640445818409D+13, B -1.23204913055982872D+13, 4.61272578084913197D+12, C -1.19655288019618160D+12, 2.05914503232410016D+11, D -2.18229277575292237D+10, 1.24700929351271032D+09/ DATA C(119), C(120)/ 1 -2.91883881222208134D+07, 1.18838426256783253D+05/ C***FIRST EXECUTABLE STATEMENT ZUNIK IF (INIT.NE.0) GO TO 40 C----------------------------------------------------------------------- C INITIALIZE ALL VARIABLES C----------------------------------------------------------------------- RFN = 1.0D0/FNU C----------------------------------------------------------------------- C OVERFLOW TEST (ZR/FNU TOO SMALL) C----------------------------------------------------------------------- TEST = D1MACH(1)*1.0D+3 AC = FNU*TEST IF (ABS(ZRR).GT.AC .OR. ABS(ZRI).GT.AC) GO TO 15 ZETA1R = 2.0D0*ABS(LOG(TEST))+FNU ZETA1I = 0.0D0 ZETA2R = FNU ZETA2I = 0.0D0 PHIR = 1.0D0 PHII = 0.0D0 RETURN 15 CONTINUE TR = ZRR*RFN TI = ZRI*RFN SR = CONER + (TR*TR-TI*TI) SI = CONEI + (TR*TI+TI*TR) CALL ZSQRT(SR, SI, SRR, SRI) STR = CONER + SRR STI = CONEI + SRI CALL ZDIV(STR, STI, TR, TI, ZNR, ZNI) CALL ZLOG(ZNR, ZNI, STR, STI, IDUM) ZETA1R = FNU*STR ZETA1I = FNU*STI ZETA2R = FNU*SRR ZETA2I = FNU*SRI CALL ZDIV(CONER, CONEI, SRR, SRI, TR, TI) SRR = TR*RFN SRI = TI*RFN CALL ZSQRT(SRR, SRI, CWRKR(16), CWRKI(16)) PHIR = CWRKR(16)*CON(IKFLG) PHII = CWRKI(16)*CON(IKFLG) IF (IPMTR.NE.0) RETURN CALL ZDIV(CONER, CONEI, SR, SI, T2R, T2I) CWRKR(1) = CONER CWRKI(1) = CONEI CRFNR = CONER CRFNI = CONEI AC = 1.0D0 L = 1 DO 20 K=2,15 SR = ZEROR SI = ZEROI DO 10 J=1,K L = L + 1 STR = SR*T2R - SI*T2I + C(L) SI = SR*T2I + SI*T2R SR = STR 10 CONTINUE STR = CRFNR*SRR - CRFNI*SRI CRFNI = CRFNR*SRI + CRFNI*SRR CRFNR = STR CWRKR(K) = CRFNR*SR - CRFNI*SI CWRKI(K) = CRFNR*SI + CRFNI*SR AC = AC*RFN TEST = ABS(CWRKR(K)) + ABS(CWRKI(K)) IF (AC.LT.TOL .AND. TEST.LT.TOL) GO TO 30 20 CONTINUE K = 15 30 CONTINUE INIT = K 40 CONTINUE IF (IKFLG.EQ.2) GO TO 60 C----------------------------------------------------------------------- C COMPUTE SUM FOR THE I FUNCTION C----------------------------------------------------------------------- SR = ZEROR SI = ZEROI DO 50 I=1,INIT SR = SR + CWRKR(I) SI = SI + CWRKI(I) 50 CONTINUE SUMR = SR SUMI = SI PHIR = CWRKR(16)*CON(1) PHII = CWRKI(16)*CON(1) RETURN 60 CONTINUE C----------------------------------------------------------------------- C COMPUTE SUM FOR THE K FUNCTION C----------------------------------------------------------------------- SR = ZEROR SI = ZEROI TR = CONER DO 70 I=1,INIT SR = SR + TR*CWRKR(I) SI = SI + TR*CWRKI(I) TR = -TR 70 CONTINUE SUMR = SR SUMI = SI PHIR = CWRKR(16)*CON(2) PHII = CWRKI(16)*CON(2) RETURN END SUBROUTINE ZUOIK (ZR, ZI, FNU, KODE, IKFLG, N, YR, YI, NUF, TOL, + ELIM, ALIM) C***BEGIN PROLOGUE ZUOIK C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESH, ZBESI and ZBESK C***LIBRARY SLATEC C***TYPE ALL (CUOIK-A, ZUOIK-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC C EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM C (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW C WHERE ALIM.LT.ELIM. IF THE MAGNITUDE, BASED ON THE LEADING C EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN C THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER C MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE C EXP(-ELIM)=SMALLEST MACHINE NUMBER*1.0E+3 AND EXP(-ALIM)= C EXP(-ELIM)/TOL C C IKFLG=1 MEANS THE I SEQUENCE IS TESTED C =2 MEANS THE K SEQUENCE IS TESTED C NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE C =-1 MEANS AN OVERFLOW WOULD OCCUR C IKFLG=1 AND NUF.GT.0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO C THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE C IKFLG=2 AND NUF.EQ.N MEANS ALL Y VALUES WERE SET TO ZERO C IKFLG=2 AND 0.LT.NUF.LT.N NOT CONSIDERED. Y MUST BE SET BY C ANOTHER ROUTINE C C***SEE ALSO ZBESH, ZBESI, ZBESK C***ROUTINES CALLED D1MACH, ZABS, ZLOG, ZUCHK, ZUNHJ, ZUNIK C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZUOIK C COMPLEX ARG,ASUM,BSUM,CWRK,CZ,CZERO,PHI,SUM,Y,Z,ZB,ZETA1,ZETA2,ZN, C *ZR DOUBLE PRECISION AARG, AIC, ALIM, APHI, ARGI, ARGR, ASUMI, ASUMR, * ASCLE, AX, AY, BSUMI, BSUMR, CWRKI, CWRKR, CZI, CZR, ELIM, FNN, * FNU, GNN, GNU, PHII, PHIR, RCZ, STR, STI, SUMI, SUMR, TOL, YI, * YR, ZBI, ZBR, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI, * ZNI, ZNR, ZR, ZRI, ZRR, D1MACH, ZABS INTEGER I, IDUM, IFORM, IKFLG, INIT, KODE, N, NN, NUF, NW DIMENSION YR(N), YI(N), CWRKR(16), CWRKI(16) EXTERNAL ZABS DATA ZEROR,ZEROI / 0.0D0, 0.0D0 / DATA AIC / 1.265512123484645396D+00 / C***FIRST EXECUTABLE STATEMENT ZUOIK NUF = 0 NN = N ZRR = ZR ZRI = ZI IF (ZR.GE.0.0D0) GO TO 10 ZRR = -ZR ZRI = -ZI 10 CONTINUE ZBR = ZRR ZBI = ZRI AX = ABS(ZR)*1.7321D0 AY = ABS(ZI) IFORM = 1 IF (AY.GT.AX) IFORM = 2 GNU = MAX(FNU,1.0D0) IF (IKFLG.EQ.1) GO TO 20 FNN = NN GNN = FNU + FNN - 1.0D0 GNU = MAX(GNN,FNN) 20 CONTINUE C----------------------------------------------------------------------- C ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE C REAL PARTS OF ZETA1, ZETA2 AND ZB. NO ATTEMPT IS MADE TO GET C THE SIGN OF THE IMAGINARY PART CORRECT. C----------------------------------------------------------------------- IF (IFORM.EQ.2) GO TO 30 INIT = 0 CALL ZUNIK(ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, PHIR, PHII, * ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI) CZR = -ZETA1R + ZETA2R CZI = -ZETA1I + ZETA2I GO TO 50 30 CONTINUE ZNR = ZRI ZNI = -ZRR IF (ZI.GT.0.0D0) GO TO 40 ZNR = -ZNR 40 CONTINUE CALL ZUNHJ(ZNR, ZNI, GNU, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) CZR = -ZETA1R + ZETA2R CZI = -ZETA1I + ZETA2I AARG = ZABS(ARGR,ARGI) 50 CONTINUE IF (KODE.EQ.1) GO TO 60 CZR = CZR - ZBR CZI = CZI - ZBI 60 CONTINUE IF (IKFLG.EQ.1) GO TO 70 CZR = -CZR CZI = -CZI 70 CONTINUE APHI = ZABS(PHIR,PHII) RCZ = CZR C----------------------------------------------------------------------- C OVERFLOW TEST C----------------------------------------------------------------------- IF (RCZ.GT.ELIM) GO TO 210 IF (RCZ.LT.ALIM) GO TO 80 RCZ = RCZ + LOG(APHI) IF (IFORM.EQ.2) RCZ = RCZ - 0.25D0*LOG(AARG) - AIC IF (RCZ.GT.ELIM) GO TO 210 GO TO 130 80 CONTINUE C----------------------------------------------------------------------- C UNDERFLOW TEST C----------------------------------------------------------------------- IF (RCZ.LT.(-ELIM)) GO TO 90 IF (RCZ.GT.(-ALIM)) GO TO 130 RCZ = RCZ + LOG(APHI) IF (IFORM.EQ.2) RCZ = RCZ - 0.25D0*LOG(AARG) - AIC IF (RCZ.GT.(-ELIM)) GO TO 110 90 CONTINUE DO 100 I=1,NN YR(I) = ZEROR YI(I) = ZEROI 100 CONTINUE NUF = NN RETURN 110 CONTINUE ASCLE = 1.0D+3*D1MACH(1)/TOL CALL ZLOG(PHIR, PHII, STR, STI, IDUM) CZR = CZR + STR CZI = CZI + STI IF (IFORM.EQ.1) GO TO 120 CALL ZLOG(ARGR, ARGI, STR, STI, IDUM) CZR = CZR - 0.25D0*STR - AIC CZI = CZI - 0.25D0*STI 120 CONTINUE AX = EXP(RCZ)/TOL AY = CZI CZR = AX*COS(AY) CZI = AX*SIN(AY) CALL ZUCHK(CZR, CZI, NW, ASCLE, TOL) IF (NW.NE.0) GO TO 90 130 CONTINUE IF (IKFLG.EQ.2) RETURN IF (N.EQ.1) RETURN C----------------------------------------------------------------------- C SET UNDERFLOWS ON I SEQUENCE C----------------------------------------------------------------------- 140 CONTINUE GNU = FNU + (NN-1) IF (IFORM.EQ.2) GO TO 150 INIT = 0 CALL ZUNIK(ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, PHIR, PHII, * ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI) CZR = -ZETA1R + ZETA2R CZI = -ZETA1I + ZETA2I GO TO 160 150 CONTINUE CALL ZUNHJ(ZNR, ZNI, GNU, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) CZR = -ZETA1R + ZETA2R CZI = -ZETA1I + ZETA2I AARG = ZABS(ARGR,ARGI) 160 CONTINUE IF (KODE.EQ.1) GO TO 170 CZR = CZR - ZBR CZI = CZI - ZBI 170 CONTINUE APHI = ZABS(PHIR,PHII) RCZ = CZR IF (RCZ.LT.(-ELIM)) GO TO 180 IF (RCZ.GT.(-ALIM)) RETURN RCZ = RCZ + LOG(APHI) IF (IFORM.EQ.2) RCZ = RCZ - 0.25D0*LOG(AARG) - AIC IF (RCZ.GT.(-ELIM)) GO TO 190 180 CONTINUE YR(NN) = ZEROR YI(NN) = ZEROI NN = NN - 1 NUF = NUF + 1 IF (NN.EQ.0) RETURN GO TO 140 190 CONTINUE ASCLE = 1.0D+3*D1MACH(1)/TOL CALL ZLOG(PHIR, PHII, STR, STI, IDUM) CZR = CZR + STR CZI = CZI + STI IF (IFORM.EQ.1) GO TO 200 CALL ZLOG(ARGR, ARGI, STR, STI, IDUM) CZR = CZR - 0.25D0*STR - AIC CZI = CZI - 0.25D0*STI 200 CONTINUE AX = EXP(RCZ)/TOL AY = CZI CZR = AX*COS(AY) CZI = AX*SIN(AY) CALL ZUCHK(CZR, CZI, NW, ASCLE, TOL) IF (NW.NE.0) GO TO 180 RETURN 210 CONTINUE NUF = -1 RETURN END SUBROUTINE ZWRSK (ZRR, ZRI, FNU, KODE, N, YR, YI, NZ, CWR, CWI, + TOL, ELIM, ALIM) C***BEGIN PROLOGUE ZWRSK C***SUBSIDIARY C***PURPOSE Subsidiary to ZBESI and ZBESK C***LIBRARY SLATEC C***TYPE ALL (CWRSK-A, ZWRSK-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C ZWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY C NORMALIZING THE I FUNCTION RATIOS FROM ZRATI BY THE WRONSKIAN C C***SEE ALSO ZBESI, ZBESK C***ROUTINES CALLED D1MACH, ZABS, ZBKNU, ZRATI C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE ZWRSK C COMPLEX CINU,CSCL,CT,CW,C1,C2,RCT,ST,Y,ZR DOUBLE PRECISION ACT, ACW, ALIM, ASCLE, CINUI, CINUR, CSCLR, CTI, * CTR, CWI, CWR, C1I, C1R, C2I, C2R, ELIM, FNU, PTI, PTR, RACT, * STI, STR, TOL, YI, YR, ZRI, ZRR, ZABS, D1MACH INTEGER I, KODE, N, NW, NZ DIMENSION YR(N), YI(N), CWR(2), CWI(2) EXTERNAL ZABS C***FIRST EXECUTABLE STATEMENT ZWRSK C----------------------------------------------------------------------- C I(FNU+I-1,Z) BY BACKWARD RECURRENCE FOR RATIOS C Y(I)=I(FNU+I,Z)/I(FNU+I-1,Z) FROM CRATI NORMALIZED BY THE C WRONSKIAN WITH K(FNU,Z) AND K(FNU+1,Z) FROM CBKNU. C----------------------------------------------------------------------- C NZ = 0 CALL ZBKNU(ZRR, ZRI, FNU, KODE, 2, CWR, CWI, NW, TOL, ELIM, ALIM) IF (NW.NE.0) GO TO 50 CALL ZRATI(ZRR, ZRI, FNU, N, YR, YI, TOL) C----------------------------------------------------------------------- C RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z), C R(FNU+J-1,Z)=Y(J), J=1,...,N C----------------------------------------------------------------------- CINUR = 1.0D0 CINUI = 0.0D0 IF (KODE.EQ.1) GO TO 10 CINUR = COS(ZRI) CINUI = SIN(ZRI) 10 CONTINUE C----------------------------------------------------------------------- C ON LOW EXPONENT MACHINES THE K FUNCTIONS CAN BE CLOSE TO BOTH C THE UNDER AND OVERFLOW LIMITS AND THE NORMALIZATION MUST BE C SCALED TO PREVENT OVER OR UNDERFLOW. CUOIK HAS DETERMINED THAT C THE RESULT IS ON SCALE. C----------------------------------------------------------------------- ACW = ZABS(CWR(2),CWI(2)) ASCLE = 1.0D+3*D1MACH(1)/TOL CSCLR = 1.0D0 IF (ACW.GT.ASCLE) GO TO 20 CSCLR = 1.0D0/TOL GO TO 30 20 CONTINUE ASCLE = 1.0D0/ASCLE IF (ACW.LT.ASCLE) GO TO 30 CSCLR = TOL 30 CONTINUE C1R = CWR(1)*CSCLR C1I = CWI(1)*CSCLR C2R = CWR(2)*CSCLR C2I = CWI(2)*CSCLR STR = YR(1) STI = YI(1) C----------------------------------------------------------------------- C CINU=CINU*(CONJG(CT)/ABS(CT))*(1.0D0/ABS(CT) PREVENTS C UNDER- OR OVERFLOW PREMATURELY BY SQUARING ABS(CT) C----------------------------------------------------------------------- PTR = STR*C1R - STI*C1I PTI = STR*C1I + STI*C1R PTR = PTR + C2R PTI = PTI + C2I CTR = ZRR*PTR - ZRI*PTI CTI = ZRR*PTI + ZRI*PTR ACT = ZABS(CTR,CTI) RACT = 1.0D0/ACT CTR = CTR*RACT CTI = -CTI*RACT PTR = CINUR*RACT PTI = CINUI*RACT CINUR = PTR*CTR - PTI*CTI CINUI = PTR*CTI + PTI*CTR YR(1) = CINUR*CSCLR YI(1) = CINUI*CSCLR IF (N.EQ.1) RETURN DO 40 I=2,N PTR = STR*CINUR - STI*CINUI CINUI = STR*CINUI + STI*CINUR CINUR = PTR STR = YR(I) STI = YI(I) YR(I) = CINUR*CSCLR YI(I) = CINUI*CSCLR 40 CONTINUE RETURN 50 CONTINUE NZ = -1 IF(NW.EQ.(-2)) NZ=-2 RETURN END