Attribute VB_Name = "Module4" '************************************************************************** '* Procedures 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/01/2005). * '* (www.jpmoreau.fr) * '************************************************************************** 'UNIT CBESS00 (subroutines ZAIRY, ZBKNU) DefDbl A-H, O-Z DefInt I-N 'LIST OF SUBROUTINES: ' Subroutine ZAIRY(ZR, ZI, ID, KODE, AIR, AII, NZ, IERR) ' Subroutine ZBKNU(ZR, ZI, FNU, KODE, N, NZ, TOL, ELIM, ALIM) ' Subroutine ZACAI(ZR, ZI, FNU, KODE, MR, N, NZ, RL, TOL, ELIM, ALIM) ' Subroutine ZKSCL(ZRR,ZRI,FNU, N, NZ, RZR,RZI,ASCLE,TOL,ELIM) Sub ZAIRY(ZR, ZI, ID, KODE, AIR, AII, NZ, IERR) '***BEGIN PROLOGUE ZAIRY '***DATE WRITTEN 830501 (YYMMDD) '***REVISION DATE 830501 (YYMMDD) '***CATEGORY NO. B5K '***KEYWORDS AIRY FUNCTION,BESSEL FUNCTIONS OF ORDER ONE THIRD '***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES '***PURPOSE TO COMPUTE AIRY FUNCTIONS AI(Z) AND DAI(Z) FOR COMPLEX Z '***DESCRIPTION ' ' ***A DOUBLE PRECISION ROUTINE*** ' ON KODE=1, ZAIRY COMPUTES THE COMPLEX AIRY FUNCTION AI(Z) OR ' ITS DERIVATIVE DAI(Z)/DZ ON ID=0 OR ID=1 RESPECTIVELY. ON ' KODE=2, A SCALING OPTION CEXP(ZTA)*AI(Z) OR CEXP(ZTA)* ' DAI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL DECAY IN ' -PI/3 < ARG(Z) < PI/3 AND THE EXPONENTIAL GROWTH IN ' PI/3 < ABS(ARG(Z)) < PI, WHERE ZTA=(2/3)*Z*CSQRT(Z). ' ' WHILE THE AIRY FUNCTIONS AI(Z) AND DAI(Z)/DZ ARE ANALYTIC IN ' THE WHOLE Z PLANE, THE CORRESPONDING SCALED FUNCTIONS DEFINED ' FOR KODE=2 HAVE A CUT ALONG THE NEGATIVE REAL AXIS. ' DEFINTIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF ' MATHEMATICAL FUNCTIONS (REF. 1). ' ' INPUT ZR,ZI ARE DOUBLE PRECISION ' ZR,ZI - Z=CMPLX(ZR,ZI) ' ID - ORDER OF DERIVATIVE, ID=0 OR ID=1 ' KODE - A PARAMETER TO INDICATE THE SCALING OPTION ' KODE= 1 RETURNS ' AI=AI(Z) ON ID=0 OR ' AI=DAI(Z)/DZ ON ID=1 ' = 2 RETURNS ' AI=CEXP(ZTA)*AI(Z) ON ID=0 OR ' AI=CEXP(ZTA)*DAI(Z)/DZ ON ID=1 WHERE 'ZTA = (2 / 3) * Z * CSQRT(Z) ' ' OUTPUT AIR,AII ARE DOUBLE PRECISION ' AIR,AII- COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND ' KODE ' NZ - UNDERFLOW INDICATOR ' NZ= 0 , NORMAL RETURN ' NZ= 1 , AI=CMPLX(0.0,0.0) DUE TO UNDERFLOW IN ' -PI/3 < ARG(Z) < PI/3 ON KODE=1. ' IERR - ERROR FLAG ' IERR=0, NORMAL RETURN - COMPUTATION COMPLETED ' IERR=1, INPUT ERROR - NO COMPUTATION ' IERR=2, OVERFLOW - NO COMPUTATION, REAL(ZTA) ' TOO LARGE ON KODE=1 ' IERR=3, CABS(Z) LARGE - COMPUTATION COMPLETED ' LOSSES OF SIGNIFCANCE BY ARGUMENT REDUCTION ' PRODUCE LESS THAN HALF OF MACHINE ACCURACY ' IERR=4, CABS(Z) TOO LARGE - NO COMPUTATION ' COMPLETE LOSS OF ACCURACY BY ARGUMENT ' REDUCTION ' IERR=5, ERROR - NO COMPUTATION, ' ALGORITHM TERMINATION CONDITION NOT MET ' '***LONG DESCRIPTION ' ' AI AND DAI ARE COMPUTED FOR CABS(Z).GT.1.0 FROM THE K BESSEL ' FUNCTIONS BY ' ' AI(Z)=C*SQRT(Z)*K(1/3,ZTA) , DAI(Z)=-C*Z*K(2/3,ZTA) ' ' with C = 1# / (PI * SQRT(3#)) ' ZTA = (2 / 3) * Z ^ (3 / 2) ' ' WITH THE POWER SERIES FOR CABS(Z) <= 1. ' ' IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- ' MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z IS LARGE, LOSSES ' OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. CONSEQUENTLY, IF ' THE MAGNITUDE OF ZETA=(2/3)*Z^1.5 EXCEEDS U1=SQR(0.5/UR), ' THEN LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR ' FLAG IERR=3 IS TRIGGERED WHERE UR=DMAX(D1MACH(4),1.0D-18) IS ' DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. ' ALSO, IF THE MAGNITUDE OF ZETA IS LARGER THAN U2=0.5/UR, THEN ' ALL SIGNIFICANCE IS LOST AND IERR=4. IN ORDER TO USE THE INT ' FUNCTION, ZETA MUST BE FURTHER RESTRICTED NOT TO EXCEED THE ' LARGEST INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF ZETA ' MUST BE RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, ' AND U3 ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE ' PRECISION ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE ' PRECISION ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMIT- ' ING IN THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT THE MAG- ' NITUDE OF Z CANNOT EXCEED 3.1E+4 IN SINGLE AND 2.1E+6 IN ' DOUBLE PRECISION ARITHMETIC. THIS ALSO MEANS THAT ONE CAN ' EXPECT TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, ' NO DIGITS IN SINGLE PRECISION AND ONLY 7 DIGITS IN DOUBLE ' PRECISION ARITHMETIC. SIMILAR CONSIDERATIONS HOLD FOR OTHER ' MACHINES. ' ' THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX ' BESSEL FUNCTION CAN BE EXPRESSED BY P*10^S WHERE P=MAX(UNIT ' ROUNDOFF,1E-18) IS THE NOMINAL PRECISION AND 10^S REPRE- ' SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE ' ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), ' ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF ' CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY ' HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN ' ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY ' SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10^K LARGER ' THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, ' 0) SIGNIFICANT DIGITS OR, STATED ANOTHER WAY, WHEN K EXCEEDS ' THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER ' COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY ' BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER ' COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE ' MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, ' THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, ' OR -PI/2+P. ' '***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ ' AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF ' COMMERCE, 1955. ' ' COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT ' AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 ' ' A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX ' ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- ' 1018, MAY, 1985 ' ' A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX ' ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. 'Math.SOFTWARE , 1986 ' '***ROUTINES CALLED ZACAI,ZBKNU,ZEXP,ZSQRT,I1MACH,D1MACH '***END PROLOGUE ZAIRY ' COMPLEX AI,CONE,CSQ,CY,S1,S2,TRM1,TRM2,Z,ZTA,Z3 } 'Labels: 40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,190,200,210,260,270,280,300 ' AA, AD, AK, ALIM, ATRM, AZ, AZ3, BK, ' CC, CK, COEF, CONEI, CONER, CSQI, CSQR, 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, ZTAI, ZTAR, Z3I, Z3R, ALAZ, BB: Double ' IFLAG, K, K1, K2, MR, NN: Integer Dim STR As Double Dim CYR(10), CYI(10) TTH = 0.666666666666667 C1 = 0.355028053887817 C2 = 0.258819403792807 COEF = 0.183776298473931 ZEROR = 0#: ZEROI = 0#: CONER = 1#: CONEI = 0# IERR = 0 NZ = 0 If (ID < 0) Or (ID > 1) Then IERR = 1 If (KODE < 1) Or (KODE > 2) Then IERR = 1 If IERR <> 0 Then GoTo 300 AZ = ZABS(ZR, ZI) TOL = DMAX(D1MACH(4), 1E-18) FID = 1# * ID If AZ > 1# Then GoTo 70 '----------------------------------------------------------------------- ' POWER SERIES FOR CABS(Z) <= 1. '----------------------------------------------------------------------- S1R = CONER S1I = CONEI S2R = CONER S2I = CONEI If AZ < TOL Then GoTo 170 AA = AZ * AZ If AA < TOL / AZ Then GoTo 40 TRM1R = CONER TRM1I = CONEI TRM2R = CONER TRM2I = CONEI ATRM = 1# 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# + FID BK = 3# - FID - FID CK = 4# - FID DK = 3# + FID + FID D1 = AK * DK D2 = BK * CK AD = DMIN(D1, D2) AK = 24# + 9# * FID BK = 30# - 9# * FID For K = 1 To 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 = DMIN(D1, D2) If ATRM < TOL * AD Then GoTo 40 AK = AK + 18# BK = BK + 18# Next K 40 If ID = 1 Then GoTo 50 AIR = S1R * C1 - C2 * (ZR * S2R - ZI * S2I) AII = S1I * C1 - C2 * (ZR * S2I + ZI * S2R) If KODE = 1 Then GoTo 300 ZSQRT ZR, ZI, STR, STI ZTAR = TTH * (ZR * STR - ZI * STI) ZTAI = TTH * (ZR * STI + ZI * STR) ZEXP ZTAR, ZTAI, STR, STI PTR = AIR * STR - AII * STI AII = AIR * STI + AII * STR AIR = PTR GoTo 300 50 AIR = -S2R * C2 AII = -S2I * C2 If AZ <= TOL Then GoTo 60 STR = ZR * S1R - ZI * S1I STI = ZR * S1I + ZI * S1R CC = C1 / (1# + FID) AIR = AIR + CC * (STR * ZR - STI * ZI) AII = AII + CC * (STR * ZI + STI * ZR) 60 If KODE = 1 Then GoTo 300 ZSQRT ZR, ZI, STR, STI ZTAR = TTH * (ZR * STR - ZI * STI) ZTAI = TTH * (ZR * STI + ZI * STR) ZEXP ZTAR, ZTAI, STR, STI PTR = STR * AIR - STI * AII AII = STR * AII + STI * AIR AIR = PTR GoTo 300 '----------------------------------------------------------------------- ' CASE FOR CABS(Z) > 1 '----------------------------------------------------------------------- 70 FNU = (1# + FID) / 3# '----------------------------------------------------------------------- ' SET PARAMETERS RELATED TO MACHINE CONSTANTS. ' TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0D-18. ' ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. ' EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND ' EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR ' UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. ' RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. ' DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10^(-DIG). '----------------------------------------------------------------------- K1 = Int(XI1MACH(15)) K2 = Int(XI1MACH(16)) R1M5 = D1MACH(5) K = IMIN(Abs(K1), Abs(K2)) ELIM = 2.303 * (K * R1M5 - 3#) K1 = Int(XI1MACH(14)) - 1 AA = R1M5 * K1 DIG = DMIN(AA, 18#) AA = AA * 2.303 ALIM = ELIM + DMAX(-AA, -41.45) RL = 1.2 * DIG + 3# ALAZ = Log(AZ) '----------------------------------------------------------------------- ' TEST FOR PROPER RANGE '----------------------------------------------------------------------- AA = 0.5 / TOL BB = 0.5 * Int(XI1MACH(9)) AA = DMIN(AA, BB) AA = AA ^ TTH If AZ > AA Then GoTo 260 AA = Sqr(AA) If AZ > AA Then IERR = 3 ZSQRT ZR, ZI, CSQR, CSQI ZTAR = TTH * (ZR * CSQR - ZI * CSQI) ZTAI = TTH * (ZR * CSQI + ZI * CSQR) '----------------------------------------------------------------------- ' RE(ZTA) <= 0 WHEN RE(Z) < 0, ESPECIALLY WHEN IM(Z) IS SMALL '----------------------------------------------------------------------- IFLAG = 0 SFAC = 1# AK = ZTAI If ZR >= 0# Then GoTo 80 BK = ZTAR CK = -Abs(BK) ZTAR = CK ZTAI = AK 80 If ZI <> 0# Then GoTo 90 If ZR > 0# Then GoTo 90 ZTAR = 0# ZTAI = AK 90 AA = ZTAR If (AA >= 0#) And (ZR > 0#) Then GoTo 110 If KODE = 2 Then GoTo 100 '----------------------------------------------------------------------- ' OVERFLOW TEST '----------------------------------------------------------------------- If AA > -ALIM Then GoTo 100 AA = -AA + 0.25 * ALAZ IFLAG = 1 SFAC = TOL If AA > ELIM Then GoTo 270 '----------------------------------------------------------------------- ' CBKNU AND CACON RETURN EXP(ZTA)*K(FNU,ZTA) ON KODE=2 '----------------------------------------------------------------------- 100 MR = 1 If ZI < 0# Then MR = -1 ZACAI ZTAR, ZTAI, FNU, KODE, MR, 1, CYR(), CYI(), NN, RL, TOL, ELIM, ALIM If NN < 0 Then GoTo 280 NZ = NZ + NN GoTo 130 110 If KODE = 2 Then GoTo 120 '----------------------------------------------------------------------- ' UNDERFLOW TEST '----------------------------------------------------------------------- If AA < ALIM Then GoTo 120 AA = -AA - 0.25 * ALAZ IFLAG = 2 SFAC = 1# / TOL If AA < -ELIM Then GoTo 210 120 ZBKNU ZTAR, ZTAI, FNU, KODE, 1, CYR(), CYI(), NZ, TOL, ELIM, ALIM 130 S1R = CYR(1) * COEF S1I = CYI(1) * COEF If IFLAG <> 0 Then GoTo 150 If ID = 1 Then GoTo 140 AIR = CSQR * S1R - CSQI * S1I AII = CSQR * S1I + CSQI * S1R GoTo 300 140 AIR = -(ZR * S1R - ZI * S1I) AII = -(ZR * S1I + ZI * S1R) GoTo 300 150 S1R = S1R * SFAC S1I = S1I * SFAC If ID = 1 Then GoTo 160 STR = S1R * CSQR - S1I * CSQI S1I = S1R * CSQI + S1I * CSQR S1R = STR AIR = S1R / SFAC AII = S1I / SFAC GoTo 300 160 STR = -(S1R * ZR - S1I * ZI) S1I = -(S1R * ZI + S1I * ZR) S1R = STR AIR = S1R / SFAC AII = S1I / SFAC GoTo 300 170 AA = 1000 * D1MACH(1) S1R = ZEROR S1I = ZEROI If ID = 1 Then GoTo 190 If AZ <= AA Then GoTo 180 S1R = C2 * ZR S1I = C2 * ZI 180 AIR = C1 - S1R AII = -S1I GoTo 300 190 AIR = -C2 AII = 0# AA = Sqr(AA) If AZ <= AA Then GoTo 200 S1R = 0.5 * (ZR * ZR - ZI * ZI) S1I = ZR * ZI 200 AIR = AIR + C1 * S1R AII = AII + C1 * S1I GoTo 300 210 NZ = 1 AIR = ZEROR AII = ZEROI GoTo 300 270 NZ = 0 IERR = 2 GoTo 300 280 If NN = -1 Then GoTo 270 NZ = 0 IERR = 5 GoTo 300 260 IERR = 4 NZ = 0 300 End Sub 'ZAIRY Sub ZACAI(ZR, ZI, FNU, KODE, MR, N, YR(), YI(), NZ, RL, TOL, ELIM, ALIM) '***BEGIN PROLOGUE ZACAI '***REFER TO ZAIRY ' ' ZACAI APPLIES THE ANALYTIC CONTINUATION FORMULA ' 'K(FNU, ZN * Exp(MP)) = K(FNU, ZN) * Exp(-MP * FNU) - MP * I(FNU, ZN) 'MP = PI * MR * CMPLX(0#, 1#) ' ' TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT ' HALF Z PLANE FOR USE WITH ZAIRY WHERE FNU=1/3 OR 2/3 AND N=1. ' ZACAI IS THE SAME AS ZACON WITH THE PARTS FOR LARGER ORDERS AND ' RECURRENCE REMOVED. A RECURSIVE CALL TO ZACON CAN RESULT IF ZACON ' IS CALLED FROM ZAIRY. ' '***ROUTINES CALLED ZASYI,ZBKNU,ZMLRI,ZSERI,ZS1S2,D1MACH,ZABS '***END PROLOGUE ZACAI ' COMPLEX CSGN,CSPN,C1,C2,Y,Z,ZN,Y 'Labels: 10,20,30,40,50,60,70,80, 100 ' ARG, ASCLE, AZ, CSGNR, CSGNI, CSPNR, CSPNI, C1R, C1I, C2R, C2I, ' DFNU, FMR, SGN, YY, ZNR, ZNI: Double ' INU, IUF, NN, NW: Integer Dim PI As Double PI = 4# * Atn(1#) NZ = 0 ZNR = -ZR ZNI = -ZI AZ = ZABS(ZR, ZI) NN = N DFNU = FNU + 1# * (N - 1) If AZ <= 2# Then GoTo 10 If AZ * AZ * 0.25 > DFNU + 1# Then GoTo 20 '----------------------------------------------------------------------- ' POWER SERIES FOR THE I FUNCTION '----------------------------------------------------------------------- 10 ZSERI ZNR, ZNI, FNU, KODE, NN, YR(), YI(), NW, TOL, ELIM, ALIM GoTo 40 20 If AZ < RL Then GoTo 30 '----------------------------------------------------------------------- ' ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION '----------------------------------------------------------------------- ZASYI ZNR, ZNI, FNU, KODE, NN, YR(), YI(), NW, RL, TOL, ELIM, ALIM If NW < 0 Then GoTo 80 GoTo 40 '----------------------------------------------------------------------- ' MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION '----------------------------------------------------------------------- 30 ZMLRI ZNR, ZNI, FNU, KODE, NN, YR(), YI(), NW, TOL If NW < 0 Then GoTo 80 '----------------------------------------------------------------------- ' ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION '----------------------------------------------------------------------- 40 ZBKNU ZNR, ZNI, FNU, KODE, 1, YR(), YI(), NW, TOL, ELIM, ALIM If NW <> 0 Then GoTo 80 FMR = 1# * MR SGN0 = -Sign(PI, FMR) CSGNR = 0# CSGNI = SGN0 If KODE = 1 Then GoTo 50 YY = -ZNI CSGNR = -CSGNI * Sin(YY) CSGNI = CSGNI * Cos(YY) '----------------------------------------------------------------------- ' CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE ' WHEN FNU IS LARGE '----------------------------------------------------------------------- 50: INU = Int(FNU) ARG = (FNU - 1# * INU) * SGN0 CSPNR = Cos(ARG) CSPNI = Sin(ARG) If (INU Mod 2) = 0 Then GoTo 60 CSPNR = -CSPNR CSPNI = -CSPNI 60 C1R = YR(1) C1I = YI(1) C2R = YR(1) C2I = YI(1) If KODE = 1 Then GoTo 70 IUF = 0 ASCLE = 1000 * D1MACH(1) / TOL ZS1S2 ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF NZ = NZ + NW 70 YR(1) = CSPNR * C1R - CSPNI * C1I + CSGNR * C2R - CSGNI * C2I YI(1) = CSPNR * C1I + CSPNI * C1R + CSGNR * C2I + CSGNI * C2R GoTo 100 80 NZ = -1 If NW = -2 Then NZ = -2 100 End Sub 'ZACAI Sub ZBKNU(ZR, ZI, FNU, KODE, N, YR(), YI(), NZ, TOL, ELIM, ALIM) '***BEGIN PROLOGUE ZBKNU '***REFER TO ZBESI,ZBESK,ZAIRY,ZBESH ' ' ZBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE. ' '***ROUTINES CALLED DGAMLN,I1MACH,D1MACH,ZKSCL,ZSHCH,ZUCHK,ZABS,ZDIV, ' ZEXP,ZLOG,ZMLT,ZSQRT '***END PROLOGUE ZBKNU } 'Labels: 10,30,40,50,60,70,80,90,100,110,120,130,140,160,170,180,200,210,215,220, ' 225,230,240,250,260,261,262,263,264, 270,280,290,300,310, 400 ' AA, AK, ASCLE, A1, A2, BB, BK, CAZ, ' CBI, CBR, CCHI, CCHR, CKI, CKR, COEFI, COEFR, CONEI, CONER, ' CRSCR, CSCLR, CSHI, CSHR, CSI, CSR, CTWOI, CTWOR, ' CZEROI, CZEROR, CZI, CZR, DNU, DNU2, DPI, ETEST, FC, FHS, ' FI, FK, FKS, FMUI, FMUR, 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, ' TTH, T1, T2, ELM, CELMR, ZDR, ZDI, AS, ALAS, HELIM: Double ' I, IFLAG, INU, K, KFLAG, KK, KMAX, KODED, IDUM, J, IC, INUB, NW: Integer Dim CC(8), CSSR(3), CSRR(3), BRY(3) Dim STR As Double ' COMPLEX: Z,Y,A,B,RZ,SMU,FU,FMU,F,FLRZ,CZ,S1,S2,CSH,CCH, ' CK,P,Q,COEF,P1,P2,CBK,PT,CZERO,CONE,CTWO,ST,EZ,CS,DK KMAX = 30 CZEROR = 0#: CZEROI = 0#: CONER = 1#: CONEI = 0# CTWOR = 2#: CTWOI = 0#: R1 = 2# DPI = 3.14159265358979: RTHPI = 1.2533141373155 SPI = 1.90985931710274: HPI = 1.5707963267949 FPI = 1.89769999331518: TTH = 0.666666666666667 CC(1) = 0.577215664901533: CC(2) = -4.20026350340952E-02 CC(3) = -4.21977345555443E-02: CC(4) = 7.2189432466631E-03 CC(5) = -2.15241674114951E-04: CC(6) = -2.01348547807882E-05 CC(7) = 1.1330272319817E-06: CC(8) = 6.11609510448142E-09 CAZ = ZABS(ZR, ZI) CSCLR = 1# / TOL CRSCR = TOL CSSR(1) = CSCLR CSSR(2) = 1# CSSR(3) = CRSCR CSRR(1) = CRSCR CSRR(2) = 1# CSRR(3) = CSCLR BRY(1) = 1000 * D1MACH(1) / TOL BRY(2) = 1# / BRY(1) BRY(3) = D1MACH(2) NZ = 0 IFLAG = 0 KODED = KODE RCAZ = 1# / CAZ STR = ZR * RCAZ STI = -ZI * RCAZ RZR = (STR + STR) * RCAZ RZI = (STI + STI) * RCAZ INU = Int(FNU + 0.5) DNU = FNU - 1# * INU If Abs(DNU) = 0.5 Then GoTo 110 DNU2 = 0# If Abs(DNU) > TOL Then DNU2 = DNU * DNU If CAZ > R1 Then GoTo 110 '----------------------------------------------------------------------- ' SERIES FOR CABS(Z) <= R1 '----------------------------------------------------------------------- FC = 1# ZLOG RZR, RZI, SMUR, SMUI, IDUM FMUR = SMUR * DNU FMUI = SMUI * DNU ZSHCH FMUR, FMUI, CSHR, CSHI, CCHR, CCHI If DNU = 0# Then GoTo 10 FC = DNU * DPI FC = FC / Sin(FC) SMUR = CSHR / DNU SMUI = CSHI / DNU 10 A2 = 1# + DNU '----------------------------------------------------------------------- ' GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU) '----------------------------------------------------------------------- T2 = Exp(-DGAMLN(A2, IDUM)) T1 = 1# / (T2 * FC) If Abs(DNU) > 0.1 Then GoTo 40 '----------------------------------------------------------------------- ' SERIES FOR F0 TO RESOLVE INDETERMINAY FOR SMALL ABS(DNU) '----------------------------------------------------------------------- AK = 1# S = CC(1) For K = 2 To 8 AK = AK * DNU2 TM = CC(K) * AK S = S + TM If Abs(TM) < TOL Then GoTo 30 Next K 30 G1 = -S GoTo 50 40 G1 = (T1 - T2) / (DNU + DNU) 50 G2 = (T1 + T2) * 0.5 FR = FC * (CCHR * G1 + SMUR * G2) FI = FC * (CCHI * G1 + SMUI * G2) ZEXP FMUR, FMUI, STR, STI PR = 0.5 * STR / T2 PI = 0.5 * STI / T2 ZDIV 0.5, 0#, STR, STI, PTR, PTI QR = PTR / T1 QI = PTI / T1 S1R = FR S1I = FI S2R = PR S2I = PI AK = 1# A1 = 1# CKR = CONER CKI = CONEI BK = 1# - DNU2 If (INU > 0) Or (N > 1) Then GoTo 80 '----------------------------------------------------------------------- ' GENERATE K(FNU,Z), 0.0 <= FNU < 0.5 AND N=1 '----------------------------------------------------------------------- If CAZ < TOL Then GoTo 70 ZMLT ZR, ZI, ZR, ZI, CZR, CZI CZR = 0.25 * CZR CZI = 0.25 * CZI T1 = 0.25 * CAZ * CAZ 60 FR = (FR * AK + PR + QR) / BK FI = (FI * AK + PI + QI) / BK STR = 1# / (AK - DNU) PR = PR * STR PI = PI * STR STR = 1# / (AK + DNU) QR = QR * STR QI = QI * STR STR = CKR * CZR - CKI * CZI RAK = 1# / 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# AK = AK + 1# If A1 > TOL Then GoTo 60 70 YR(1) = S1R YI(1) = S1I If KODED = 1 Then GoTo 400 ZEXP ZR, ZI, STR, STI ZMLT S1R, S1I, STR, STI, YR(1), YI(1) GoTo 400 '----------------------------------------------------------------------- ' GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE '----------------------------------------------------------------------- 80 If CAZ < TOL Then GoTo 100 ZMLT ZR, ZI, ZR, ZI, CZR, CZI CZR = 0.25 * CZR CZI = 0.25 * CZI T1 = 0.25 * CAZ * CAZ 90 FR = (FR * AK + PR + QR) / BK FI = (FI * AK + PI + QI) / BK STR = 1# / (AK - DNU) PR = PR * STR PI = PI * STR STR = 1# / (AK + DNU) QR = QR * STR QI = QI * STR STR = CKR * CZR - CKI * CZI RAK = 1# / 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# AK = AK + 1# If A1 > TOL Then GoTo 90 100: KFLAG = 2 A1 = FNU + 1# AK = A1 * Abs(SMUR) If AK > ALIM Then KFLAG = 3 STR = CSSR(KFLAG) P2R = S2R * STR P2I = S2I * STR ZMLT P2R, P2I, RZR, RZI, S2R, S2I S1R = S1R * STR S1I = S1I * STR If KODED = 1 Then GoTo 210 ZEXP ZR, ZI, FR, FI ZMLT S1R, S1I, FR, FI, S1R, S1I ZMLT S2R, S2I, FR, FI, S2R, S2I GoTo 210 '----------------------------------------------------------------------- ' IFLAG=0 MEANS NO UNDERFLOW OCCURRED ' IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH ' KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD ' RECURSION '----------------------------------------------------------------------- 110 ZSQRT ZR, ZI, STR, STI ZDIV RTHPI, CZEROI, STR, STI, COEFR, COEFI KFLAG = 2 If KODED = 2 Then GoTo 120 If ZR > ALIM Then GoTo 290 STR = Exp(-ZR) * CSSR(KFLAG) STI = -STR * Sin(ZI) STR = STR * Cos(ZI) ZMLT COEFR, COEFI, STR, STI, COEFR, COEFI 120 If Abs(DNU) = 0.5 Then GoTo 300 '----------------------------------------------------------------------- ' MILLER ALGORITHM FOR CABS(Z) > R1 '----------------------------------------------------------------------- AK = Cos(DPI * DNU) AK = Abs(AK) If AK = CZEROR Then GoTo 300 FHS = Abs(0.25 - DNU2) If FHS = CZEROR Then GoTo 300 '----------------------------------------------------------------------- ' COMPUTE R2=F(E). IF CABS(Z) >= R2, USE FORWARD RECURRENCE TO ' DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON ' 12.LE.E.LE.60. E IS COMPUTED FROM 2^(-E)=B^(1-XI1MACH(14))= ' TOL, WHERE B IS THE BASE OF THE ARITHMETIC. '----------------------------------------------------------------------- T1 = 1# * (XI1MACH(14) - 1) T1 = T1 * D1MACH(5) * 3.321928094 T1 = DMAX(T1, 12#) T1 = DMIN(T1, 60#) T2 = TTH * T1 - 6# If ZR <> 0# Then GoTo 130 T1 = HPI GoTo 140 130 T1 = Atn(ZI / ZR) T1 = Abs(T1) 140 If T2 > CAZ Then GoTo 170 '----------------------------------------------------------------------- ' FORWARD RECURRENCE LOOP WHEN CABS(Z) >= R2 '----------------------------------------------------------------------- ETEST = AK / (DPI * CAZ * TOL) FK = CONER If ETEST < CONER Then GoTo 180 FKS = CTWOR CKR = CAZ + CAZ + CTWOR P1R = CZEROR P2R = CONER For I = 1 To 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 < STR Then GoTo 160 Next I GoTo 310 160 FK = FK + SPI * T1 * Sqr(T2 / CAZ) FHS = Abs(0.25 - DNU2) GoTo 180 '----------------------------------------------------------------------- ' COMPUTE BACKWARD INDEX K FOR CABS(Z) < R2 '----------------------------------------------------------------------- 170 A2 = Sqr(CAZ) AK = FPI * AK / (TOL * Sqr(A2)) AA = 3# * T1 / (1# + CAZ) BB = 14.7 * T1 / (28# + CAZ) AK = (Log(AK) + CAZ * Cos(AA) / (1# + 0.008 * CAZ)) / Cos(BB) FK = 0.12125 * AK * AK / CAZ + 1.5 '----------------------------------------------------------------------- ' BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM '----------------------------------------------------------------------- 180 K = Int(FK) FK = 1# * K FKS = FK * FK P1R = CZEROR P1I = CZEROI P2R = TOL P2I = CZEROI CSR = P2R CSI = P2I For I = 1 To K A1 = FKS - FK AK = (FKS + FK) / (A1 + FHS) RAK = 2# / (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 Next I '----------------------------------------------------------------------- ' COMPUTE (P2/CS)=(P2/CABS(CS))*(CONJG(CS)/CABS(CS)) FOR BETTER ' SCALING '----------------------------------------------------------------------- TM = ZABS(CSR, CSI) PTR = 1# / TM S1R = P2R * PTR S1I = P2I * PTR CSR = CSR * PTR CSI = -CSI * PTR ZMLT COEFR, COEFI, S1R, S1I, STR, STI ZMLT STR, STI, CSR, CSI, S1R, S1I If (INU > 0) Or (N > 1) Then GoTo 200 ZDR = ZR ZDI = ZI If IFLAG = 1 Then GoTo 270 GoTo 240 '----------------------------------------------------------------------- ' COMPUTE P1/P2=(P1/CABS(P2)*CONJG(P2)/CABS(P2) FOR SCALING '----------------------------------------------------------------------- 200 TM = ZABS(P2R, P2I) PTR = 1# / TM P1R = P1R * PTR P1I = P1I * PTR P2R = P2R * PTR P2I = -P2I * PTR ZMLT P1R, P1I, P2R, P2I, PTR, PTI STR = DNU + 0.5 - PTR STI = -PTI ZDIV STR, STI, ZR, ZI, STR, STI STR = STR + 1# ZMLT STR, STI, S1R, S1I, S2R, S2I '----------------------------------------------------------------------- ' FORWARD RECURSION ON THE THREE TERM RECURSION WITH RELATION WITH ' SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3 '----------------------------------------------------------------------- 210 STR = DNU + 1# CKR = STR * RZR CKI = STR * RZI If N = 1 Then INU = INU - 1 If INU > 0 Then GoTo 220 If N > 1 Then GoTo 215 S1R = S2R S1I = S2I 215 ZDR = ZR ZDI = ZI If IFLAG = 1 Then GoTo 270 GoTo 240 220 INUB = 1 If IFLAG = 1 Then GoTo 261 225 P1R = CSRR(KFLAG) ASCLE = BRY(KFLAG) For I = INUB To 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 >= 3 Then GoTo 230 P2R = S2R * P1R P2I = S2I * P1R STR = Abs(P2R) STI = Abs(P2I) P2M = DMAX(STR, STI) If P2M <= ASCLE Then GoTo 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 Next I If N <> 1 Then GoTo 240 S1R = S2R S1I = S2I 240 STR = CSRR(KFLAG) YR(1) = S1R * STR YI(1) = S1I * STR If N = 1 Then GoTo 400 YR(2) = S2R * STR YI(2) = S2I * STR If N = 2 Then GoTo 400 KK = 2 250 KK = KK + 1 If KK > N Then GoTo 400 P1R = CSRR(KFLAG) ASCLE = BRY(KFLAG) For I = KK To 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 >= 3 Then GoTo 260 STR = Abs(P2R) STI = Abs(P2I) P2M = DMAX(STR, STI) If P2M <= ASCLE Then GoTo 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 Next I GoTo 400 '----------------------------------------------------------------------- ' IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW '----------------------------------------------------------------------- 261 HELIM = 0.5 * ELIM ELM = Exp(-ELIM) CELMR = ELM ASCLE = BRY(1) ZDR = ZR ZDI = ZI IC = -1 J = 2 For I = 1 To 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 AS0 = ZABS(S2R, S2I) ALAS = Log(AS0) P2R = -ZDR + ALAS If P2R < -ELIM Then GoTo 263 ZLOG S2R, S2I, STR, STI, IDUM P2R = -ZDR + STR P2I = -ZDI + STI P2M = Exp(P2R) / TOL P1R = P2M * Cos(P2I) P1I = P2M * Sin(P2I) ZUCHK P1R, P1I, NW, ASCLE, TOL If NW <> 0 Then GoTo 263 J = 3 - J YR(J) = P1R YI(J) = P1I If IC = I - 1 Then GoTo 264 IC = I GoTo 262 263 If ALAS < HELIM Then GoTo 262 ZDR = ZDR - ELIM S1R = S1R * CELMR S1I = S1I * CELMR S2R = S2R * CELMR S2I = S2I * CELMR 262 Next I If N <> 1 Then GoTo 270 S1R = S2R S1I = S2I GoTo 270 264 KFLAG = 1 INUB = I + 1 S2R = YR(J) S2I = YI(J) J = 3 - J S1R = YR(J) S1I = YI(J) If INUB <= INU Then GoTo 225 If N <> 1 Then GoTo 240 S1R = S2R S1I = S2I GoTo 240 270 YR(1) = S1R YI(1) = S1I If N = 1 Then GoTo 280 YR(2) = S2R YI(2) = S2I 280 ASCLE = BRY(1) ZKSCL ZDR, ZDI, FNU, N, YR(), YI(), NZ, RZR, RZI, ASCLE, TOL, ELIM INU = N - NZ If INU <= 0 Then GoTo 400 KK = NZ + 1 S1R = YR(KK) S1I = YI(KK) YR(KK) = S1R * CSRR(1) YI(KK) = S1I * CSRR(1) If INU = 1 Then GoTo 400 KK = NZ + 2 S2R = YR(KK) S2I = YI(KK) YR(KK) = S2R * CSRR(1) YI(KK) = S2I * CSRR(1) If INU = 2 Then GoTo 400 T2 = FNU + 1# * (KK - 1) CKR = T2 * RZR CKI = T2 * RZI KFLAG = 1 GoTo 250 '----------------------------------------------------------------------- ' SCALE BY DEXP(Z), IFLAG = 1 CASES '----------------------------------------------------------------------- 290 KODED = 2 IFLAG = 1 KFLAG = 2 GoTo 120 '----------------------------------------------------------------------- ' FNU=HALF ODD INTEGER CASE, DNU=-0.5 '----------------------------------------------------------------------- 300 S1R = COEFR S1I = COEFI S2R = COEFR S2I = COEFI GoTo 210 310 NZ = -2 400 End Sub 'ZBKNU Sub ZKSCL(ZRR, ZRI, FNU, N, YR(), YI(), NZ, RZR, RZI, ASCLE, TOL, ELIM) '***BEGIN PROLOGUE ZKSCL '***REFER TO ZBESK ' ' SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE ' ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN ' RETURN WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL. ' '***ROUTINES CALLED ZUCHK,ZABS,ZLOG '***END PROLOGUE ZKSCL ' COMPLEX CK,CS,Y,CZERO,RZ,S1,S2,Y,ZR,ZD,CELM 'Labels: 10,20,25,30,40,45,50 ' ACS, AS, CKI, CKR, CSI, CSR, FN, STR, S1I, S1R, S2I, S2R, ' ZEROI, ZEROR, ZDR, ZDI, CELMR, ELM, HELIM, ALAS: Double ' I, IC, IDUM, KK, NN, NW: Integer Dim STR As Double ZEROR = 0#: ZEROI = 0# NZ = 0 IC = 0 NN = IMIN(2, N) For I = 1 To NN S1R = YR(I) S1I = YI(I) YR(I) = S1R YI(I) = S1I AS0 = ZABS(S1R, S1I) ACS = -ZRR + Log(AS0) NZ = NZ + 1 YR(I) = ZEROR YI(I) = ZEROI If ACS < -ELIM Then GoTo 10 ZLOG S1R, S1I, CSR, CSI, IDUM CSR = CSR - ZRR CSI = CSI - ZRI STR = Exp(CSR) / TOL CSR = STR * Cos(CSI) CSI = STR * Sin(CSI) ZUCHK CSR, CSI, NW, ASCLE, TOL If NW <> 0 Then GoTo 10 YR(I) = CSR YI(I) = CSI IC = I NZ = NZ - 1 10 Next I If N = 1 Then GoTo 50 If IC > 1 Then GoTo 20 YR(1) = ZEROR YI(1) = ZEROI NZ = 2 20 If N = 2 Then GoTo 50 If NZ = 0 Then GoTo 50 FN = FNU + 1# CKR = FN * RZR CKI = FN * RZI S1R = YR(1) S1I = YI(1) S2R = YR(2) S2I = YI(2) HELIM = 0.5 * ELIM ELM = Exp(-ELIM) CELMR = ELM ZDR = ZRR ZDI = ZRI ' FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF ' S2 GETS LARGER THAN EXP(ELIM/2) For I = 3 To 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 AS0 = ZABS(S2R, S2I) ALAS = Log(AS0) ACS = -ZDR + ALAS NZ = NZ + 1 YR(I) = ZEROR YI(I) = ZEROI If ACS < -ELIM Then GoTo 25 ZLOG S2R, S2I, CSR, CSI, IDUM CSR = CSR - ZDR CSI = CSI - ZDI STR = Exp(CSR) / TOL CSR = STR * Cos(CSI) CSI = STR * Sin(CSI) ZUCHK CSR, CSI, NW, ASCLE, TOL If NW <> 0 Then GoTo 25 YR(I) = CSR YI(I) = CSI NZ = NZ - 1 If IC = KK - 1 Then GoTo 40 IC = KK GoTo 30 25 If ALAS < HELIM Then GoTo 30 ZDR = ZDR - ELIM S1R = S1R * CELMR S1I = S1I * CELMR S2R = S2R * CELMR S2I = S2I * CELMR 30 Next I NZ = N If IC = N Then NZ = N - 1 GoTo 45 40 NZ = KK - 2 45 For I = 1 To NZ YR(I) = ZEROR YI(I) = ZEROI Next I 50 End Sub 'ZKSCL 'end of module CBess00