Attribute VB_Name = "Module3" DefDbl A-H, O-Z DefInt I-N '************************************************************************** '* Subroutines and Functions used By programs TZBESJ, TZBESK, TZBESY * '* (Evalute Bessel Functions with complex argument, 1st to 3rd kind) * '* ---------------------------------------------------------------------- * '* Reference: AMOS, DONALD E., SANDIA NATIONAL LABORATORIES, 1983. * '* * '* Visual Basic Release By J-P Moreau, Paris (07/13/2005). * '* (www.jpmoreau.fr) * '************************************************************************** 'MODULE CBESS0 ' Subroutines or functions using module COMPLEX 'LIST OF SUBROUTINES OR FUNCTIONS ' Function Sign(a,b) ' Subroutine ZUCHK(YR, YI, NZ, ASCLE, TOL) ' Function DGAMLN(Z, IERR) ' Subroutine ZUNIK(ZRR, ZRI, FNU, IKFLG, IPMTR, TOL, ' INIT, PHIR, PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, ' SUMR, SUMI, TMPR, TMPI) ' Subroutine ZUNHJ(ZR, ZI, FNU, IPMTR, TOL, PHIR, PHII, ARGR, ARGI, ' ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) ' Subroutine ZUOIK(ZR, ZI, FNU, KODE, IKFLG, N, YR, YI, NUF, TOL, ELIM, ALIM) ' Subroutine ZS1S2(ZRR, ZRI, S1R, S1I, S2R, S2I, NZ, ASCLE, ALIM, IUF) ' Subroutine ZSERI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM) ' Subroutine ZASYI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM, ALIM) ' Subroutine ZMLRI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL) ' Subroutine ZSHCH(ZR, ZI, CSHR, CSHI, CCHR, CCHI) ' Subroutine ZRATI(ZR, ZI, FNU, N, CYR, CYI, TOL) Function Sign(a, b) If b < 0# Then Sign = -Abs(a) Else Sign = Abs(a) End If End Function Sub ZUCHK(YR, YI, NZ, ASCLE, TOL) '***BEGIN PROLOGUE ZUCHK '***REFER TO ZSERI,ZUOIK,ZUNK1,ZUNK2,ZUNI1,ZUNI2,ZKSCL ' ' Y ENTERS AS A SCALED QUANTITY WHOSE MAGNITUDE IS GREATER THAN ' EXP(-ALIM)=ASCLE=1000*D1MACH(1)/TOL. THE TEST IS MADE TO SEE ' IF THE MAGNITUDE OF THE REAL OR IMAGINARY PART WOULD UNDERFLOW ' WHEN Y IS SCALED (BY TOL) TO ITS PROPER VALUE. Y IS ACCEPTED ' IF THE UNDERFLOW IS AT LEAST ONE PRECISION BELOW THE MAGNITUDE ' OF THE LARGEST COMPONENT OTHERWISE THE PHASE ANGLE DOES NOT HAVE ' ABSOLUTE ACCURACY AND AN UNDERFLOW IS ASSUMED. ' '***ROUTINES CALLED (NONE) '***END PROLOGUE ZUCHK ' ' COMPLEX Y ' SS, ST, WR, WI: double NZ = 0 WR = Abs(YR) WI = Abs(YI) ST = DMIN(WR, WI) If ST > ASCLE Then GoTo 10 SS = DMAX(WR, WI) ST = ST / TOL If SS < ST Then NZ = 1 10 End Sub 'ZUCHK Function DGAMLN(Z, IERR) '***BEGIN PROLOGUE DGAMLN '***DATE WRITTEN 830501 (YYMMDD) '***REVISION DATE 830501 (YYMMDD) '***CATEGORY NO. B5F '***KEYWORDS GAMMA FUNCTION,LOGARITHM OF GAMMA FUNCTION '***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES '***PURPOSE TO COMPUTE THE LOGARITHM OF THE GAMMA FUNCTION '***DESCRIPTION ' ' **** A DOUBLE PRECISION ROUTINE **** ' DGAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR ' Z.GT.0. THE ASYMPTOTIC EXPANSION IS USED TO GENERATE VALUES ' GREATER THAN ZMIN WHICH ARE ADJUSTED BY THE RECURSION ' G(Z+1)=Z*G(Z) FOR Z <= ZMIN. THE FUNCTION WAS MADE AS ' PORTABLE AS POSSIBLE BY COMPUTIMG ZMIN FROM THE NUMBER OF BASE ' 10 DIGITS IN A WORD, RLN=MAX(-ALOG10(D1MACH(4)),0.5E-18) ' LIMITED TO 18 DIGITS OF (RELATIVE) ACCURACY. ' ' SINCE INTEGER ARGUMENTS ARE COMMON, A TABLE LOOK UP ON 100 ' VALUES IS USED FOR SPEED OF EXECUTION. ' ' DESCRIPTION OF ARGUMENTS ' ' INPUT Z IS D0UBLE PRECISION ' Z - ARGUMENT, Z > 0.0 ' ' OUTPUT DGAMLN IS DOUBLE PRECISION ' DGAMLN - NATURAL LOG OF THE GAMMA FUNCTION AT Z.NE.0.0D0 ' IERR - ERROR FLAG ' IERR=0, NORMAL RETURN, COMPUTATION COMPLETED ' IERR=1, Z <= 0.0, NO COMPUTATION ' ' '***REFERENCES COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT ' BY D. E. AMOS, SAND83-0083, MAY, 1983. '***ROUTINES CALLED I1MACH,D1MACH (of unit COMPLEX). '***END PROLOGUE DGAMLN 'Labels: 1410, 1420, 1440, 1460, 1470 ' CON, FLN, FZ, RLN, S, TLG, TRM, TST, T1, WDTOL, ZDMY, ' ZINC, ZM, ZMIN, ZP, ZSQ: Double ' I, I1M, K, MZ, NZ: Integer Dim CF(22) As Double Dim GLN(100) As Double ' LNGAMMA(N), N=1,100 GLN(1) = 0#: GLN(2) = 0# GLN(3) = 0.693147180559945: GLN(4) = 1.79175946922805 GLN(5) = 3.17805383034795: GLN(6) = 4.78749174278205 GLN(7) = 6.5792512120101: GLN(8) = 8.52516136106541 GLN(9) = 10.6046029027453: GLN(10) = 12.8018274800815 GLN(11) = 15.1044125730755: GLN(12) = 17.5023078458739 GLN(13) = 19.9872144956619: GLN(14) = 22.5521638531234 GLN(15) = 25.1912211827387: GLN(16) = 27.8992713838409 GLN(17) = 30.6718601060807: GLN(18) = 33.5050734501369 GLN(19) = 36.3954452080331: GLN(20) = 39.3398841871995 GLN(21) = 42.3356164607535: GLN(22) = 45.3801388984769 GLN(23) = 48.4711813518352: GLN(24) = 51.6066755677644 GLN(25) = 54.7847293981123: GLN(26) = 58.0036052229805 GLN(27) = 61.261701761002: GLN(28) = 64.5575386270063 GLN(29) = 67.8897431371815: GLN(30) = 71.257038967168 GLN(31) = 74.6582363488302: GLN(32) = 78.0922235533153 GLN(33) = 81.557959456115: GLN(34) = 85.0544670175815 GLN(35) = 88.5808275421977: GLN(36) = 92.1361756036871 GLN(37) = 95.7196945421432: GLN(38) = 99.3306124547874 GLN(39) = 102.968198614514: GLN(40) = 106.631760260644 GLN(41) = 110.320639714757: GLN(42) = 114.034211781462 GLN(43) = 117.771881399745: GLN(44) = 121.533081515439 GLN(45) = 125.317271149357: GLN(46) = 129.123933639127 GLN(47) = 132.952575035616: GLN(48) = 136.802722637326 GLN(49) = 140.673923648234: GLN(50) = 144.565743946345 GLN(51) = 148.477766951773: GLN(52) = 152.409592584497 GLN(53) = 156.360836303079: GLN(54) = 160.331128216631 GLN(55) = 164.320112263195: GLN(56) = 168.327445448428 GLN(57) = 172.352797139163: GLN(58) = 176.395848406997 GLN(59) = 180.456291417544: GLN(60) = 184.53382886145 GLN(61) = 188.628173423672: GLN(62) = 192.739047287845 GLN(63) = 196.86618167289: GLN(64) = 201.009316399281 GLN(65) = 205.168199482641: GLN(66) = 209.342586752537 GLN(67) = 213.532241494563: GLN(68) = 217.736934113954 GLN(69) = 221.95644181913: GLN(70) = 226.190548323728 GLN(71) = 230.439043565777: GLN(72) = 234.701723442818 GLN(73) = 238.978389561834: GLN(74) = 243.268849002983 GLN(75) = 247.572914096187: GLN(76) = 251.890402209723 GLN(77) = 256.221135550009: GLN(78) = 260.564940971863 GLN(79) = 264.921649798553: GLN(80) = 269.29109765102 GLN(81) = 273.673124285694: GLN(82) = 278.067573440366 GLN(83) = 282.47429268763: GLN(84) = 286.893133295427 GLN(85) = 291.32395009427: GLN(86) = 295.766601350761 GLN(87) = 300.220948647014: GLN(88) = 304.686856765669 GLN(89) = 309.164193580147: GLN(90) = 313.652829949879 GLN(91) = 318.152639620209: GLN(92) = 322.663499126726 GLN(93) = 327.185287703775: GLN(94) = 331.717887196928 GLN(95) = 336.261181979199: GLN(96) = 340.815058870799 GLN(97) = 345.379407062267: GLN(98) = 349.95411804077 GLN(99) = 354.539085519441: GLN(100) = 359.134205369575 ' COEFFICIENTS OF ASYMPTOTIC EXPANSION CF(1) = 8.33333333333333E-02: CF(2) = -2.77777777777778E-03 CF(3) = 7.93650793650794E-04: CF(4) = -5.95238095238095E-04 CF(5) = 8.41750841750842E-04: CF(6) = -1.91752691752692E-03 CF(7) = 6.41025641025641E-03: CF(8) = -2.95506535947712E-02 CF(9) = 0.179644372368831: CF(10) = -1.3924322169059 CF(11) = 13.4028640441684: CF(12) = -156.848284626002 CF(13) = 2193.10333333333: CF(14) = -36108.771253725 CF(15) = 691472.268851313: CF(16) = -15238221.5394074 CF(17) = 382900751.391414: CF(18) = -10882266035.7844 CF(19) = 347320283765.002: CF(20) = -12369602142269.3 CF(21) = 488788064793079#: CF(22) = -2.13203339609194E+16 ' LN(2*PI) CON = 1.83787706640935 IERR = 0 If Z <= 0# Then GoTo 1470 If Z > 101# Then GoTo 1410 NZ = Int(Z) FZ = Z - 1# * NZ If FZ > 0# Then GoTo 1410 If NZ > 100 Then GoTo 1410 DGAMLN = GLN(NZ) GoTo 1480 1410 WDTOL = D1MACH(4) WDTOL = DMAX(WDTOL, 5E-19) RLN = D1MACH(5) * XI1MACH(14) FLN = DMIN(RLN, 20#) FLN = DMAX(FLN, 3#) FLN = FLN - 3# ZM = 1.8 + 0.3875 * FLN MZ = Int(ZM) + 1 ZMIN = 1# * MZ ZDMY = Z ZINC = 0# If Z >= ZMIN Then GoTo 1420 ZINC = ZMIN - 1# * NZ ZDMY = Z + ZINC 1420 ZP = 1# / ZDMY T1 = CF(1) * ZP S = T1 If ZP < WDTOL Then GoTo 1440 ZSQ = ZP * ZP TST = T1 * WDTOL For K = 2 To 22 ZP = ZP * ZSQ TRM = CF(K) * ZP If Abs(TRM) < TST Then GoTo 1440 S = S + TRM Next K 1440 If ZINC <> 0# Then GoTo 1460 TLG = Log(Z) DGAMLN = Z * (TLG - 1#) + 0.5 * (CON - TLG) + S GoTo 1480 1460 ZP = 1# NZ = Int(ZINC) For I = 1 To NZ ZP = ZP * (Z + 1# * (I - 1)) Next I TLG = Log(ZDMY) DGAMLN = ZDMY * (TLG - 1#) - Log(ZP) + 0.5 * (CON - TLG) + S GoTo 1480 1470 IERR = 1 1480 End Function 'DGAMLN Sub ZUNIK(ZRR, ZRI, FNU, IKFLG, IPMTR, TOL, INIT, PHIR, PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, TMPR(), TMPI()) '***BEGIN PROLOGUE ZUNIK '***REFER TO ZBESI,ZBESK ' ' ZUNIK COMPUTES PARAMETERS FOR THE UNIFORM ASYMPTOTIC ' EXPANSIONS OF THE I AND K FUNCTIONS ON IKFLG= 1 OR 2 ' RESPECTIVELY BY: ' ' W(FNU,ZR) = PHI*EXP(ZETA)*SUM ' ' WHERE ZETA=-ZETA1 + ZETA2 OR ' ZETA1 - ZETA2 ' ' THE FIRST CALL MUST HAVE INIT=0. SUBSEQUENT CALLS WITH THE ' SAME ZR AND FNU WILL RETURN THE I OR K FUNCTION ON IKFLG= ' 1 OR 2 WITH NO CHANGE IN INIT. TMPR IS A COMPLEX WORK ' ARRAY. IPMTR=0 COMPUTES ALL PARAMETERS. IPMTR=1 COMPUTES PHI, ' ZETA1,ZETA2. ' '***ROUTINES CALLED ZDIV,ZLOG,ZSQRT '***END PROLOGUE ZUNIK ' COMPLEX CFN,CON,CONE,CRFN,TMPR,CZERO,P,HI,S,SR,SUM,T,T2,ZETA1, ' ZETA2,ZN,ZR 'Labels: 1510,1520,1530 ' AC, CONEI, CONER, CRFNI, CRFNR, RFN, SI, SR, SRI, SRR, STI, STR, ' TEST, TI, TR, T2I, T2R, ZEROI, ZEROR, ZNI, ZNR: Double ' I, IDUM, J, K, L: Integer Dim C(120) As Double Dim CON(2) As Double Dim STR As Double ZEROR = 0#: ZEROI = 0#: CONER = 1#: CONEI = 0# CON(1) = 0.398942280401433 CON(2) = 1.2533141373155 'Initialize table C C(1) = 1#: C(2) = -0.208333333333333 C(3) = 0.125: C(4) = 0.334201388888889 C(5) = -0.401041666666667: C(6) = 0.0703125 C(7) = -1.02581259645062: C(8) = 1.84646267361111 C(9) = -0.8912109375: C(10) = 0.0732421875 C(11) = 4.66958442342625: C(12) = -11.207002616223 C(13) = 8.78912353515625: C(14) = -2.3640869140625 C(15) = 0.112152099609375: C(16) = -28.2120725582002 C(17) = 84.6362176746007: C(18) = -91.81824154324 C(19) = 42.5349987453885: C(20) = -7.36879435947963 C(21) = 0.227108001708984: C(22) = 212.570130039217 C(23) = -765.252468141182: C(24) = 1059.990452528 C(25) = -699.579627376133: C(26) = 218.190511744212 C(27) = -26.4914304869516: C(28) = 0.572501420974731 C(29) = -1919.45766231841: C(30) = 8061.72218173731 C(31) = -13586.5500064341: C(32) = 11655.3933368645 C(33) = -5305.6469786134: C(34) = 1200.90291321635 C(35) = -108.090919788395: C(36) = 1.72772750258446 C(37) = 20204.2913309661: C(38) = -96980.5983886375 C(39) = 192547.001232532: C(40) = -203400.177280416 C(41) = 122200.464983017: C(42) = -41192.6549688976 C(43) = 7109.51430248936: C(44) = -493.915304773088 C(45) = 6.07404200127348: C(46) = -242919.187900551 C(47) = 1311763.61466298: C(48) = -2998015.91853811 C(49) = 3763271.2976564: C(50) = -2813563.22658653 C(51) = 1268365.27332162: C(52) = -331645.172484564 C(53) = 45218.7689813627: C(54) = -2499.83048181121 C(55) = 24.3805296995561: C(56) = 3284469.85307204 C(57) = -19706819.1184322: C(58) = 50952602.4926646 C(59) = -74105148.2115327: C(60) = 66344512.274729 C(61) = -37567176.6607634: C(62) = 13288767.1664218 C(63) = -2785618.12808645: C(64) = 308186.404612662 C(65) = -13886.089753717: C(66) = 110.017140269247 C(67) = -49329253.66451: C(68) = 325573074.185766 C(69) = -939462359.681578: C(70) = 1553596899.57058 C(71) = -1621080552.10834: C(72) = 1106842816.82301 C(73) = -495889784.27503: C(74) = 142062907.797533 C(75) = -24474062.7257387: C(76) = 2243768.17792245 C(77) = -84005.4336030241: C(78) = 551.335896122021 C(79) = 814789096.118312: C(80) = -5866481492.05185 C(81) = 18688207509.2958: C(82) = -34632043388.1588 C(83) = 41280185579.754: C(84) = -33026599749.8007 C(85) = 17954213731.1556: C(86) = -6563293792.61928 C(87) = 1559279864.87926: C(88) = -225105661.889415 C(89) = 17395107.5539782: C(90) = -549842.327572289 C(91) = 3038.09051092238: C(92) = -14679261247.6956 C(93) = 114498237732.026: C(94) = -399096175224.466 C(95) = 819218669548.577: C(96) = -1098375156081.22 C(97) = 1008158106865.38: C(98) = -645364869245.376 C(99) = 287900649906.151: C(100) = -87867072178.0233 C(101) = 17634730606.835: C(102) = -2167164983.22379 C(103) = 143157876.718889: C(104) = -3871833.44257261 C(105) = 18257.7554742932: C(106) = 286464035717.679 C(107) = -2406297900028.5: C(108) = 9109341185239.9 C(109) = -20516899410934.4: C(110) = 30565125519935.3 C(111) = -31667088584785.2: C(112) = 23348364044581.8 C(113) = -12320491305598.3: C(114) = 4612725780849.13 C(115) = -1196552880196.18: C(116) = 205914503232.41 C(117) = -21822927757.5292: C(118) = 1247009293.51271 C(119) = -29188388.1222208: C(120) = 118838.426256783 If INIT <> 0 Then GoTo 1520 '----------------------------------------------------------------------- ' INITIALIZE ALL VARIABLES '----------------------------------------------------------------------- RFN = 1# / FNU TR = ZRR * RFN TI = ZRI * RFN SR = CONER + (TR * TR - TI * TI) SI = CONEI + (TR * TI + TI * TR) ZSQRT SR, SI, SRR, SRI STR = CONER + SRR STI = CONEI + SRI ZDIV STR, STI, TR, TI, ZNR, ZNI ZLOG ZNR, ZNI, STR, STI, IDUM ZETA1R = FNU * STR ZETA1I = FNU * STI ZETA2R = FNU * SRR ZETA2I = FNU * SRI ZDIV CONER, CONEI, SRR, SRI, TR, TI SRR = TR * RFN SRI = TI * RFN ZSQRT SRR, SRI, TMPR(16), TMPI(16) PHIR = TMPR(16) * CON(IKFLG) PHII = TMPI(16) * CON(IKFLG) If IPMTR <> 0 Then GoTo 1540 ZDIV CONER, CONEI, SR, SI, T2R, T2I TMPR(1) = CONER TMPI(1) = CONEI CRFNR = CONER CRFNI = CONEI AC = 1# L = 1 For K = 2 To 15 SR = ZEROR SI = ZEROI For J = 1 To K L = L + 1 STR = SR * T2R - SI * T2I + C(L) SI = SR * T2I + SI * T2R SR = STR Next J STR = CRFNR * SRR - CRFNI * SRI CRFNI = CRFNR * SRI + CRFNI * SRR CRFNR = STR TMPR(K) = CRFNR * SR - CRFNI * SI TMPI(K) = CRFNR * SI + CRFNI * SR AC = AC * RFN TEST = Abs(TMPR(K)) + Abs(TMPI(K)) If AC < TOL And TEST < TOL Then GoTo 1510 Next K K = 15 1510 INIT = K 1520 If IKFLG = 2 Then GoTo 1530 '----------------------------------------------------------------------- ' COMPUTE SUM FOR THE I FUNCTION '----------------------------------------------------------------------- SR = ZEROR SI = ZEROI For I = 1 To INIT SR = SR + TMPR(I) SI = SI + TMPI(I) Next I SUMR = SR SUMI = SI PHIR = TMPR(16) * CON(1) PHII = TMPI(16) * CON(1) GoTo 1540 '----------------------------------------------------------------------- ' COMPUTE SUM FOR THE K FUNCTION '----------------------------------------------------------------------- 1530 SR = ZEROR SI = ZEROI TR = CONER For I = 1 To INIT SR = SR + TR * TMPR(I) SI = SI + TR * TMPI(I) TR = -TR Next I SUMR = SR SUMI = SI PHIR = TMPR(16) * CON(2) PHII = TMPI(16) * CON(2) 1540 End Sub 'ZUNIK Sub ZUNHJ(ZR, ZI, FNU, IPMTR, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) '***BEGIN PROLOGUE ZUNHJ '***REFER TO ZBESI,ZBESK ' ' REFERENCES ' HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ AND I.A. ' STEGUN, AMS55, NATIONAL BUREAU OF STANDARDS, 1965, CHAPTER 9. ' ' ASYMPTOTICS AND SPECIAL FUNCTIONS BY F.W.J. OLVER, ACADEMIC ' PRESS, N.Y., 1974, PAGE 420 ' ' ABSTRACT ' ZUNHJ COMPUTES PARAMETERS FOR BESSEL FUNCTIONS C(FNU,Z) = ' J(FNU,Z), Y(FNU,Z) OR H(I,FNU,Z) I=1,2 FOR LARGE ORDERS FNU ' BY MEANS OF THE UNIFORM ASYMPTOTIC EXPANSION ' ' C(FNU,Z)=C1*PHI*( ASUM*AIRY(ARG) + C2*BSUM*DAIRY(ARG) ) ' ' FOR PROPER CHOICES OF C1, C2, AIRY AND DAIRY WHERE AIRY IS ' AN AIRY FUNCTION AND DAIRY IS ITS DERIVATIVE. ' ' (2/3)*FNU*ZETA^1.5 = ZETA1-ZETA2, ' ' ZETA1=0.5*FNU*CLOG((1+W)/(1-W)), ZETA2=FNU*W FOR SCALING ' PURPOSES IN AIRY FUNCTIONS FROM CAIRY OR CBIRY. ' ' MCONJ=SIGN OF AIMAG(Z), BUT IS AMBIGUOUS WHEN Z IS REAL AND ' MUST BE SPECIFIED. IPMTR=0 RETURNS ALL PARAMETERS. IPMTR=1 ' COMPUTES ALL EXCEPT ASUM AND BSUM. ' '***ROUTINES CALLED ZABS,ZDIV,ZLOG,ZSQRT '***END PROLOGUE ZUNHJ ' COMPLEX ARG,ASUM,BSUM,CFNU,CONE,CR,CZERO,DR,P,PHI,PRZTH,PTFN, ' RFN13,RTZTA,RZTH,SUMA,SUMB,TFN,T2,UP,W,W2,Z,ZA,ZB,ZC,ZETA,ZETA1, ' ZETA2,ZTH 'Labels: 1555, 1558,1560,1562,1564,1566,1568,1570,1572,1574,1576,1580 ' ANG, ATOL, AW2, AZTH, BTOL, CONEI, CONER, EX1, EX2, FN13, FN23, GPI, HPI, ' PP, PRZTHI, PRZTHR, PTFNI, PTFNR, RAW, RAW2, RAZTH, RFNU, RFNU2, RFN13, ' RTZTI, RTZTR, RZTHI, RZTHR, STI, STR, SUMAI, SUMAR, SUMBI, SUMBR, TEST, ' TFNI, TFNR, THPI, TZAI, TZAR, T2I, T2R, WI, WR, W2I, W2R, ZAI, ZAR, ZBI, ' ZBR, ZCI, ZCR, ZEROI, ZEROR, ZETAI, ZETAR, ZTHI, ZTHR: Double ' IAS, IBS, IS, J, JR, JU, K, KMAX, KP1, KS, L, LR, LRP1, L1, L2, M, IDUM: Integer Dim AR1(16), BR1(16), UPR(16), UPI(16), CRR(16), CRI(16), DRR(16), DRI(16) Dim C(105) Dim ALFA(180), BETA(210) Dim GAMA(30), AP(30), PR(30), PI(30) Dim STR As Double 'Initialize tables AR1(1) = 1#: AR1(2) = 0.104166666666667 AR1(3) = 8.35503472222222E-02: AR1(4) = 0.128226574556327 AR1(5) = 0.29184902646414: AR1(6) = 0.881627267443758 AR1(7) = 3.32140828186277: AR1(8) = 14.9957629868626 AR1(9) = 78.9230130115865: AR1(10) = 474.451538868264 AR1(11) = 3207.49009089066: AR1(12) = 24086.549640874 AR1(13) = 198923.11916951: AR1(14) = 1791902.00777534 AR1(15) = 0#: AR1(16) = 0# BR1(1) = 1#: BR1(2) = -0.145833333333333 BR1(3) = -9.87413194444444E-02: BR1(4) = -0.143312053915895 BR1(5) = -0.317227202678414: BR1(6) = -0.94242914795712 BR1(7) = -3.51120304082635: BR1(8) = -15.727263620368 BR1(9) = -82.2814390971859: BR1(10) = -492.355370523671 BR1(11) = -3316.21856854797: BR1(12) = -24827.6742452086 BR1(13) = -204526.58731513: BR1(14) = -1838444.91706821 BR1(15) = 0#: BR1(16) = 0# C(1) = 1#: C(2) = -0.208333333333333 C(3) = 0.125: C(4) = 0.334201388888889 C(5) = -0.401041666666667: C(6) = 0.0703125 C(7) = -1.02581259645062: C(8) = 1.84646267361111 C(9) = -0.8912109375: C(10) = 0.0732421875 C(11) = 4.66958442342625: C(12) = -11.207002616223 C(13) = 8.78912353515625: C(14) = -2.3640869140625 C(15) = 0.112152099609375: C(16) = -28.2120725582002 C(17) = 84.6362176746007: C(18) = -91.81824154324 C(19) = 42.5349987453885: C(20) = -7.36879435947963 C(21) = 0.227108001708984: C(22) = 212.570130039217 C(23) = -765.252468141182: C(24) = 1059.990452528 C(25) = -699.579627376133: C(26) = 218.190511744212 C(27) = -26.4914304869516: C(28) = 0.572501420974731 C(29) = -1919.45766231841: C(30) = 8061.72218173731 C(31) = -13586.5500064341: C(32) = 11655.3933368645 C(33) = -5305.6469786134: C(34) = 1200.90291321635 C(35) = -108.090919788395: C(36) = 1.72772750258446 C(37) = 20204.2913309661: C(38) = -96980.5983886375 C(39) = 192547.001232532: C(40) = -203400.177280416 C(41) = 122200.464983017: C(42) = -41192.6549688976 C(43) = 7109.51430248936: C(44) = -493.915304773088 C(45) = 6.07404200127348: C(46) = -242919.187900551 C(47) = 1311763.61466298: C(48) = -2998015.91853811 C(49) = 3763271.2976564: C(50) = -2813563.22658653 C(51) = 1268365.27332162: C(52) = -331645.172484564 C(53) = 45218.7689813627: C(54) = -2499.83048181121 C(55) = 24.3805296995561: C(56) = 3284469.85307204 C(57) = -19706819.1184322: C(58) = 50952602.4926646 C(59) = -74105148.2115327: C(60) = 66344512.274729 C(61) = -37567176.6607634: C(62) = 13288767.1664218 C(63) = -2785618.12808645: C(64) = 308186.404612662 C(65) = -13886.089753717: C(66) = 110.017140269247 C(67) = -49329253.66451: C(68) = 325573074.185766 C(69) = -939462359.681578: C(70) = 1553596899.57058 C(71) = -1621080552.10834: C(72) = 1106842816.82301 C(73) = -495889784.27503: C(74) = 142062907.797533 C(75) = -24474062.7257387: C(76) = 2243768.17792245 C(77) = -84005.4336030241: C(78) = 551.335896122021 C(79) = 814789096.118312: C(80) = -5866481492.05185 C(81) = 18688207509.2958: C(82) = -34632043388.1588 C(83) = 41280185579.754: C(84) = -33026599749.8007 C(85) = 17954213731.1556: C(86) = -6563293792.61928 C(87) = 1559279864.87926: C(88) = -225105661.889415 C(89) = 17395107.5539782: C(90) = -549842.327572289 C(91) = 3038.09051092238: C(92) = -14679261247.6956 C(93) = 114498237732.026: C(94) = -399096175224.466 C(95) = 819218669548.577: C(96) = -1098375156081.22 C(97) = 1008158106865.38: C(98) = -645364869245.376 C(99) = 287900649906.151: C(100) = -87867072178.0233 C(101) = 17634730606.835: C(102) = -2167164983.22379 C(103) = 143157876.718889: C(104) = -3871833.44257261 C(105) = 18257.7554742932: ALFA(1) = -4.44444444444444E-03: ALFA(2) = -9.22077922077922E-04 ALFA(3) = -8.84892884892885E-05: ALFA(4) = 1.6592768783245E-04 ALFA(5) = 2.46691372741793E-04: ALFA(6) = 2.65995589346255E-04 ALFA(7) = 2.61824297061501E-04: ALFA(8) = 2.48730437344656E-04 ALFA(9) = 2.32721040083232E-04: ALFA(10) = 2.16362485712365E-04 ALFA(11) = 2.00738858762752E-04: ALFA(12) = 1.86267636637545E-04 ALFA(13) = 1.73060775917876E-04: ALFA(14) = 1.61091705929016E-04 ALFA(15) = 1.50274774160908E-04: ALFA(16) = 1.4050349739127E-04 ALFA(17) = 1.31668816545923E-04: ALFA(18) = 1.23667445598253E-04 ALFA(19) = 1.16405271474738E-04: ALFA(20) = 1.09798298372713E-04 ALFA(21) = 1.03772410422993E-04: ALFA(22) = 9.82626078369363E-05 ALFA(23) = 9.32120517249503E-05: ALFA(24) = 8.85710852478712E-05 ALFA(25) = 8.429631057157E-05: ALFA(26) = 8.03497548407791E-05 ALFA(27) = 7.66981345359207E-05: ALFA(28) = 7.33122157481778E-05 ALFA(29) = 7.01662625163141E-05: ALFA(30) = 6.7237563379016E-05 ALFA(31) = 6.93735541354589E-04: ALFA(32) = 2.32241745182922E-04 ALFA(33) = -1.41986273556691E-05: ALFA(34) = -1.16444931672049E-04 ALFA(35) = -1.50803558053049E-04: ALFA(36) = -1.55121924918096E-04 ALFA(37) = -1.46809756646466E-04: ALFA(38) = -1.33815503867491E-04 ALFA(39) = -1.19744975684254E-04: ALFA(40) = -1.06184319207974E-04 ALFA(41) = -9.37699549891194E-05: ALFA(42) = -8.26923045588193E-05 ALFA(43) = -7.29374348155221E-05: ALFA(44) = -6.44042357721016E-05 ALFA(45) = -5.69611566009369E-05: ALFA(46) = -5.04731044303562E-05 ALFA(47) = -4.48134868008883E-05: ALFA(48) = -3.98688727717599E-05 ALFA(49) = -3.55400532972043E-05: ALFA(50) = -3.17414256609023E-05 ALFA(51) = -2.83996793904175E-05: ALFA(52) = -2.54522720634871E-05 ALFA(53) = -2.28459297164725E-05: ALFA(54) = -2.05352753106481E-05 ALFA(55) = -1.84816217627666E-05: ALFA(56) = -1.66519330021394E-05 ALFA(57) = -1.50179412980119E-05: ALFA(58) = -1.35554031379041E-05 ALFA(59) = -1.22434746473858E-05: ALFA(60) = -1.10641884811308E-05 ALFA(61) = -3.54211971457744E-04: ALFA(62) = -1.56161263945159E-04 ALFA(63) = 3.04465503594936E-05: ALFA(64) = 1.30198655773243E-04 ALFA(65) = 1.67471106699712E-04: ALFA(66) = 1.70222587683593E-04 ALFA(67) = 1.56501427608595E-04: ALFA(68) = 1.36339170977445E-04 ALFA(69) = 1.14886692029825E-04: ALFA(70) = 9.45869093034688E-05 ALFA(71) = 7.64498419250898E-05: ALFA(72) = 6.07570334965197E-05 ALFA(73) = 4.74394299290509E-05: ALFA(74) = 3.62757512005344E-05 ALFA(75) = 2.69939714979225E-05: ALFA(76) = 1.93210938247939E-05 ALFA(77) = 1.30056674793963E-05: ALFA(78) = 7.82620866744497E-06 ALFA(79) = 3.59257485819352E-06: ALFA(80) = 1.44040049814252E-07 ALFA(81) = -2.65396769697939E-06: ALFA(82) = -4.91346867098486E-06 ALFA(83) = -6.72739296091248E-06: ALFA(84) = -8.17269379678658E-06 ALFA(85) = -9.31304715093561E-06: ALFA(86) = -1.02011418798016E-05 ALFA(87) = -1.08805962510593E-05: ALFA(88) = -1.13875481509604E-05 ALFA(89) = -1.17519675674556E-05: ALFA(90) = -1.19987364870944E-05 ALFA(91) = 3.78194199201773E-04: ALFA(92) = 2.02471952761816E-04 ALFA(93) = -6.37938506318862E-05: ALFA(94) = -2.38598230603006E-04 ALFA(95) = -3.10916256027362E-04: ALFA(96) = -3.13680115247576E-04 ALFA(97) = -2.78950273791323E-04: ALFA(98) = -2.28564082619141E-04 ALFA(99) = -1.75245280340847E-04: ALFA(100) = -1.2554406306069E-04 ALFA(101) = -8.22982872820208E-05: ALFA(102) = -4.62860730588116E-05 ALFA(103) = -1.72334302366962E-05: ALFA(104) = 5.60690482304602E-06 ALFA(105) = 2.31395443148287E-05: ALFA(106) = 3.62642745856794E-05 ALFA(107) = 4.58006124490189E-05: ALFA(108) = 5.24595294959114E-05 ALFA(109) = 5.68396208545815E-05: ALFA(110) = 5.94349820393104E-05 ALFA(111) = 6.06478527578422E-05: ALFA(112) = 6.08023907788436E-05 ALFA(113) = 6.0157789453946E-05: ALFA(114) = 5.89199657344698E-05 ALFA(115) = 5.72515823777593E-05: ALFA(116) = 5.52804375585853E-05 ALFA(117) = 5.3106377380288E-05: ALFA(118) = 5.08069302012326E-05 ALFA(119) = 4.84418647620095E-05: ALFA(120) = 4.60568581607475E-05 ALFA(121) = -6.91141397288294E-04: ALFA(122) = -4.29976633058872E-04 ALFA(123) = 1.83067735980039E-04: ALFA(124) = 6.60088147542014E-04 ALFA(125) = 8.75964969951186E-04: ALFA(126) = 8.77335235958236E-04 ALFA(127) = 7.49369585378991E-04: ALFA(128) = 5.63832329756981E-04 ALFA(129) = 3.68059319971443E-04: ALFA(130) = 1.88464535514456E-04 ALFA(131) = 3.70663057664904E-05: ALFA(132) = -8.28520220232137E-05 ALFA(133) = -1.72751952869173E-04: ALFA(134) = -2.36314873605873E-04 ALFA(135) = -2.77966150694907E-04: ALFA(136) = -3.02079514155457E-04 ALFA(137) = -3.1259471264382E-04: ALFA(138) = -3.12872558758067E-04 ALFA(139) = -3.05678038466324E-04: ALFA(140) = -2.93226470614557E-04 ALFA(141) = -2.77255655582935E-04: ALFA(142) = -2.59103928467032E-04 ALFA(143) = -2.3978401439648E-04: ALFA(144) = -2.20048260045423E-04 ALFA(145) = -2.00443911094971E-04: ALFA(146) = -1.81358692210971E-04 ALFA(147) = -1.63057674478657E-04: ALFA(148) = -1.45712672175206E-04 ALFA(149) = -1.29425421983925E-04: ALFA(150) = -1.14245691942446E-04 ALFA(151) = 1.92821964248776E-03: ALFA(152) = 1.35592576302022E-03 ALFA(153) = -7.17858090421303E-04: ALFA(154) = -2.5808480257527E-03 ALFA(155) = -3.49271130826168E-03: ALFA(156) = -3.46986299340961E-03 ALFA(157) = -2.8228523335131E-03: ALFA(158) = -1.88103076404891E-03 ALFA(159) = -8.89531718383948E-04: ALFA(160) = 3.87912102631035E-06 ALFA(161) = 7.28688540119691E-04: ALFA(162) = 1.26566373053458E-03 ALFA(163) = 1.62518158372674E-03: ALFA(164) = 1.83203153216373E-03 ALFA(165) = 1.91588388990528E-03: ALFA(166) = 1.90588846755546E-03 ALFA(167) = 1.82798982421826E-03: ALFA(168) = 1.70389506421122E-03 ALFA(169) = 1.55097127171098E-03: ALFA(170) = 1.38261421852276E-03 ALFA(171) = 1.20881424230065E-03: ALFA(172) = 1.03676532638345E-03 ALFA(173) = 8.71437918068619E-04: ALFA(174) = 7.16080155297701E-04 ALFA(175) = 5.72637002558129E-04: ALFA(176) = 4.42089819465802E-04 ALFA(177) = 3.24724948503091E-04: ALFA(178) = 2.20342042730247E-04 ALFA(179) = 1.28412898401354E-04: ALFA(180) = 4.82005924552095E-05 BETA(1) = 1.79988721413553E-02: BETA(2) = 5.59964911064388E-03 BETA(3) = 2.88501402231133E-03: BETA(4) = 1.80096606761054E-03 BETA(5) = 1.24753110589199E-03: BETA(6) = 9.22878876572938E-04 BETA(7) = 7.14430421727287E-04: BETA(8) = 5.71787281789705E-04 BETA(9) = 4.69431007606482E-04: BETA(10) = 3.93232835462917E-04 BETA(11) = 3.34818889318298E-04: BETA(12) = 2.88952148495752E-04 BETA(13) = 2.52211615549573E-04: BETA(14) = 2.22280580798883E-04 BETA(15) = 1.97541838033063E-04: BETA(16) = 1.76836855019718E-04 BETA(17) = 1.59316899661821E-04: BETA(18) = 1.44347930197334E-04 BETA(19) = 1.31448068119965E-04: BETA(20) = 1.20245444949303E-04 BETA(21) = 1.10449144504599E-04: BETA(22) = 1.01828770740567E-04 BETA(23) = 9.41998224204238E-05: BETA(24) = 8.74130545753834E-05 BETA(25) = 8.13466262162801E-05: BETA(26) = 7.59002269646219E-05 BETA(27) = 7.09906300634154E-05: BETA(28) = 6.65482874842468E-05 BETA(29) = 6.25146958969275E-05: BETA(30) = 5.88403394426252E-05 BETA(31) = -1.49282953213429E-03: BETA(32) = -8.78204709546389E-04 BETA(33) = -5.02916549572035E-04: BETA(34) = -2.94822138512746E-04 BETA(35) = -1.75463996970783E-04: BETA(36) = -1.04008550460816E-04 BETA(37) = -5.96141953046458E-05: BETA(38) = -3.12038929076098E-05 BETA(39) = -1.2608973598023E-05: BETA(40) = -2.4289260857573E-07 BETA(41) = 8.05996165414274E-06: BETA(42) = 1.36507009262147E-05 BETA(43) = 1.73964125472926E-05: BETA(44) = 1.98672978842134E-05 BETA(45) = 2.14463263790823E-05: BETA(46) = 2.23954659232457E-05 BETA(47) = 2.28967783814713E-05: BETA(48) = 2.30785389811178E-05 BETA(49) = 2.30321976080909E-05: BETA(50) = 2.28236073720349E-05 BETA(51) = 2.25005881105292E-05: BETA(52) = 2.20981015361991E-05 BETA(53) = 2.16418427448104E-05: BETA(54) = 2.11507649256221E-05 BETA(55) = 2.06388749782171E-05: BETA(56) = 2.01165241997082E-05 BETA(57) = 1.95913450141179E-05: BETA(58) = 1.90689367910437E-05 BETA(59) = 1.85533719641637E-05: BETA(60) = 1.80475722259674E-05 BETA(61) = 5.52213076721293E-04: BETA(62) = 4.47932581552385E-04 BETA(63) = 2.79520653992021E-04: BETA(64) = 1.52468156198447E-04 BETA(65) = 6.93271105657044E-05: BETA(66) = 1.76258683069991E-05 BETA(67) = -1.35744996343269E-05: BETA(68) = -3.17972413350427E-05 BETA(69) = -4.18861861696693E-05: BETA(70) = -4.69004889379141E-05 BETA(71) = -4.87665447413787E-05: BETA(72) = -4.87010031186735E-05 BETA(73) = -4.74755620890087E-05: BETA(74) = -4.55813058138628E-05 BETA(75) = -4.33309644511266E-05: BETA(76) = -4.0923019315775E-05 BETA(77) = -3.84822638603221E-05: BETA(78) = -3.60857167535411E-05 BETA(79) = -3.37793306123367E-05: BETA(80) = -3.1588856077211E-05 BETA(81) = -2.95269561750807E-05: BETA(82) = -2.75978914828336E-05 BETA(83) = -2.58006174666884E-05: BETA(84) = -2.4130835676128E-05 BETA(85) = -2.25823509518346E-05: BETA(86) = -2.11479656768913E-05 BETA(87) = -1.98200638885295E-05: BETA(88) = -1.85909870801065E-05 BETA(89) = -1.7453269984421E-05: BETA(90) = -1.63997823854498E-05 BETA(91) = -4.7461779655996E-04: BETA(92) = -4.77864567147322E-04 BETA(93) = -3.20390228067038E-04: BETA(94) = -1.61105016119962E-04 BETA(95) = -4.25778101285435E-05: BETA(96) = 3.44571294294967E-05 BETA(97) = 7.97092684075675E-05: BETA(98) = 1.03138236708272E-04 BETA(99) = 1.12466775262204E-04: BETA(100) = 1.13103642108481E-04 BETA(101) = 1.08651634848774E-04: BETA(102) = 1.01437951597662E-04 BETA(103) = 9.29298396593364E-05: BETA(104) = 8.4029313301609E-05 BETA(105) = 7.52727991349134E-05: BETA(106) = 6.69632521975731E-05 BETA(107) = 5.92564547323195E-05: BETA(108) = 5.22169308826976E-05 BETA(109) = 4.58539485165361E-05: BETA(110) = 4.01445513891487E-05 BETA(111) = 3.50481730031328E-05: BETA(112) = 3.05157995034347E-05 BETA(113) = 2.64956119950516E-05: BETA(114) = 2.29363633690998E-05 BETA(115) = 1.97893056664022E-05: BETA(116) = 1.70091984636413E-05 BETA(117) = 1.45547428261524E-05: BETA(118) = 1.23886640995878E-05 BETA(119) = 1.04775876076583E-05: BETA(120) = 8.79179954978479E-06 BETA(121) = 7.36465810572578E-04: BETA(122) = 8.72790805146194E-04 BETA(123) = 6.22614862573135E-04: BETA(124) = 2.85998154194304E-04 BETA(125) = 3.84737672879366E-06: BETA(126) = -1.87906003636972E-04 BETA(127) = -2.97603646594555E-04: BETA(128) = -3.45998126832656E-04 BETA(129) = -3.53382470916038E-04: BETA(130) = -3.35715635775049E-04 BETA(131) = -3.0432112478904E-04: BETA(132) = -2.66722723047613E-04 BETA(133) = -2.2765421412282E-04: BETA(134) = -1.89922611854562E-04 BETA(135) = -1.55058918599094E-04: BETA(136) = -1.23778240761874E-04 BETA(137) = -9.62926147717644E-05: BETA(138) = -7.25178327714425E-05 BETA(139) = -5.22070028895634E-05: BETA(140) = -3.50347750511901E-05 BETA(141) = -2.06489761035552E-05: BETA(142) = -8.70106096849767E-06 BETA(143) = 1.136986866751E-06: BETA(144) = 9.16426474122779E-06 BETA(145) = 1.56477785428873E-05: BETA(146) = 2.08223629482467E-05 BETA(147) = 2.48923381004595E-05: BETA(148) = 2.80340509574146E-05 BETA(149) = 3.03987774629862E-05: BETA(150) = 3.21156731406701E-05 BETA(151) = -1.80182191963886E-03: BETA(152) = -2.43402962938043E-03 BETA(153) = -1.83422663549857E-03: BETA(154) = -7.6220459635401E-04 BETA(155) = 2.39079475256927E-04: BETA(156) = 9.49266117176881E-04 BETA(157) = 1.3446744970154E-03: BETA(158) = 1.48457495259449E-03 BETA(159) = 1.44732339830618E-03: BETA(160) = 1.30268261285657E-03 BETA(161) = 1.10351597375643E-03: BETA(162) = 8.86047440419792E-04 BETA(163) = 6.73073208165665E-04: BETA(164) = 4.77603872856582E-04 BETA(165) = 3.05991926358789E-04: BETA(166) = 1.60315694594722E-04 BETA(167) = 4.00749555270613E-05: BETA(168) = -5.66607461635252E-05 BETA(169) = -1.32506186772983E-04: BETA(170) = -1.90296187989614E-04 BETA(171) = -2.32811450376937E-04: BETA(172) = -2.62628811464669E-04 BETA(173) = -2.82050469867599E-04: BETA(174) = -2.93081563192861E-04 BETA(175) = -2.97435962176317E-04: BETA(176) = -2.96557334239348E-04 BETA(177) = -2.91647363312091E-04: BETA(178) = -2.83696203837734E-04 BETA(179) = -2.73512317095673E-04: BETA(180) = -2.61750155806769E-04 BETA(181) = 6.38585891212051E-03: BETA(182) = 9.62374215806378E-03 BETA(183) = 7.61878061207001E-03: BETA(184) = 2.83219055545628E-03 BETA(185) = -2.0984135201272E-03: BETA(186) = -5.73826764216626E-03 BETA(187) = -7.70804244495415E-03: BETA(188) = -8.21011692264844E-03 BETA(189) = -7.65824520346905E-03: BETA(190) = -6.47209729391045E-03 BETA(191) = -4.99132412004966E-03: BETA(192) = -3.45612289713133E-03 BETA(193) = -2.01785580014171E-03: BETA(194) = -7.59430686781961E-04 BETA(195) = 2.84173631523859E-04: BETA(196) = 1.10891667586337E-03 BETA(197) = 1.72901493872729E-03: BETA(198) = 2.16812590802685E-03 BETA(199) = 2.4535771049454E-03: BETA(200) = 2.61281821058335E-03 BETA(201) = 2.67141039656277E-03: BETA(202) = 2.6520307339598E-03 BETA(203) = 2.57411652877287E-03: BETA(204) = 2.45389126236094E-03 BETA(205) = 2.30460058071796E-03: BETA(206) = 2.13684837686713E-03 BETA(207) = 1.95896528478871E-03: BETA(208) = 1.77737008679454E-03 BETA(209) = 1.59690280765839E-03: BETA(210) = 1.42111975664439E-03 GAMA(1) = 0.629960524947437: GAMA(2) = 0.251984209978975 GAMA(3) = 0.154790300415656: GAMA(4) = 0.110713062416159 GAMA(5) = 8.57309395527395E-02: GAMA(6) = 6.97161316958684E-02 GAMA(7) = 5.86085671893714E-02: GAMA(8) = 5.04698873536311E-02 GAMA(9) = 4.42600580689155E-02: GAMA(10) = 0.039372066154351 GAMA(11) = 3.54283195924455E-02: GAMA(12) = 3.21818857502098E-02 GAMA(13) = 2.94646240791158E-02: GAMA(14) = 2.71581677112934E-02 GAMA(15) = 2.51768272973862E-02: GAMA(16) = 2.34570755306079E-02 GAMA(17) = 2.19508390134907E-02: GAMA(18) = 2.06210828235646E-02 GAMA(19) = 1.94388240897881E-02: GAMA(20) = 1.83810633800683E-02 GAMA(21) = 1.74293213231963E-02: GAMA(22) = 1.65685837786612E-02 GAMA(23) = 1.57865285987918E-02: GAMA(24) = 1.50729501494096E-02 GAMA(25) = 1.44193250839955E-02: GAMA(26) = 1.38184805735342E-02 GAMA(27) = 1.32643378994277E-02: GAMA(28) = 1.27517121970499E-02 GAMA(29) = 1.22761545318763E-02: GAMA(30) = 1.18338262398482E-02 EX1 = 0.333333333333333: EX2 = 0.666666666666667 HPI = 1.5707963267949: GPI = 3.14159265358979 THPI = 4.71238898038469 ZEROR = 0#: ZEROI = 0#: CONER = 1#: CONEI = 0# RFNU = 1# / FNU ZBR = ZR * RFNU ZBI = ZI * RFNU RFNU2 = RFNU * RFNU '----------------------------------------------------------------------- ' COMPUTE IN THE FOURTH QUADRANT '----------------------------------------------------------------------- FN13 = FNU ^ EX1 FN23 = FN13 * FN13 RFN13 = 1# / FN13 W2R = CONER - ZBR * ZBR + ZBI * ZBI W2I = CONEI - ZBR * ZBI - ZBR * ZBI AW2 = ZABS(W2R, W2I) If AW2 > 0.25 Then GoTo 1570 '----------------------------------------------------------------------- ' POWER SERIES FOR CABS(W2) <= 0.25 '----------------------------------------------------------------------- K = 1 PR(1) = CONER PI(1) = CONEI SUMAR = GAMA(1) SUMAI = ZEROI AP(1) = 1# If AW2 < TOL Then GoTo 1555 For K = 2 To 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) < TOL Then GoTo 1555 Next K K = 30 1555 KMAX = K ZETAR = W2R * SUMAR - W2I * SUMAI ZETAI = W2R * SUMAI + W2I * SUMAR ARGR = ZETAR * FN23 ARGI = ZETAI * FN23 ZSQRT SUMAR, SUMAI, ZAR, ZAI 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 ZSQRT ZAR, ZAI, STR, STI PHIR = STR * RFN13 PHII = STI * RFN13 If IPMTR = 1 Then GoTo 1568 '----------------------------------------------------------------------- ' SUM SERIES FOR ASUM AND BSUM '----------------------------------------------------------------------- SUMBR = ZEROR SUMBI = ZEROI For K = 1 To KMAX SUMBR = SUMBR + PR(K) * BETA(K) SUMBI = SUMBI + PI(K) * BETA(K) Next K ASUMR = ZEROR ASUMI = ZEROI BSUMR = SUMBR BSUMI = SUMBI L1 = 0 L2 = 30 BTOL = TOL * (Abs(BSUMR) + Abs(BSUMI)) ATOL = TOL PP = 1# IAS = 0 IBS = 0 If RFNU2 < TOL Then GoTo 1566 For IS0 = 2 To 7 ATOL = ATOL / RFNU2 PP = PP * RFNU2 If IAS = 1 Then GoTo 1560 SUMAR = ZEROR SUMAI = ZEROI For K = 1 To KMAX M = L1 + K SUMAR = SUMAR + PR(K) * ALFA(M) SUMAI = SUMAI + PI(K) * ALFA(M) If AP(K) < ATOL Then GoTo 1558 Next K 1558 ASUMR = ASUMR + SUMAR * PP ASUMI = ASUMI + SUMAI * PP If PP < TOL Then IAS = 1 1560 If IBS = 1 Then GoTo 1564 SUMBR = ZEROR SUMBI = ZEROI For K = 1 To KMAX M = L2 + K SUMBR = SUMBR + PR(K) * BETA(M) SUMBI = SUMBI + PI(K) * BETA(M) If AP(K) < ATOL Then GoTo 1562 Next K 1562 BSUMR = BSUMR + SUMBR * PP BSUMI = BSUMI + SUMBI * PP If PP < BTOL Then IBS = 1 1564 If (IAS = 1) And (IBS = 1) Then GoTo 1566 L1 = L1 + 30 L2 = L2 + 30 Next IS0 1566 ASUMR = ASUMR + CONER PP = RFNU * RFN13 BSUMR = BSUMR * PP BSUMI = BSUMI * PP 1568 GoTo 1590 '----------------------------------------------------------------------- ' CABS(W2) > 0.25 '----------------------------------------------------------------------- 1570 ZSQRT W2R, W2I, WR, WI If WR < 0# Then WR = 0# If WI < 0# Then WI = 0# STR = CONER + WR STI = WI ZDIV STR, STI, ZBR, ZBI, ZAR, ZAI ZLOG ZAR, ZAI, ZCR, ZCI, IDUM If ZCI < 0# Then ZCI = 0# If ZCI > HPI Then ZCI = HPI If ZCR < 0# Then ZCR = 0# ZTHR = (ZCR - WR) * 1.5 ZTHI = (ZCI - WI) * 1.5 ZETA1R = ZCR * FNU ZETA1I = ZCI * FNU ZETA2R = WR * FNU ZETA2I = WI * FNU AZTH = ZABS(ZTHR, ZTHI) ANG = THPI If (ZTHR >= 0#) And (ZTHI < 0#) Then GoTo 1572 ANG = HPI If ZTHR = 0# Then GoTo 1572 ANG = Atn(ZTHI / ZTHR) If ZTHR < 0# Then ANG = ANG + GPI 1572 PP = AZTH ^ EX2 ANG = ANG * EX2 ZETAR = PP * Cos(ANG) ZETAI = PP * Sin(ANG) If ZETAI < 0# Then ZETAI = 0# ARGR = ZETAR * FN23 ARGI = ZETAI * FN23 ZDIV ZTHR, ZTHI, ZETAR, ZETAI, RTZTR, RTZTI ZDIV RTZTR, RTZTI, WR, WI, ZAR, ZAI TZAR = ZAR + ZAR TZAI = ZAI + ZAI ZSQRT TZAR, TZAI, STR, STI PHIR = STR * RFN13 PHII = STI * RFN13 If IPMTR = 1 Then GoTo 1568 RAW = 1# / Sqr(AW2) STR = WR * RAW STI = -WI * RAW TFNR = STR * RFNU * RAW TFNI = STI * RFNU * RAW RAZTH = 1# / AZTH STR = ZTHR * RAZTH STI = -ZTHI * RAZTH RZTHR = STR * RAZTH * RFNU RZTHI = STI * RAZTH * RFNU ZCR = RZTHR * AR1(2) ZCI = RZTHI * AR1(2) RAW2 = 1# / 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 < TOL Then GoTo 1580 PRZTHR = RZTHR PRZTHI = RZTHI PTFNR = TFNR PTFNI = TFNI UPR(1) = CONER UPI(1) = CONEI PP = 1# BTOL = TOL * (Abs(BSUMR) + Abs(BSUMI)) KS = 0 KP1 = 2 L = 3 IAS = 0 IBS = 0 LR = 2 While LR <= 12 LRP1 = LR + 1 '----------------------------------------------------------------------- ' COMPUTE TWO ADDITIONAL CR, DR, AND UP FOR TWO MORE TERMS IN ' NEXT SUMA AND SUMB '----------------------------------------------------------------------- For K = LR To LRP1 KS = KS + 1 KP1 = KP1 + 1 L = L + 1 ZAR = C(L) ZAI = ZEROI For J = 2 To KP1 L = L + 1 STR = ZAR * T2R - T2I * ZAI + C(L) ZAI = ZAR * T2I + ZAI * T2R ZAR = STR Next J 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 * BR1(KS + 1) CRI(KS) = PRZTHI * BR1(KS + 1) STR = PRZTHR * RZTHR - PRZTHI * RZTHI PRZTHI = PRZTHR * RZTHI + PRZTHI * RZTHR PRZTHR = STR DRR(KS) = PRZTHR * AR1(KS + 2) DRI(KS) = PRZTHI * AR1(KS + 2) Next K PP = PP * RFNU2 If IAS = 1 Then GoTo 1574 SUMAR = UPR(LRP1) SUMAI = UPI(LRP1) JU = LRP1 For JR = 1 To LR JU = JU - 1 SUMAR = SUMAR + CRR(JR) * UPR(JU) - CRI(JR) * UPI(JU) SUMAI = SUMAI + CRR(JR) * UPI(JU) + CRI(JR) * UPR(JU) Next JR ASUMR = ASUMR + SUMAR ASUMI = ASUMI + SUMAI TEST = Abs(SUMAR) + Abs(SUMAI) If (PP < TOL) And (TEST < TOL) Then IAS = 1 1574 If IBS = 1 Then GoTo 1576 SUMBR = UPR(LR + 2) + UPR(LRP1) * ZCR - UPI(LRP1) * ZCI SUMBI = UPI(LR + 2) + UPR(LRP1) * ZCI + UPI(LRP1) * ZCR JU = LRP1 For JR = 1 To LR JU = JU - 1 SUMBR = SUMBR + DRR(JR) * UPR(JU) - DRI(JR) * UPI(JU) SUMBI = SUMBI + DRR(JR) * UPI(JU) + DRI(JR) * UPR(JU) Next JR BSUMR = BSUMR + SUMBR BSUMI = BSUMI + SUMBI TEST = Abs(SUMBR) + Abs(SUMBI) If (PP < BTOL) And (TEST < BTOL) Then IBS = 1 1576 If (IAS = 1) And (IBS = 1) Then GoTo 1580 LR = LR + 2 Wend 1580 ASUMR = ASUMR + CONER STR = -BSUMR * RFN13 STI = -BSUMI * RFN13 ZDIV STR, STI, RTZTR, RTZTI, BSUMR, BSUMI GoTo 1568 1590 End Sub 'ZUNHJ Sub ZUOIK(ZR, ZI, FNU, KODE, IKFLG, N, YR(), YI(), NUF, TOL, ELIM, ALIM) '***BEGIN PROLOGUE ZUOIK '***REFER TO ZBESI,ZBESK,ZBESH ' ' ZUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC ' EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM ' (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW ' WHERE ALIM < ELIM. IF THE MAGNITUDE, BASED ON THE LEADING ' EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN ' THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER ' MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE ' EXP(-ELIM)=SMALLEST MACHINE NUMBER*1E+03 AND EXP(-ALIM)= ' EXP(-ELIM)/TOL. ' ' IKFLG=1 MEANS THE I SEQUENCE IS TESTED ' =2 MEANS THE K SEQUENCE IS TESTED ' NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE ' =-1 MEANS AN OVERFLOW WOULD OCCUR ' IKFLG=1 AND NUF > 0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO ' THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE ' IKFLG=2 AND NUF = N MEANS ALL Y VALUES WERE SET TO ZERO ' IKFLG=2 AND 0 < NUF < N NOT CONSIDERED. Y MUST BE SET BY ' ANOTHER ROUTINE. ' '***ROUTINES CALLED ZUCHK,ZUNHJ,ZUNIK,D1MACH,ZABS,ZLOG '***END PROLOGUE ZUOIK ' COMPLEX ARG,ASUM,BSUM,TMP,CZ,CZERO,PHI,SUM,Y,Z,ZB,ZETA1,ZETA2,ZN,ZR 'Labels 1602,1604,1606,1608,1610,1612,1614,1616,1618,1620,1622, ' 1624,1626,1628,1630,1632,1634,1638,1640,1645 ' AARG, AIC, APHI, ARGI, ARGR, ASUMI, ASUMR, ASCLE, AX, AY, ' BSUMI, BSUMR, CZI, CZR, FNN, GNN, GNU, PHII, PHIR, RCZ, ' STR, STI, SUMI, SUMR, ZBI, ZBR, ZEROI, ZEROR, ZETA1I, ' ZETA1R, ZETA2I, ZETA2R, ZNI, ZNR, ZRI, ZRR: Double ' I, IDUM, IFORM, INIT, NN, NW: Integer Dim TMPR(16), TMPI(16) Dim STR As Double ZEROR = 0#: ZEROI = 0# AIC = 1.26551212348465 NUF = 0 NN = N ZRR = ZR ZRI = ZI If ZR >= 0# Then GoTo 1602 ZRR = -ZR ZRI = -ZI 1602 ZBR = ZRR ZBI = ZRI AX = Abs(ZR) * 1.7321 AY = Abs(ZI) IFORM = 1 If AY > AX Then IFORM = 2 GNU = DMAX(FNU, 1#) If IKFLG = 1 Then GoTo 1604 FNN = 1# * NN GNN = FNU + FNN - 1# GNU = DMAX(GNN, FNN) '----------------------------------------------------------------------- ' ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE ' REAL PARTS OF ZETA1, ZETA2 AND ZB. NO ATTEMPT IS MADE TO GET ' THE SIGN OF THE IMAGINARY PART CORRECT. '----------------------------------------------------------------------- 1604 If IFORM = 2 Then GoTo 1606 INIT = 0 ZUNIK ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, PHIR, PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, TMPR(), TMPI() CZR = -ZETA1R + ZETA2R CZI = -ZETA1I + ZETA2I GoTo 1610 1606 ZNR = ZRI ZNI = -ZRR If ZI > 0# Then GoTo 1608 ZNR = -ZNR 1608 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) 1610 If KODE = 1 Then GoTo 1612 CZR = CZR - ZBR CZI = CZI - ZBI 1612 If IKFLG = 1 Then GoTo 1614 CZR = -CZR CZI = -CZI 1614 APHI = ZABS(PHIR, PHII) RCZ = CZR '----------------------------------------------------------------------- ' OVERFLOW TEST '----------------------------------------------------------------------- If RCZ > ELIM Then GoTo 1645 If RCZ < ALIM Then GoTo 1616 RCZ = RCZ + Log(APHI) If IFORM = 2 Then RCZ = RCZ - 0.25 * Log(AARG) - AIC If RCZ > ELIM Then GoTo 1645 GoTo 1624 '----------------------------------------------------------------------- ' UNDERFLOW TEST '----------------------------------------------------------------------- 1616 If RCZ < -ELIM Then GoTo 1618 If RCZ > -ALIM Then GoTo 1624 RCZ = RCZ + Log(APHI) If IFORM = 2 Then RCZ = RCZ - 0.25 * Log(AARG) - AIC If RCZ > -ELIM Then GoTo 1620 1618 For I = 1 To NN YR(I) = ZEROR YI(I) = ZEROI Next I NUF = NN GoTo 1647 1620 ASCLE = 1000 * D1MACH(1) / TOL ZLOG PHIR, PHII, STR, STI, IDUM CZR = CZR + STR CZI = CZI + STI If IFORM = 1 Then GoTo 1622 ZLOG ARGR, ARGI, STR, STI, IDUM CZR = CZR - 0.25 * STR - AIC CZI = CZI - 0.25 * STI 1622 AX = Exp(RCZ) / TOL AY = CZI CZR = AX * Cos(AY) CZI = AX * Sin(AY) ZUCHK CZR, CZI, NW, ASCLE, TOL If NW <> 0 Then GoTo 1618 1624 If IKFLG = 2 Then GoTo 1647 If N = 1 Then GoTo 1647 '----------------------------------------------------------------------- ' SET UNDERFLOWS ON I SEQUENCE '----------------------------------------------------------------------- 1626 GNU = FNU + 1# * (NN - 1) If IFORM = 2 Then GoTo 1628 INIT = 0 ZUNIK ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, PHIR, PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, TMPR(), TMPI() CZR = -ZETA1R + ZETA2R CZI = -ZETA1I + ZETA2I GoTo 1630 1628 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) 1630 If KODE = 1 Then GoTo 1632 CZR = CZR - ZBR CZI = CZI - ZBI 1632 APHI = ZABS(PHIR, PHII) RCZ = CZR If RCZ < -ELIM Then GoTo 1634 If RCZ > -ALIM Then GoTo 1647 RCZ = RCZ + Log(APHI) If IFORM = 2 Then RCZ = RCZ - 0.25 * Log(AARG) - AIC If RCZ > -ELIM Then GoTo 1638 1634 YR(NN) = ZEROR YI(NN) = ZEROI NN = NN - 1 NUF = NUF + 1 If NN = 0 Then GoTo 1647 GoTo 1626 1638 ASCLE = 1000 * D1MACH(1) / TOL ZLOG PHIR, PHII, STR, STI, IDUM CZR = CZR + STR CZI = CZI + STI If IFORM = 1 Then GoTo 1640 ZLOG ARGR, ARGI, STR, STI, IDUM CZR = CZR - 0.25 * STR - AIC CZI = CZI - 0.25 * STI 1640 AX = Exp(RCZ) / TOL AY = CZI CZR = AX * Cos(AY) CZI = AX * Sin(AY) ZUCHK CZR, CZI, NW, ASCLE, TOL If NW <> 0 Then GoTo 1634 GoTo 1647 1645 NUF = -1 1647 End Sub 'ZUOIK Sub ZS1S2(ZRR, ZRI, S1R, S1I, S2R, S2I, NZ, ASCLE, ALIM, IUF) '***BEGIN PROLOGUE ZS1S2 '***REFER TO ZBESK,ZAIRY ' ' ZS1S2 TESTS FOR A POSSIBLE UNDERFLOW RESULTING FROM THE ' ADDITION OF THE I AND K FUNCTIONS IN THE ANALYTIC CON- ' TINUATION FORMULA WHERE S1=K FUNCTION AND S2=I FUNCTION. ' ON KODE=1 THE I AND K FUNCTIONS ARE DIFFERENT ORDERS OF ' MAGNITUDE, BUT FOR KODE=2 THEY CAN BE OF THE SAME ORDER ' OF MAGNITUDE AND THE MAXIMUM MUST BE AT LEAST ONE ' PRECISION ABOVE THE UNDERFLOW LIMIT. ' '***ROUTINES CALLED ZABS,ZEXP,ZLOG '***END PROLOGUE ZS1S2 ' COMPLEX CZERO,C1,S1,S1D,S2,ZR 'Label: 1660 ' AA, ALN, AS1, AS2, C1I, C1R, S1DI, S1DR, ZEROI, ZEROR: Double ' IDUM: Integer ZEROR = 0#: ZEROI = 0# NZ = 0 AS1 = ZABS(S1R, S1I) AS2 = ZABS(S2R, S2I) If (S1R = 0#) And (S1I = 0#) Then GoTo 1660 If AS1 = 0# Then GoTo 1660 ALN = -ZRR - ZRR + Log(AS1) S1DR = S1R S1DI = S1I S1R = ZEROR S1I = ZEROI AS1 = ZEROR If ALN < -ALIM Then GoTo 1660 ZLOG S1DR, S1DI, C1R, C1I, IDUM C1R = C1R - ZRR - ZRR C1I = C1I - ZRI - ZRI ZEXP C1R, C1I, S1R, S1I AS1 = ZABS(S1R, S1I) IUF = IUF + 1 1660 AA = DMAX(AS1, AS2) If AA > ASCLE Then GoTo 1670 S1R = ZEROR S1I = ZEROI S2R = ZEROR S2I = ZEROI NZ = 1 IUF = 0 1670 End Sub 'ZS1S2 Sub ZSERI(ZR, ZI, FNU, KODE, N, YR(), YI(), NZ, TOL, ELIM, ALIM) '***BEGIN PROLOGUE ZSERI '***REFER TO ZBESI,ZBESK ' ' ZSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z) >= 0 BY ' MEANS OF THE POWER SERIES FOR LARGE CABS(Z) IN THE ' REGION CABS(Z) <= 2*SQR(FNU+1). NZ=0 IS A NORMAL RETURN. ' NZ > 0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO ' DUE TO UNDERFLOW. NZ < 0 MEANS UNDERFLOW OCCURRED, BUT THE ' CONDITION CABS(Z) <= 2*SQR(FNU+1) WAS VIOLATED AND THE ' COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ). ' '***ROUTINES CALLED DGAMLN,D1MACH,ZUCHK,ZABS,ZDIV,ZLOG,ZMLT '***END PROLOGUE ZSERI ' COMPLEX AK1,CK,COEF,CONE,CRSC,CSCL,CZ,CZERO,HZ,RZ,S1,S2,Y,Z 'Labels: 1702,1704,1706,1708,1710,1712,1714,1716,1718,1720,1722,1724,1726,1728,1730,1740 ' AA, ACZ, AK, AK1I, AK1R, ARM, ASCLE, ATOL,AZ, CKI, CKR, COEFI, ' COEFR, CONEI, CONER, CRSCR, CZI, CZR, DFNU, FNUP, HZI, HZR, ' RAZ, RS, RTR1, RZI, RZR, S, SS, STI, STR, S1I, S1R, S2I, S2R, ' ZEROI, ZEROR: Double ' I, IB, IDUM, IFLAG, IL, K, L, M, NN, NW: Integer Dim WR(10), WI(10) Dim STR As Double ZEROR = 0#: ZEROI = 0#: CONER = 1#: CONEI = 0# NZ = 0 AZ = ZABS(ZR, ZI) If AZ = 0# Then GoTo 1728 ARM = 1000 * D1MACH(1) RTR1 = Sqr(ARM) CRSCR = 1# IFLAG = 0 If AZ < ARM Then GoTo 1726 HZR = 0.5 * ZR HZI = 0.5 * ZI CZR = ZEROR CZI = ZEROI If AZ <= RTR1 Then GoTo 1702 ZMLT HZR, HZI, HZR, HZI, CZR, CZI 1702 ACZ = ZABS(CZR, CZI) NN = N ZLOG HZR, HZI, CKR, CKI, IDUM 1704 DFNU = FNU + 1# * (NN - 1) FNUP = DFNU + 1# '----------------------------------------------------------------------- ' UNDERFLOW TEST '----------------------------------------------------------------------- AK1R = CKR * DFNU AK1I = CKI * DFNU AK = DGAMLN(FNUP, IDUM) AK1R = AK1R - AK If KODE = 2 Then AK1R = AK1R - ZR If AK1R > -ELIM Then GoTo 1708 1706 NZ = NZ + 1 YR(NN) = ZEROR YI(NN) = ZEROI If ACZ > DFNU Then GoTo 1740 NN = NN - 1 If NN = 0 Then GoTo 1745 GoTo 1704 1708 If AK1R > -ALIM Then GoTo 1710 IFLAG = 1 SS = 1# / TOL CRSCR = TOL ASCLE = ARM * SS 1710 AA = Exp(AK1R) If IFLAG = 1 Then AA = AA * SS COEFR = AA * Cos(AK1I) COEFI = AA * Sin(AK1I) ATOL = TOL * ACZ / FNUP IL = IMIN(2, NN) For I = 1 To IL DFNU = FNU + 1# * (NN - I) FNUP = DFNU + 1# S1R = CONER S1I = CONEI If ACZ < TOL * FNUP Then GoTo 1714 AK1R = CONER AK1I = CONEI AK = FNUP + 2# S = FNUP AA = 2# 1712 RS = 1# / 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# AA = AA * ACZ * RS If AA > ATOL Then GoTo 1712 1714 S2R = S1R * COEFR - S1I * COEFI S2I = S1R * COEFI + S1I * COEFR WR(I) = S2R WI(I) = S2I If IFLAG = 0 Then GoTo 1716 ZUCHK S2R, S2I, NW, ASCLE, TOL If NW <> 0 Then GoTo 1706 1716 M = NN - I + 1 YR(M) = S2R * CRSCR YI(M) = S2I * CRSCR If I = IL Then GoTo 1718 ZDIV COEFR, COEFI, HZR, HZI, STR, STI COEFR = STR * DFNU COEFI = STI * DFNU 1718 Next I If NN <= 2 Then GoTo 1745 K = NN - 2 AK = 1# * K RAZ = 1# / AZ STR = ZR * RAZ STI = -ZI * RAZ RZR = (STR + STR) * RAZ RZI = (STI + STI) * RAZ If IFLAG = 1 Then GoTo 1722 ib = 3 1720 For I = ib To 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# K = K - 1 Next I GoTo 1745 '----------------------------------------------------------------------- ' RECUR BACKWARD WITH SCALED VALUES '----------------------------------------------------------------------- '----------------------------------------------------------------------- ' EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE ' UNDERFLOW LIMIT = ASCLE = D1MACH(1)*SS*1000. '----------------------------------------------------------------------- 1722 S1R = WR(1) S1I = WI(1) S2R = WR(2) S2I = WI(2) For L = 3 To 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# K = K - 1 If ZABS(CKR, CKI) > ASCLE Then GoTo 1724 Next L GoTo 1745 1724 ib = L + 1 If ib > NN Then GoTo 1745 GoTo 1720 1726 NZ = N If FNU = 0# Then NZ = NZ - 1 1728 YR(1) = ZEROR YI(1) = ZEROI If FNU <> 0# Then GoTo 1730 YR(1) = CONER YI(1) = CONEI 1730 If N = 1 Then GoTo 1745 For I = 2 To N YR(I) = ZEROR YI(I) = ZEROI Next I GoTo 1745 '----------------------------------------------------------------------- ' RETURN WITH NZ < 0 IF CABS(Z*Z/4) > FNU+N-NZ-1 COMPLETE ' THE CALCULATION IN CBINU WITH N=N-ABS(NZ) '----------------------------------------------------------------------- 1740 NZ = -NZ 1745 End Sub 'ZSERI Sub ZASYI(ZR, ZI, FNU, KODE, N, YR(), YI(), NZ, RL, TOL, ELIM, ALIM) '***BEGIN PROLOGUE ZASYI '***REFER TO ZBESI,ZBESK ' ZASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z) >= 0 BY ' MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE ' REGION CABS(Z) > MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN. ' NZ < 0 INDICATES AN OVERFLOW ON KODE=1. ' '***ROUTINES CALLED D1MACH,ZABS,ZDIV,ZEXP,ZMLT,ZSQRT '***END PROLOGUE ZASYI ' COMPLEX AK1,CK,CONE,CS1,CS2,CZ,CZERO,DK,EZ,P1,RZ,S2,Y,Z 'Labels: 10, 20, 30, 50, 60, 100, 110 ' 1752,1754,1756,1758,1760,1762,1770 ' AA, AEZ, AK, AK1I, AK1R, ARG, ARM, ATOL, AZ, BB, BK, CKI, CKR, ' CONEI, CONER, CS1I, CS1R, CS2I, CS2R, CZI, CZR, DFNU, DKI, ' DKR, DNU2, EZI, EZR, FDN, PI, P1I, P1R, RAZ, RTPI, RTR1, RZI, ' RZR, S, SGN0, SQK, STI, STR, S2I, S2R, TZI, TZR, ZEROI, ZEROR: Double ' I, IB, IL, INU, J, JL, K, KODED, M, NN: Integer Dim STR As Double RTPI = 0.159154943091895 ZEROR = 0#: ZEROI = 0#: CONER = 1#: CONEI = 0# PI = 4# * Atn(1#) NZ = 0 AZ = ZABS(ZR, ZI) ARM = 1000 * D1MACH(1) RTR1 = Sqr(ARM) IL = IMIN(2, N) DFNU = FNU + 1# * (N - IL) '----------------------------------------------------------------------- ' OVERFLOW TEST '----------------------------------------------------------------------- RAZ = 1# / AZ STR = ZR * RAZ STI = -ZI * RAZ AK1R = RTPI * STR * RAZ AK1I = RTPI * STI * RAZ ZSQRT AK1R, AK1I, AK1R, AK1I CZR = ZR CZI = ZI If KODE <> 2 Then GoTo 1752 CZR = ZEROR CZI = ZI 1752 If Abs(CZR) > ELIM Then GoTo 1762 DNU2 = DFNU + DFNU KODED = 1 If Abs(CZR) > ALIM And N > 2 Then GoTo 1754 KODED = 0 ZEXP CZR, CZI, STR, STI ZMLT AK1R, AK1I, STR, STI, AK1R, AK1I 1754 FDN = 0# If DNU2 > RTR1 Then FDN = DNU2 * DNU2 EZR = ZR * 8# EZI = ZI * 8# '----------------------------------------------------------------------- ' WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE ' FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE ' EXPANSION FOR THE IMAGINARY PART. '----------------------------------------------------------------------- AEZ = 8# * AZ S = TOL / AEZ JL = Int(RL + RL) + 2 P1R = ZEROR P1I = ZEROI If ZI = 0# Then GoTo 1756 '----------------------------------------------------------------------- ' CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF ' SIGNIFICANCE WHEN FNU OR N IS LARGE '----------------------------------------------------------------------- INU = Int(FNU) ARG = (FNU - 1# * INU) * PI INU = INU + N - IL AK = -Sin(ARG) BK = Cos(ARG) If ZI < 0# Then BK = -BK P1R = AK P1I = BK If (INU Mod 2) = 0 Then GoTo 1756 P1R = -P1R P1I = -P1I 1756 For K = 1 To IL SQK = FDN - 1# ATOL = S * Abs(SQK) SGN0 = 1# CS1R = CONER CS1I = CONEI CS2R = CONER CS2I = CONEI CKR = CONER CKI = CONEI AK = 0# AA = 1# BB = AEZ DKR = EZR DKI = EZI For J = 1 To JL ZDIV CKR, CKI, DKR, DKI, STR, STI CKR = STR * SQK CKI = STI * SQK CS2R = CS2R + CKR CS2I = CS2I + CKI SGN0 = -SGN0 CS1R = CS1R + CKR * SGN0 CS1I = CS1I + CKI * SGN0 DKR = DKR + EZR DKI = DKI + EZI AA = AA * Abs(SQK) / BB BB = BB + AEZ AK = AK + 8# SQK = SQK - AK If AA <= ATOL Then GoTo 1758 Next J GoTo 1770 1758 S2R = CS1R S2I = CS1I If ZR + ZR >= ELIM Then GoTo 1760 TZR = ZR + ZR TZI = ZI + ZI ZEXP -TZR, -TZI, STR, STI ZMLT STR, STI, P1R, P1I, STR, STI ZMLT STR, STI, CS2R, CS2I, STR, STI S2R = S2R + STR S2I = S2I + STI 1760 FDN = FDN + 8# * DFNU + 4# P1R = -P1R P1I = -P1I M = N - IL + K YR(M) = S2R * AK1R - S2I * AK1I YI(M) = S2R * AK1I + S2I * AK1R Next K If N <= 2 Then GoTo 1780 NN = N K = NN - 2 AK = 1# * K STR = ZR * RAZ STI = -ZI * RAZ RZR = (STR + STR) * RAZ RZI = (STI + STI) * RAZ ib = 3 For I = ib To 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# K = K - 1 If KODED = 0 Then GoTo 1780 Next I ZEXP CZR, CZI, CKR, CKI For I = 1 To NN STR = YR(I) * CKR - YI(I) * CKI YI(I) = YR(I) * CKI + YI(I) * CKR YR(I) = STR Next I GoTo 1780 1762 NZ = -1 GoTo 1780 1770 NZ = -2 1780 End Sub 'ZASYI Sub ZMLRI(ZR, ZI, FNU, KODE, N, YR(), YI(), NZ, TOL) '***BEGIN PROLOGUE ZMLRI '***REFER TO ZBESI,ZBESK ' ' ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z) >= 0 BY THE ' MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES. ' '***ROUTINES CALLED DGAMLN,D1MACH,ZABS,ZEXP,ZLOG,ZMLT '***END PROLOGUE ZMLRI ' COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z 'Labels: 1805, 1810, 1815, 1820, 1825, 1830, 1840 ' ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI, ' CNORMR, CONEI, CONER, FKAP, FKK, FLAM, FNF, PTI, PTR, P1I, ' P1R, P2I, P2R, RAZ, RHO, RHO2, RZI, RZR, SCLE, STI, STR, SUMI, ' SUMR, TFNF, TST, ZEROI, ZEROR: Double ' I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, M: Integer Dim STR As Double ZEROR = 0#: ZEROI = 0#: CONER = 1#: CONEI = 0# SCLE = D1MACH(1) / TOL NZ = 0 AZ = ZABS(ZR, ZI) IAZ = Int(AZ) IFNU = Int(FNU) INU = IFNU + N - 1 AT = 1# * IAZ + 1# RAZ = 1# / 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#) * RAZ RHO = ACK + Sqr(ACK * ACK - 1#) RHO2 = RHO * RHO TST = (RHO2 + RHO2) / ((RHO2 - 1#) * (RHO - 1#)) TST = TST / TOL '----------------------------------------------------------------------- ' COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES '----------------------------------------------------------------------- AK = AT For I = 1 To 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 > TST * AK * AK Then GoTo 1805 AK = AK + 1# Next I GoTo 1840 1805 I = I + 1 K = 0 If INU < IAZ Then GoTo 1815 '----------------------------------------------------------------------- ' COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS '----------------------------------------------------------------------- P1R = ZEROR P1I = ZEROI P2R = CONER P2I = CONEI AT = 1# * INU + 1# STR = ZR * RAZ STI = -ZI * RAZ CKR = STR * AT * RAZ CKI = STI * AT * RAZ ACK = AT * RAZ TST = Sqr(ACK / TOL) ITIME = 1 For K = 1 To 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 < TST Then GoTo 1810 If ITIME = 2 Then GoTo 1815 ACK = ZABS(CKR, CKI) FLAM = ACK + Sqr(ACK * ACK - 1#) FKAP = AP / ZABS(P1R, P1I) RHO = DMIN(FLAM, FKAP) TST = TST * Sqr(RHO / (RHO * RHO - 1#)) ITIME = 2 1810 Next K GoTo 1840 '----------------------------------------------------------------------- ' BACKWARD RECURRENCE AND SUM NORMALIZING RELATION '----------------------------------------------------------------------- 1815 K = K + 1 KK = IMAX(I + IAZ, K + INU) FKK = 1# * KK P1R = ZEROR P1I = ZEROI '----------------------------------------------------------------------- ' SCALE P2 AND SUM BY SCLE '----------------------------------------------------------------------- P2R = SCLE P2I = ZEROI FNF = FNU - 1# * IFNU TFNF = FNF + FNF BK = DGAMLN(FKK + TFNF + 1#, IDUM) - DGAMLN(FKK + 1#, IDUM) - DGAMLN(TFNF + 1#, IDUM) BK = Exp(BK) SUMR = ZEROR SUMI = ZEROI KM = KK - INU For I = 1 To 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# - TFNF / (FKK + TFNF) ACK = BK * AK SUMR = SUMR + (ACK + BK) * P1R SUMI = SUMI + (ACK + BK) * P1I BK = ACK FKK = FKK - 1# Next I YR(N) = P2R YI(N) = P2I If N = 1 Then GoTo 1825 For I = 2 To 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# - TFNF / (FKK + TFNF) ACK = BK * AK SUMR = SUMR + (ACK + BK) * P1R SUMI = SUMI + (ACK + BK) * P1I BK = ACK FKK = FKK - 1# M = N - I + 1 YR(M) = P2R YI(M) = P2I Next I 1825 If IFNU <= 0 Then GoTo 1830 For I = 1 To 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# - TFNF / (FKK + TFNF) ACK = BK * AK SUMR = SUMR + (ACK + BK) * P1R SUMI = SUMI + (ACK + BK) * P1I BK = ACK FKK = FKK - 1# Next I 1830 PTR = ZR PTI = ZI If KODE = 2 Then PTR = ZEROR ZLOG RZR, RZI, STR, STI, IDUM P1R = -FNF * STR + PTR P1I = -FNF * STI + PTI AP = DGAMLN(1# + FNF, IDUM) PTR = P1R - AP PTI = P1I '----------------------------------------------------------------------- ' THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW ' IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES '----------------------------------------------------------------------- P2R = P2R + SUMR P2I = P2I + SUMI AP = ZABS(P2R, P2I) P1R = 1# / AP ZEXP PTR, PTI, STR, STI CKR = STR * P1R CKI = STI * P1R PTR = P2R * P1R PTI = -P2I * P1R ZMLT CKR, CKI, PTR, PTI, CNORMR, CNORMI For I = 1 To N STR = YR(I) * CNORMR - YI(I) * CNORMI YI(I) = YR(I) * CNORMI + YI(I) * CNORMR YR(I) = STR Next I GoTo 1845 1840 NZ = -2 1845 End Sub 'ZMLRI Sub ZSHCH(ZR, ZI, CSHR, CSHI, CCHR, CCHI) '{***BEGIN PROLOGUE ZSHCH '***REFER TO ZBESK,ZBESH ' ' ZSHCH COMPUTES THE COMPLEX HYPERBOLIC FUNCTIONS CSH=SINH(X+I*Y) ' AND CCH=COSH(X+I*Y), WHERE I^2=-1. ' '***ROUTINES CALLED (NONE) '***END PROLOGUE ZSHCH ' CH, CN, SH, SN: Double SH = SINH(ZR) CH = COSH(ZR) SN = Sin(ZI) CN = Cos(ZI) CSHR = SH * CN CSHI = CH * SN CCHR = CH * CN CCHI = SH * SN End Sub Function COSH(ZR) COSH = (Exp(ZR) + Exp(-ZR)) / 2# End Function Function SINH(ZR) SINH = (Exp(ZR) - Exp(-ZR)) / 2# End Function Sub ZRATI(ZR, ZI, FNU, N, CYR(), CYI(), TOL) '***BEGIN PROLOGUE ZRATI '***REFER TO ZBESI,ZBESK,ZBESH ' ' ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD ' RECURRENCE. THE STARTING INDEX IS DETERMINED BY FORWARD ' RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B, ' MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973, ' BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER, ' BY D. J. SOOKNE. ' '***ROUTINES CALLED ZABS,ZDIV '***END PROLOGUE ZRATI ' COMPLEX Z,CY(1),CONE,CZERO,P1,P2,T1,RZ,PT,CDFNU 'Labels: 10, 20, 40, 50 ' 1910,1920,1940,1950 ' AK, AMAGZ, AP1, AP2, ARG, AZ, CDFNUI, CDFNUR, CONEI, CONER, ' CZEROI, CZEROR, DFNU, FDNU, FLAM, FNUP, PTI, PTR, P1I, P1R, ' P2I, P2R, RAK, RAP1, RHO, RT2, RZI, RZR, TEST, TEST1, TTI, ' TTR, T1I, T1R: Double ' I,ID, IDNU, INU, ITIME, K, KK, MAGZ: Integer CZEROR = 0#: CZEROI = 0#: CONER = 1#: CONEI = 0# RT2 = 1.4142135623731 AZ = ZABS(ZR, ZI) INU = Int(FNU) IDNU = INU + N - 1 MAGZ = Int(AZ) AMAGZ = 1# * (MAGZ + 1) FDNU = 1# * IDNU FNUP = DMAX(AMAGZ, FDNU) ID = IDNU - MAGZ - 1 ITIME = 1 K = 1 PTR = 1# / 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 > 0 Then ID = 0 AP2 = ZABS(P2R, P2I) AP1 = ZABS(P1R, P1I) '----------------------------------------------------------------------- ' THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU ' GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT ' P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR ' PREMATURELY. '----------------------------------------------------------------------- ARG = (AP2 + AP2) / (AP1 * TOL) TEST1 = Sqr(ARG) TEST = TEST1 RAP1 = 1# / AP1 P1R = P1R * RAP1 P1I = P1I * RAP1 P2R = P2R * RAP1 P2I = P2I * RAP1 AP2 = AP2 * RAP1 1910 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 <= TEST Then GoTo 1910 If ITIME = 2 Then GoTo 1920 AK = ZABS(T1R, T1I) * 0.5 FLAM = AK + Sqr(AK * AK - 1#) RHO = DMIN(AP2 / AP1, FLAM) TEST = TEST1 * Sqr(RHO / (RHO * RHO - 1#)) ITIME = 2 GoTo 1910 1920 KK = K + 1 - ID AK = 1# * KK T1R = AK T1I = CZEROI DFNU = FNU + 1# * (N - 1) P1R = 1# / AP2 P1I = CZEROI P2R = CZEROR P2I = CZEROI For I = 1 To 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 Next I If P1R <> CZEROR Or P1I <> CZEROI Then GoTo 1940 P1R = TOL P1I = TOL 1940 ZDIV P2R, P2I, P1R, P1I, CYR(N), CYI(N) If N = 1 Then GoTo 1955 K = N - 1 AK = 1# * K T1R = AK T1I = CZEROI CDFNUR = FNU * RZR CDFNUI = FNU * RZI For I = 2 To 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 <> CZEROR Then GoTo 1950 PTR = TOL PTI = TOL AK = TOL * RT2 1950 RAK = CONER / AK CYR(K) = RAK * PTR * RAK CYI(K) = -RAK * PTI * RAK T1R = T1R - CONER K = K - 1 Next I 1955 End Sub 'ZRATI 'end of file Cbess0.bas