Attribute VB_Name = "Module7" '************************************************************************** '* 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/16/2005). * '* (www.jpmoreau.fr) * '************************************************************************** 'MODULE CBESS3 DefDbl A-H, O-Z DefInt I-N 'LIST OF SUBROUTINES: ' ZBUNK(ZR, ZI, FNU, KODE, MR, N, YR(), YI(), NZ, TOL, ELIM, ALIM) ' ZACON(ZR, ZI, FNU, KODE, MR, N, YR(), YI(), NZ, RL, FNUL, TOL, ELIM, ALIM) ' ZUNK1(ZR, ZI, FNU, KODE, MR, N, YR(), YI(), NZ, TOL, ELIM, ALIM) ' ZUNK2(ZR, ZI, FNU, KODE, MR, N, YR(), YI(), NZ, TOL, ELIM, ALIM) Sub ZBUNK(ZR, ZI, FNU, KODE, MR, N, YR(), YI(), NZ, TOL, ELIM, ALIM) '***BEGIN PROLOGUE ZBUNK '***REFER TO ZBESK,ZBESH ' ' ZBUNK COMPUTES THE K BESSEL FUNCTION FOR FNU.GT.FNUL. ' ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR K(FNU,Z) ' IN ZUNK1 AND THE EXPANSION FOR H(2,FNU,Z) IN ZUNK2 ' '***ROUTINES CALLED ZUNK1,ZUNK2 '***END PROLOGUE ZBUNK ' COMPLEX Y,Z 'Labels: 10, 20 ' AX, AY: Double NZ = 0 AX = Abs(ZR) * 1.7321 AY = Abs(ZI) If AY > AX Then GoTo 10 '----------------------------------------------------------------------- ' ASYMPTOTIC EXPANSION FOR K(FNU,Z) FOR LARGE FNU APPLIED IN ' -PI/3 <= ARG(Z) <= PI/3 '----------------------------------------------------------------------- ZUNK1 ZR, ZI, FNU, KODE, MR, N, YR(), YI(), NZ, TOL, ELIM, ALIM GoTo 20 '----------------------------------------------------------------------- ' ASYMPTOTIC EXPANSION FOR H(2,FNU,Z*EXP(M*HPI)) FOR LARGE FNU ' APPLIED IN PI/3 < ABS(ARG(Z)) <= PI/2 WHERE M=+I OR -I ' AND HPI=PI/2 '----------------------------------------------------------------------- 10 ZUNK2 ZR, ZI, FNU, KODE, MR, N, YR(), YI(), NZ, TOL, ELIM, ALIM 20 End Sub 'ZBUNK Sub ZUNK1(ZR, ZI, FNU, KODE, MR, N, YR(), YI(), NZ, TOL, ELIM, ALIM) '***BEGIN PROLOGUE ZUNK1 '***REFER TO ZBESK ' ' ZUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE ' RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE ' UNIFORM ASYMPTOTIC EXPANSION. ' MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION. ' NZ=-1 MEANS AN OVERFLOW WILL OCCUR ' '***ROUTINES CALLED ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,ZABS '***END PROLOGUE ZUNK1 ' COMPLEX CFN,CK,CONE,CRSC,CS,CSCL,CSGN,CSPN,CSR,CSS,CWRK,CY,CZERO, ' C1,C2,PHI,PHID,RZ,SUM,SUMD,S1,S2,Y,Z,ZETA1,ZETA1D,ZETA2,ZETA2D,ZR 'Labels: 10,20,30,40,50,60,70,75,80,90,95,100,120,160,170,172,175,180, ' 200,210,220,230,250,255,260,270,275,280,290,300, 400 ' ANG, APHI, ASC, ASCLE, CKI, CKR, CONEI, CONER, CRSC, CSCL, CSGNI, ' CSPNI, CSPNR, CSR, C1I, C1R, C2I, C2M, C2R, FMR, FN, FNF, PHIDI, ' PHIDR, RAST, RAZR, RS1, RZI, RZR, SGN, STI, STR, SUMDI, SUMDR, ' S1I, S1R, S2I, S2R, ZEROI, ZEROR, ZET1DI, ZET1DR, ZET2DI, ZET2DR, ' ZRI, ZRR: Double ' I, IB, IFLAG, IFN, II, IL, INU, IUF, K, KDFLG, KFLAG, KK, NW, INITD, ' IC, IPARD, J, M: Integer Dim BRY(3), SUMR(10), SUMI(10), ZETA1R(10), ZETA1I(10), ZETA2R(10), ZETA2I(10) Dim CYR(10), CYI(10) Dim INIT(3) Dim CWRKR(16, 3), CWRKI(16, 3) Dim TMPR(16), TMPI(16) Dim CSSR(3), CSRR(3), PHIR(3), PHII(3) Dim STR As Double ZEROR = 0#: ZEROI = 0#: CONER = 1#: CONEI = 0# KDFLG = 1 NZ = 0 '----------------------------------------------------------------------- ' EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN ' THE UNDERFLOW LIMIT '----------------------------------------------------------------------- CSCL = 1# / TOL CRSC = TOL CSSR(1) = CSCL CSSR(2) = CONER CSSR(3) = CRSC CSRR(1) = CRSC CSRR(2) = CONER CSRR(3) = CSCL BRY(1) = 1000 * D1MACH(1) / TOL BRY(2) = 1# / BRY(1) BRY(3) = D1MACH(2) ZRR = ZR ZRI = ZI If ZR >= 0# Then GoTo 10 ZRR = -ZR ZRI = -ZI 10 J = 2 For I = 1 To N '----------------------------------------------------------------------- ' J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J '----------------------------------------------------------------------- J = 3 - J FN = FNU + 1# * (I - 1) INIT(J) = 0 For II = 1 To 16 TMPR(II) = CWRKR(II, J) TMPI(II) = CWRKI(II, J) Next II ZUNIK ZRR, ZRI, FN, 2, 0, TOL, INIT(J), PHIR(J), PHII(J), ZETA1R(J), ZETA1I(J), ZETA2R(J), ZETA2I(J), SUMR(J), SUMI(J), TMPR(), TMPI() For II = 1 To 16 CWRKR(II, J) = TMPR(II) CWRKI(II, J) = TMPI(II) Next II If KODE = 1 Then GoTo 20 STR = ZRR + ZETA2R(J) STI = ZRI + ZETA2I(J) RAST = FN / ZABS(STR, STI) STR = STR * RAST * RAST STI = -STI * RAST * RAST S1R = ZETA1R(J) - STR S1I = ZETA1I(J) - STI GoTo 30 20 S1R = ZETA1R(J) - ZETA2R(J) S1I = ZETA1I(J) - ZETA2I(J) 30 RS1 = S1R '----------------------------------------------------------------------- ' TEST FOR UNDERFLOW AND OVERFLOW '----------------------------------------------------------------------- If Abs(RS1) > ELIM Then GoTo 60 If KDFLG = 1 Then KFLAG = 2 If Abs(RS1) < ALIM Then GoTo 40 '----------------------------------------------------------------------- ' REFINE TEST AND SCALE '----------------------------------------------------------------------- APHI = ZABS(PHIR(J), PHII(J)) RS1 = RS1 + Log(APHI) If Abs(RS1) > ELIM Then GoTo 60 If KDFLG = 1 Then KFLAG = 1 If RS1 < 0# Then GoTo 40 If KDFLG = 1 Then KFLAG = 3 '----------------------------------------------------------------------- ' SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR ' EXPONENT EXTREMES '----------------------------------------------------------------------- 40 S2R = PHIR(J) * SUMR(J) - PHII(J) * SUMI(J) S2I = PHIR(J) * SUMI(J) + PHII(J) * SUMR(J) STR = Exp(S1R) * CSSR(KFLAG) S1R = STR * Cos(S1I) S1I = STR * Sin(S1I) STR = S2R * S1R - S2I * S1I S2I = S1R * S2I + S2R * S1I S2R = STR If KFLAG <> 1 Then GoTo 50 ZUCHK S2R, S2I, NW, BRY(1), TOL If NW <> 0 Then GoTo 60 50 CYR(KDFLG) = S2R CYI(KDFLG) = S2I YR(I) = S2R * CSRR(KFLAG) YI(I) = S2I * CSRR(KFLAG) If KDFLG = 2 Then GoTo 75 KDFLG = 2 GoTo 70 60 If RS1 > 0# Then GoTo 300 '----------------------------------------------------------------------- ' FOR ZR < 0, THE I FUNCTION TO BE ADDED WILL OVERFLOW '----------------------------------------------------------------------- If ZR < 0# Then GoTo 300 KDFLG = 1 YR(I) = ZEROR YI(I) = ZEROI NZ = NZ + 1 If I = 1 Then GoTo 70 If (YR(I - 1) = ZEROR) And (YI(I - 1) = ZEROI) Then GoTo 70 YR(I - 1) = ZEROR YI(I - 1) = ZEROI NZ = NZ + 1 70 Next I I = N 75 RAZR = 1# / ZABS(ZRR, ZRI) STR = ZRR * RAZR STI = -ZRI * RAZR RZR = (STR + STR) * RAZR RZI = (STI + STI) * RAZR CKR = FN * RZR CKI = FN * RZI ib = I + 1 If N < ib Then GoTo 160 '----------------------------------------------------------------------- ' TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO ' ON UNDERFLOW. '----------------------------------------------------------------------- FN = FNU + 1# * (N - 1) IPARD = 1 If MR <> 0 Then IPARD = 0 INITD = 0 For II = 1 To 16 TMPR(II) = CWRKR(II, 3) TMPI(II) = CWRKI(II, 3) Next II ZUNIK ZRR, ZRI, FN, 2, IPARD, TOL, INITD, PHIDR, PHIDI, ZET1DR, ZET1DI, ZET2DR, ZET2DI, SUMDR, SUMDI, TMPR(), TMPI() For II = 1 To 16 CWRKR(II, 3) = TMPR(II) CWRKI(II, 3) = TMPI(II) Next II If KODE = 1 Then GoTo 80 STR = ZRR + ZET2DR STI = ZRI + ZET2DI RAST = FN / ZABS(STR, STI) STR = STR * RAST * RAST STI = -STI * RAST * RAST S1R = ZET1DR - STR S1I = ZET1DI - STI GoTo 90 80 S1R = ZET1DR - ZET2DR S1I = ZET1DI - ZET2DI 90 RS1 = S1R If Abs(RS1) > ELIM Then GoTo 95 If Abs(RS1) < ALIM Then GoTo 100 '----------------------------------------------------------------------- ' REFINE ESTIMATE AND TEST '----------------------------------------------------------------------- APHI = ZABS(PHIDR, PHIDI) RS1 = RS1 + Log(APHI) If Abs(RS1) < ELIM Then GoTo 100 95 If Abs(RS1) > 0# Then GoTo 300 '----------------------------------------------------------------------- ' FOR ZR < 0, THE I FUNCTION TO BE ADDED WILL OVERFLOW '----------------------------------------------------------------------- If ZR < 0# Then GoTo 300 NZ = N For I = 1 To N YR(I) = ZEROR YI(I) = ZEROI Next I GoTo 400 '----------------------------------------------------------------------- ' FORWARD RECUR FOR REMAINDER OF THE SEQUENCE '----------------------------------------------------------------------- 100 S1R = CYR(1) S1I = CYI(1) S2R = CYR(2) S2I = CYI(2) C1R = CSRR(KFLAG) ASCLE = BRY(KFLAG) For I = ib To N C2R = S2R C2I = S2I S2R = CKR * C2R - CKI * C2I + S1R S2I = CKR * C2I + CKI * C2R + S1I S1R = C2R S1I = C2I CKR = CKR + RZR CKI = CKI + RZI C2R = S2R * C1R C2I = S2I * C1R YR(I) = C2R YI(I) = C2I If KFLAG >= 3 Then GoTo 120 STR = Abs(C2R) STI = Abs(C2I) C2M = DMAX(STR, STI) If C2M <= ASCLE Then GoTo 120 KFLAG = KFLAG + 1 ASCLE = BRY(KFLAG) S1R = S1R * C1R S1I = S1I * C1R S2R = C2R S2I = C2I S1R = S1R * CSSR(KFLAG) S1I = S1I * CSSR(KFLAG) S2R = S2R * CSSR(KFLAG) S2I = S2I * CSSR(KFLAG) C1R = CSRR(KFLAG) 120 Next I 160 If MR = 0 Then GoTo 400 '----------------------------------------------------------------------- ' ANALYTIC CONTINUATION FOR RE(Z) < 0 '----------------------------------------------------------------------- NZ = 0 FMR = 1# * MR SGN0 = -Sign(PI, FMR) '----------------------------------------------------------------------- ' CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP. '----------------------------------------------------------------------- CSGNI = SGN0 INU = Int(FNU) FNF = FNU - 1# * INU IFN = INU + N - 1 ANG = FNF * SGN0 CSPNR = Cos(ANG) CSPNI = Sin(ANG) If (IFN Mod 2) = 0 Then GoTo 170 CSPNR = -CSPNR CSPNI = -CSPNI 170 Asc0 = BRY(1) IUF = 0 KK = N KDFLG = 1 ib = ib - 1 IC = ib - 1 For K = 1 To N FN = FNU + 1# * (KK - 1) '----------------------------------------------------------------------- ' LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K ' FUNCTION ABOVE '----------------------------------------------------------------------- M = 3 If N > 2 Then GoTo 175 172 INITD = INIT(J) PHIDR = PHIR(J) PHIDI = PHII(J) ZET1DR = ZETA1R(J) ZET1DI = ZETA1I(J) ZET2DR = ZETA2R(J) ZET2DI = ZETA2I(J) SUMDR = SUMR(J) SUMDI = SUMI(J) M = J J = 3 - J GoTo 180 175 If (KK = N) And (ib < N) Then GoTo 180 If (KK = ib) Or (KK = IC) Then GoTo 172 INITD = 0 180 For II = 1 To 16 TMPR(II) = CWRKR(II, M) TMPI(II) = CWRKI(II, M) Next II ZUNIK ZRR, ZRI, FN, 1, 0, TOL, INITD, PHIDR, PHIDI, ZET1DR, ZET1DI, ZET2DR, ZET2DI, SUMDR, SUMDI, TMPR(), TMPI() For II = 1 To 16 CWRKR(II, M) = TMPR(II) CWRKI(II, M) = TMPI(II) Next II If KODE = 1 Then GoTo 200 STR = ZRR + ZET2DR STI = ZRI + ZET2DI RAST = FN / ZABS(STR, STI) STR = STR * RAST * RAST STI = -STI * RAST * RAST S1R = -ZET1DR + STR S1I = -ZET1DI + STI GoTo 210 200 S1R = -ZET1DR + ZET2DR S1I = -ZET1DI + ZET2DI '----------------------------------------------------------------------- ' TEST FOR UNDERFLOW AND OVERFLOW '----------------------------------------------------------------------- 210 RS1 = S1R If Abs(RS1) > ELIM Then GoTo 260 If KDFLG = 1 Then IFLAG = 2 If Abs(RS1) < ALIM Then GoTo 220 '----------------------------------------------------------------------- ' REFINE TEST AND SCALE '----------------------------------------------------------------------- APHI = ZABS(PHIDR, PHIDI) RS1 = RS1 + Log(APHI) If Abs(RS1) > ELIM Then GoTo 260 If KDFLG = 1 Then IFLAG = 1 If RS1 < 0# Then GoTo 220 If KDFLG = 1 Then IFLAG = 3 220 STR = PHIDR * SUMDR - PHIDI * SUMDI STI = PHIDR * SUMDI + PHIDI * SUMDR S2R = -CSGNI * STI S2I = CSGNI * STR STR = Exp(S1R) * CSSR(IFLAG) S1R = STR * Cos(S1I) S1I = STR * Sin(S1I) STR = S2R * S1R - S2I * S1I S2I = S2R * S1I + S2I * S1R S2R = STR If IFLAG <> 1 Then GoTo 230 ZUCHK S2R, S2I, NW, BRY(1), TOL If NW = 0 Then GoTo 230 S2R = ZEROR S2I = ZEROI 230 CYR(KDFLG) = S2R CYI(KDFLG) = S2I C2R = S2R C2I = S2I S2R = S2R * CSRR(IFLAG) S2I = S2I * CSRR(IFLAG) '----------------------------------------------------------------------- ' ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N '----------------------------------------------------------------------- S1R = YR(KK) S1I = YI(KK) If KODE = 1 Then GoTo 250 ZS1S2 ZRR, ZRI, S1R, S1I, S2R, S2I, NW, Asc0, ALIM, IUF NZ = NZ + NW 250 YR(KK) = S1R * CSPNR - S1I * CSPNI + S2R YI(KK) = CSPNR * S1I + CSPNI * S1R + S2I KK = KK - 1 CSPNR = -CSPNR CSPNI = -CSPNI If (C2R <> 0#) Or (C2I <> 0#) Then GoTo 255 KDFLG = 1 GoTo 270 255 If KDFLG = 2 Then GoTo 275 KDFLG = 2 GoTo 270 260 If RS1 > 0# Then GoTo 300 S2R = ZEROR S2I = ZEROI GoTo 230 270 Next K K = N 275: IL = N - K If IL = 0 Then GoTo 400 '----------------------------------------------------------------------- ' RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE ' K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP ' INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES. '----------------------------------------------------------------------- S1R = CYR(1) S1I = CYI(1) S2R = CYR(2) S2I = CYI(2) CSR = CSRR(IFLAG) ASCLE = BRY(IFLAG) FN = 1# * (INU + IL) For I = 1 To IL C2R = S2R C2I = S2I S2R = S1R + (FN + FNF) * (RZR * C2R - RZI * C2I) S2I = S1I + (FN + FNF) * (RZR * C2I + RZI * C2R) S1R = C2R S1I = C2I FN = FN - 1# C2R = S2R * CSR C2I = S2I * CSR CKR = C2R CKI = C2I C1R = YR(KK) C1I = YI(KK) If KODE = 1 Then GoTo 280 ZS1S2 ZRR, ZRI, C1R, C1I, C2R, C2I, NW, Asc0, ALIM, IUF NZ = NZ + NW 280 YR(KK) = C1R * CSPNR - C1I * CSPNI + C2R YI(KK) = C1R * CSPNI + C1I * CSPNR + C2I KK = KK - 1 CSPNR = -CSPNR CSPNI = -CSPNI If IFLAG >= 3 Then GoTo 290 C2R = Abs(CKR) C2I = Abs(CKI) C2M = DMAX(C2R, C2I) If C2M <= ASCLE Then GoTo 290 IFLAG = IFLAG + 1 ASCLE = BRY(IFLAG) S1R = S1R * CSR S1I = S1I * CSR S2R = CKR S2I = CKI S1R = S1R * CSSR(IFLAG) S1I = S1I * CSSR(IFLAG) S2R = S2R * CSSR(IFLAG) S2I = S2I * CSSR(IFLAG) CSR = CSRR(IFLAG) 290 Next I GoTo 400 300 NZ = -1 400 End Sub 'ZUNK1 Sub ZUNK2(ZR, ZI, FNU, KODE, MR, N, YR(), YI(), NZ, TOL, ELIM, ALIM) '***BEGIN PROLOGUE ZUNK2 '***REFER TO ZBESK ' ' ZUNK2 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE ' RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE ' UNIFORM ASYMPTOTIC EXPANSIONS FOR H(KIND,FNU,ZN) AND J(FNU,ZN) ' WHERE ZN IS IN THE RIGHT HALF PLANE, KIND=(3-MR)/2, MR=+1 OR ' -1. HERE ZN=ZR*I OR -ZR*I WHERE ZR=Z IF Z IS IN THE RIGHT ' HALF PLANE OR ZR=-Z IF Z IS IN THE LEFT HALF PLANE. MR INDIC- ' ATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION. ' NZ=-1 MEANS AN OVERFLOW WILL OCCUR ' '***ROUTINES CALLED ZAIRY,ZKSCL,ZS1S2,ZUCHK,ZUNHJ,D1MACH,ZABS '***END PROLOGUE ZUNK2 ' COMPLEX AI,ARG,ARGD,ASUM,ASUMD,BSUM,BSUMD,CFN,CI,CIP,CK,CONE,CRSC, ' CR1,CR2,CS,CSCL,CSGN,CSPN,CSR,CSS,CY,CZERO,C1,C2,DAI,PHI,PHID,RZ, ' S1,S2,Y,Z,ZB,ZETA1,ZETA1D,ZETA2,ZETA2D,ZN,ZR 'Labels: 10,20,30,40,50,60,70,80,85,90,100,105,120,130,172,175,180,190,210, ' 220,230,240,250,255,270,280,290,295,300,310,320, 400 ' AARG, AIC, AII, AIR, ANG, APHI, ARGDI, ARGDR, ASC, ASCLE, ASUMDI, ' ASUMDR, BSUMDI, BSUMDR,CAR, CKI, CKR, CONEI, CONER, CRSC, CR1I, CR1R, ' CR2I, CR2R, CSCL, CSGNI, CSI, CSPNI, CSPNR, CSR, C1I, C1R, C2I, C2M, ' C2R, DAII, DAIR, FMR, FN, FNF, HPI, PHIDI, PHIDR, PTI, PTR, RAST, ' RAZR, RS1, RZI, RZR, SAR, SGN, STI, STR, S1I, S1R, S2I, S2R, YY, ' ZBI, ZBR, ZEROI, ZEROR, ZET1DI, ZET1DR, ZET2DI, ZET2DR, ZNI, ZNR, ' ZRI, ZRR: Double ' I, IB, IFLAG, IFN, IL, IN0, INU, IUF, K, KDFLG, KFLAG, KK, NAI, ' NDAI, NW, IDUM, J, IPARD, IC: Integer Dim BRY(3), ASUMR(10), ASUMI(10), BSUMR(10), BSUMI(10), PHIR(10), PHII(10), ARGR(10), ARGI(10) Dim ZETA1R(10), ZETA1I(10), ZETA2R(10), ZETA2I(10), CIPR(4), CIPI(4), CSSR(3), CSRR(3) Dim CYR(10), CYI(10) Dim STR As Double ZEROR = 0#: ZEROI = 0#: CONER = 1#: CONEI = 0# CR1R = 1#: CR1I = 1.73205080756888 CR2R = -0.5: CR2I = -0.866025403784439 HPI = 1.5707963267949 AIC = 1.26551212348465 CIPR(1) = 1#: CIPI(1) = 0#: CIPR(2) = 0#: CIPI(2) = -1# CIPR(3) = -1#: CIPI(3) = 0#: CIPR(4) = 0#: CIPI(4) = 1# KDFLG = 1 NZ = 0 '----------------------------------------------------------------------- ' EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN ' THE UNDERFLOW LIMIT '----------------------------------------------------------------------- CSCL = 1# / TOL CRSC = TOL CSSR(1) = CSCL CSSR(2) = CONER CSSR(3) = CRSC CSRR(1) = CRSC CSRR(2) = CONER CSRR(3) = CSCL BRY(1) = 1000# * D1MACH(1) / TOL BRY(2) = 1# / BRY(1) BRY(3) = D1MACH(2) ZRR = ZR ZRI = ZI If ZR >= 0# Then GoTo 10 ZRR = -ZR ZRI = -ZI 10 YY = ZRI ZNR = ZRI ZNI = -ZRR ZBR = ZRR ZBI = ZRI INU = Int(FNU) FNF = FNU - 1# * INU ANG = -HPI * FNF CAR = Cos(ANG) SAR = Sin(ANG) C2R = HPI * SAR C2I = -HPI * CAR KK = (INU Mod 4) + 1 STR = C2R * CIPR(KK) - C2I * CIPI(KK) STI = C2R * CIPI(KK) + C2I * CIPR(KK) CSR = CR1R * STR - CR1I * STI CSI = CR1R * STI + CR1I * STR If YY > 0# Then GoTo 20 ZNR = -ZNR ZBI = -ZBI '----------------------------------------------------------------------- ' K(FNU,Z) IS COMPUTED FROM H(2,FNU,-I*Z) WHERE Z IS IN THE FIRST ' QUADRANT. FOURTH QUADRANT VALUES (YY <= 0) ARE COMPUTED BY ' CONJUGATION SINCE THE K FUNCTION IS REAL ON THE POSITIVE REAL AXIS '----------------------------------------------------------------------- 20 J = 2 For I = 1 To N '----------------------------------------------------------------------- ' J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J '----------------------------------------------------------------------- J = 3 - J FN = FNU + 1# * (I - 1) ZUNHJ ZNR, ZNI, FN, 0, TOL, PHIR(J), PHII(J), ARGR(J), ARGI(J), ZETA1R(J), ZETA1I(J), ZETA2R(J), ZETA2I(J), ASUMR(J), ASUMI(J), BSUMR(J), BSUMI(J) If KODE = 1 Then GoTo 30 STR = ZBR + ZETA2R(J) STI = ZBI + ZETA2I(J) RAST = FN / ZABS(STR, STI) STR = STR * RAST * RAST STI = -STI * RAST * RAST S1R = ZETA1R(J) - STR S1I = ZETA1I(J) - STI GoTo 40 30 S1R = ZETA1R(J) - ZETA2R(J) S1I = ZETA1I(J) - ZETA2I(J) '----------------------------------------------------------------------- ' TEST FOR UNDERFLOW AND OVERFLOW '----------------------------------------------------------------------- 40 RS1 = S1R If Abs(RS1) > ELIM Then GoTo 70 If KDFLG = 1 Then KFLAG = 2 If Abs(RS1) < ALIM Then GoTo 50 '----------------------------------------------------------------------- ' REFINE TEST AND SCALE '----------------------------------------------------------------------- APHI = ZABS(PHIR(J), PHII(J)) AARG = ZABS(ARGR(J), ARGI(J)) RS1 = RS1 + Log(APHI) - 0.25 * Log(AARG) - AIC If Abs(RS1) > ELIM Then GoTo 70 If KDFLG = 1 Then KFLAG = 1 If RS1 < 0# Then GoTo 50 If KDFLG = 1 Then KFLAG = 3 '----------------------------------------------------------------------- ' SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR ' EXPONENT EXTREMES '----------------------------------------------------------------------- 50 C2R = ARGR(J) * CR2R - ARGI(J) * CR2I C2I = ARGR(J) * CR2I + ARGI(J) * CR2R ZAIRY C2R, C2I, 0, 2, AIR, AII, NAI, IDUM ZAIRY C2R, C2I, 1, 2, DAIR, DAII, NDAI, IDUM STR = DAIR * BSUMR(J) - DAII * BSUMI(J) STI = DAIR * BSUMI(J) + DAII * BSUMR(J) PTR = STR * CR2R - STI * CR2I PTI = STR * CR2I + STI * CR2R STR = PTR + (AIR * ASUMR(J) - AII * ASUMI(J)) STI = PTI + (AIR * ASUMI(J) + AII * ASUMR(J)) PTR = STR * PHIR(J) - STI * PHII(J) PTI = STR * PHII(J) + STI * PHIR(J) S2R = PTR * CSR - PTI * CSI S2I = PTR * CSI + PTI * CSR STR = Exp(S1R) * CSSR(KFLAG) S1R = STR * Cos(S1I) S1I = STR * Sin(S1I) STR = S2R * S1R - S2I * S1I S2I = S1R * S2I + S2R * S1I S2R = STR If KFLAG <> 1 Then GoTo 60 ZUCHK S2R, S2I, NW, BRY(1), TOL If NW <> 0 Then GoTo 70 60 If YY <= 0# Then S2I = -S2I CYR(KDFLG) = S2R CYI(KDFLG) = S2I YR(I) = S2R * CSRR(KFLAG) YI(I) = S2I * CSRR(KFLAG) STR = CSI CSI = -CSR CSR = STR If KDFLG = 2 Then GoTo 85 KDFLG = 2 GoTo 80 70 If RS1 > 0# Then GoTo 320 '----------------------------------------------------------------------- ' FOR ZR < 0, THE I FUNCTION TO BE ADDED WILL OVERFLOW '----------------------------------------------------------------------- If ZR < 0# Then GoTo 320 KDFLG = 1 YR(I) = ZEROR YI(I) = ZEROI NZ = NZ + 1 STR = CSI CSI = -CSR CSR = STR If I = 1 Then GoTo 80 If (YR(I - 1) = ZEROR) And (YI(I - 1) = ZEROI) Then GoTo 80 YR(I - 1) = ZEROR YI(I - 1) = ZEROI NZ = NZ + 1 80 Next I I = N 85 RAZR = 1# / ZABS(ZRR, ZRI) STR = ZRR * RAZR STI = -ZRI * RAZR RZR = (STR + STR) * RAZR RZI = (STI + STI) * RAZR CKR = FN * RZR CKI = FN * RZI ib = I + 1 If N < ib Then GoTo 180 '----------------------------------------------------------------------- ' TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO ' ON UNDERFLOW. '----------------------------------------------------------------------- FN = FNU + 1# * (N - 1) IPARD = 1 If MR <> 0 Then IPARD = 0 ZUNHJ ZNR, ZNI, FN, IPARD, TOL, PHIDR, PHIDI, ARGDR, ARGDI, ZET1DR, ZET1DI, ZET2DR, ZET2DI, ASUMDR, ASUMDI, BSUMDR, BSUMDI If KODE = 1 Then GoTo 90 STR = ZBR + ZET2DR STI = ZBI + ZET2DI RAST = FN / ZABS(STR, STI) STR = STR * RAST * RAST STI = -STI * RAST * RAST S1R = ZET1DR - STR S1I = ZET1DI - STI GoTo 100 90 S1R = ZET1DR - ZET2DR S1I = ZET1DI - ZET2DI 100 RS1 = S1R If Abs(RS1) > ELIM Then GoTo 105 If Abs(RS1) < ALIM Then GoTo 120 '----------------------------------------------------------------------- ' REFINE ESTIMATE AND TEST '----------------------------------------------------------------------- APHI = ZABS(PHIDR, PHIDI) RS1 = RS1 + Log(APHI) If Abs(RS1) < ELIM Then GoTo 120 105 If RS1 > 0# Then GoTo 320 '----------------------------------------------------------------------- ' FOR ZR < 0, THE I FUNCTION TO BE ADDED WILL OVERFLOW '----------------------------------------------------------------------- If ZR < 0# Then GoTo 320 NZ = N For I = 1 To N YR(I) = ZEROR YI(I) = ZEROI Next I GoTo 400 120 S1R = CYR(1) S1I = CYI(1) S2R = CYR(2) S2I = CYI(2) C1R = CSRR(KFLAG) ASCLE = BRY(KFLAG) For I = ib To N C2R = S2R C2I = S2I S2R = CKR * C2R - CKI * C2I + S1R S2I = CKR * C2I + CKI * C2R + S1I S1R = C2R S1I = C2I CKR = CKR + RZR CKI = CKI + RZI C2R = S2R * C1R C2I = S2I * C1R YR(I) = C2R YI(I) = C2I If KFLAG >= 3 Then GoTo 130 STR = Abs(C2R) STI = Abs(C2I) C2M = DMAX(STR, STI) If C2M <= ASCLE Then GoTo 130 KFLAG = KFLAG + 1 ASCLE = BRY(KFLAG) S1R = S1R * C1R S1I = S1I * C1R S2R = C2R S2I = C2I S1R = S1R * CSSR(KFLAG) S1I = S1I * CSSR(KFLAG) S2R = S2R * CSSR(KFLAG) S2I = S2I * CSSR(KFLAG) C1R = CSRR(KFLAG) 130 Next I 180 If MR = 0 Then GoTo 400 '----------------------------------------------------------------------- ' ANALYTIC CONTINUATION FOR RE(Z) < 0 '----------------------------------------------------------------------- NZ = 0 FMR = 1# * MR SGN0 = -Sign(PI, FMR) '----------------------------------------------------------------------- ' CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP. '----------------------------------------------------------------------- CSGNI = SGN0 If YY <= 0# Then CSGNI = -CSGNI IFN = INU + N - 1 ANG = FNF * SGN0 CSPNR = Cos(ANG) CSPNI = Sin(ANG) If (IFN Mod 2) = 0 Then GoTo 190 CSPNR = -CSPNR CSPNI = -CSPNI '----------------------------------------------------------------------- ' CS=COEFF OF THE J FUNCTION TO GET THE I FUNCTION. I(FNU,Z) IS ' COMPUTED FROM EXP(I*FNU*HPI)*J(FNU,-I*Z) WHERE Z IS IN THE FIRST ' QUADRANT. FOURTH QUADRANT VALUES (YY <= 0) ARE COMPUTED BY ' CONJUGATION SINCE THE I FUNCTION IS REAL ON THE POSITIVE REAL AXIS '----------------------------------------------------------------------- 190: CSR = SAR * CSGNI CSI = CAR * CSGNI IN0 = (IFN Mod 4) + 1 C2R = CIPR(IN0) C2I = CIPI(IN0) STR = CSR * C2R + CSI * C2I CSI = -CSR * C2I + CSI * C2R CSR = STR Asc0 = BRY(1) IUF = 0 KK = N KDFLG = 1 ib = ib - 1 IC = ib - 1 For K = 1 To N FN = FNU + 1# * (KK - 1) '----------------------------------------------------------------------- ' LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K ' FUNCTION ABOVE '----------------------------------------------------------------------- If N > 2 Then GoTo 175 172 PHIDR = PHIR(J) PHIDI = PHII(J) ARGDR = ARGR(J) ARGDI = ARGI(J) ZET1DR = ZETA1R(J) ZET1DI = ZETA1I(J) ZET2DR = ZETA2R(J) ZET2DI = ZETA2I(J) ASUMDR = ASUMR(J) ASUMDI = ASUMI(J) BSUMDR = BSUMR(J) BSUMDI = BSUMI(J) J = 3 - J GoTo 210 175 If (KK = N) And (ib < N) Then GoTo 210 If (KK = ib) Or (KK = IC) Then GoTo 172 ZUNHJ ZNR, ZNI, FN, 0, TOL, PHIDR, PHIDI, ARGDR, ARGDI, ZET1DR, ZET1DI, ZET2DR, ZET2DI, ASUMDR, ASUMDI, BSUMDR, BSUMDI 210 If KODE = 1 Then GoTo 220 STR = ZBR + ZET2DR STI = ZBI + ZET2DI RAST = FN / ZABS(STR, STI) STR = STR * RAST * RAST STI = -STI * RAST * RAST S1R = -ZET1DR + STR S1I = -ZET1DI + STI GoTo 230 220 S1R = -ZET1DR + ZET2DR S1I = -ZET1DI + ZET2DI '----------------------------------------------------------------------- ' TEST FOR UNDERFLOW AND OVERFLOW '----------------------------------------------------------------------- 230 RS1 = S1R If Abs(RS1) > ELIM Then GoTo 280 If KDFLG = 1 Then IFLAG = 2 If Abs(RS1) < ALIM Then GoTo 240 '----------------------------------------------------------------------- ' REFINE TEST AND SCALE '----------------------------------------------------------------------- APHI = ZABS(PHIDR, PHIDI) AARG = ZABS(ARGDR, ARGDI) RS1 = RS1 + Log(APHI) - 0.25 * Log(AARG) - AIC If Abs(RS1) > ELIM Then GoTo 280 If KDFLG = 1 Then IFLAG = 1 If RS1 < 0# Then GoTo 240 If KDFLG = 1 Then IFLAG = 3 240 ZAIRY ARGDR, ARGDI, 0, 2, AIR, AII, NAI, IDUM ZAIRY ARGDR, ARGDI, 1, 2, DAIR, DAII, NDAI, IDUM STR = DAIR * BSUMDR - DAII * BSUMDI STI = DAIR * BSUMDI + DAII * BSUMDR STR = STR + (AIR * ASUMDR - AII * ASUMDI) STI = STI + (AIR * ASUMDI + AII * ASUMDR) PTR = STR * PHIDR - STI * PHIDI PTI = STR * PHIDI + STI * PHIDR S2R = PTR * CSR - PTI * CSI S2I = PTR * CSI + PTI * CSR STR = Exp(S1R) * CSSR(IFLAG) S1R = STR * Cos(S1I) S1I = STR * Sin(S1I) STR = S2R * S1R - S2I * S1I S2I = S2R * S1I + S2I * S1R S2R = STR If IFLAG <> 1 Then GoTo 250 ZUCHK S2R, S2I, NW, BRY(1), TOL If NW = 0 Then GoTo 250 S2R = ZEROR S2I = ZEROI 250 If YY <= 0# Then S2I = -S2I CYR(KDFLG) = S2R CYI(KDFLG) = S2I C2R = S2R C2I = S2I S2R = S2R * CSRR(IFLAG) S2I = S2I * CSRR(IFLAG) '----------------------------------------------------------------------- ' ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N '----------------------------------------------------------------------- S1R = YR(KK) S1I = YI(KK) If KODE = 1 Then GoTo 270 ZS1S2 ZRR, ZRI, S1R, S1I, S2R, S2I, NW, Asc0, ALIM, IUF NZ = NZ + NW 270 YR(KK) = S1R * CSPNR - S1I * CSPNI + S2R YI(KK) = S1R * CSPNI + S1I * CSPNR + S2I KK = KK - 1 CSPNR = -CSPNR CSPNI = -CSPNI STR = CSI CSI = -CSR CSR = STR If (C2R <> 0#) Or (C2I <> 0#) Then GoTo 255 KDFLG = 1 GoTo 290 255: If KDFLG = 2 Then GoTo 295 KDFLG = 2 GoTo 290 280: If RS1 > 0# Then GoTo 320 S2R = ZEROR S2I = ZEROI GoTo 250 290 Next K K = N 295 IL = N - K If IL = 0 Then GoTo 400 '----------------------------------------------------------------------- ' RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE ' K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP ' INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES. '----------------------------------------------------------------------- S1R = CYR(1) S1I = CYI(1) S2R = CYR(2) S2I = CYI(2) CSR = CSRR(IFLAG) ASCLE = BRY(IFLAG) FN = 1# * (INU + IL) For I = 1 To IL C2R = S2R C2I = S2I S2R = S1R + (FN + FNF) * (RZR * C2R - RZI * C2I) S2I = S1I + (FN + FNF) * (RZR * C2I + RZI * C2R) S1R = C2R S1I = C2I FN = FN - 1# C2R = S2R * CSR C2I = S2I * CSR CKR = C2R CKI = C2I C1R = YR(KK) C1I = YI(KK) If KODE = 1 Then GoTo 300 ZS1S2 ZRR, ZRI, C1R, C1I, C2R, C2I, NW, Asc0, ALIM, IUF NZ = NZ + NW 300: YR(KK) = C1R * CSPNR - C1I * CSPNI + C2R YI(KK) = C1R * CSPNI + C1I * CSPNR + C2I KK = KK - 1 CSPNR = -CSPNR CSPNI = -CSPNI If IFLAG >= 3 Then GoTo 310 C2R = Abs(CKR) C2I = Abs(CKI) C2M = DMAX(C2R, C2I) If C2M <= ASCLE Then GoTo 310 IFLAG = IFLAG + 1 ASCLE = BRY(IFLAG) S1R = S1R * CSR S1I = S1I * CSR S2R = CKR S2I = CKI S1R = S1R * CSSR(IFLAG) S1I = S1I * CSSR(IFLAG) S2R = S2R * CSSR(IFLAG) S2I = S2I * CSSR(IFLAG) CSR = CSRR(IFLAG) 310 Next I GoTo 400 320 NZ = -1 400 End Sub 'ZUNK2 Sub ZACON(ZR, ZI, FNU, KODE, MR, N, YR(), YI(), NZ, RL, FNUL, TOL, ELIM, ALIM) '***BEGIN PROLOGUE ZACON '***REFER TO ZBESK,ZBESH ' ' ZACON 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.0,1.0) ' ' TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT ' HALF Z PLANE ' '***ROUTINES CALLED ZBINU,ZBKNU,ZS1S2,D1MACH,ZABS,ZMLT '***END PROLOGUE ZACON ' COMPLEX CK,CONE,CSCL,CSCR,CSGN,CSPN,CY,CZERO,C1,C2,RZ,SC1,SC2,ST, ' S1,S2,Y,Z,ZN 'Labels: 10,20,30,40,50,60,70,80,90, 100 ' ARG, ASCLE, AS2, AZN, BSCLE, CKI, CKR, CONEI, CONER, CPN, CSCL, CSCR, ' CSGNI, CSGNR, CSPNI, CSPNR, CSR, C1I, C1M, C1R, C2I, C2R, FMR, FN, ' PTI, PTR, RAZN, RZI, RZR, SC1I, SC1R, SC2I, SC2R, SGN, SPN, STI, STR, ' S1I, S1R, S2I, S2R, YY, ZEROI, ZEROR, ZNI, ZNR: Double ' I, INU, IUF, KFLAG, NN, NW: Integer Dim CYR(10), CYI(10), CSSR(4), CSRR(4), BRY(3) Dim STR As Double PI=4#*Atn(1#) ZEROR = 0#: ZEROI = 0#: CONER = 1#: CONEI = 0# NZ = 0 ZNR = -ZR ZNI = -ZI NN = N ZBINU ZNR, ZNI, FNU, KODE, NN, YR(), YI(), NW, RL, FNUL, TOL, ELIM, ALIM If NW < 0 Then GoTo 90 '----------------------------------------------------------------------- ' ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION '----------------------------------------------------------------------- NN = IMIN(2, N) ZBKNU ZNR, ZNI, FNU, KODE, NN, CYR(), CYI(), NW, TOL, ELIM, ALIM If NW <> 0 Then GoTo 90 S1R = CYR(1) S1I = CYI(1) FMR = 1# * MR SGN0 = -Sign(PI, FMR) CSGNR = ZEROR CSGNI = SGN0 If KODE = 1 Then GoTo 10 YY = -ZNI CPN = Cos(YY) SPN = Sin(YY) ZMLT CSGNR, CSGNI, CPN, SPN, CSGNR, CSGNI '----------------------------------------------------------------------- ' CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE ' WHEN FNU IS LARGE '----------------------------------------------------------------------- 10 INU = Int(FNU) ARG = (FNU - 1# * INU) * SGN0 CPN = Cos(ARG) SPN = Sin(ARG) CSPNR = CPN CSPNI = SPN If (INU Mod 2) = 0 Then GoTo 20 CSPNR = -CSPNR CSPNI = -CSPNI 20 IUF = 0 C1R = S1R C1I = S1I C2R = YR(1) C2I = YI(1) ASCLE = 1000# * D1MACH(1) / TOL If KODE = 1 Then GoTo 30 ZS1S2 ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF NZ = NZ + NW SC1R = C1R SC1I = C1I 30 ZMLT CSPNR, CSPNI, C1R, C1I, STR, STI ZMLT CSGNR, CSGNI, C2R, C2I, PTR, PTI YR(1) = STR + PTR YI(1) = STI + PTI If N = 1 Then GoTo 100 CSPNR = -CSPNR CSPNI = -CSPNI S2R = CYR(2) S2I = CYI(2) C1R = S2R C1I = S2I C2R = YR(2) C2I = YI(2) If KODE = 1 Then GoTo 40 ZS1S2 ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF NZ = NZ + NW SC2R = C1R SC2I = C1I 40 ZMLT CSPNR, CSPNI, C1R, C1I, STR, STI ZMLT CSGNR, CSGNI, C2R, C2I, PTR, PTI YR(2) = STR + PTR YI(2) = STI + PTI If N = 2 Then GoTo 100 CSPNR = -CSPNR CSPNI = -CSPNI AZN = ZABS(ZNR, ZNI) RAZN = 1# / AZN STR = ZNR * RAZN STI = -ZNI * RAZN RZR = (STR + STR) * RAZN RZI = (STI + STI) * RAZN FN = FNU + 1# CKR = FN * RZR CKI = FN * RZI '----------------------------------------------------------------------- ' SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS '----------------------------------------------------------------------' CSCL = 1# / TOL CSCR = TOL CSSR(1) = CSCL CSSR(2) = CONER CSSR(3) = CSCR CSRR(1) = CSCR CSRR(2) = CONER CSRR(3) = CSCL BRY(1) = ASCLE BRY(2) = 1# / ASCLE BRY(3) = D1MACH(2) AS2 = ZABS(S2R, S2I) KFLAG = 2 If AS2 > BRY(1) Then GoTo 50 KFLAG = 1 GoTo 60 50 If AS2 < BRY(2) Then GoTo 60 KFLAG = 3 60 BSCLE = BRY(KFLAG) S1R = S1R * CSSR(KFLAG) S1I = S1I * CSSR(KFLAG) S2R = S2R * CSSR(KFLAG) S2I = S2I * CSSR(KFLAG) CSR = CSRR(KFLAG) For I = 3 To N STR = S2R STI = S2I S2R = CKR * STR - CKI * STI + S1R S2I = CKR * STI + CKI * STR + S1I S1R = STR S1I = STI C1R = S2R * CSR C1I = S2I * CSR STR = C1R STI = C1I C2R = YR(I) C2I = YI(I) If KODE = 1 Then GoTo 70 If IUF < 0 Then GoTo 70 ZS1S2 ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF NZ = NZ + NW SC1R = SC2R SC1I = SC2I SC2R = C1R SC2I = C1I If IUF <> 3 Then GoTo 70 IUF = -4 S1R = SC1R * CSSR(KFLAG) S1I = SC1I * CSSR(KFLAG) S2R = SC2R * CSSR(KFLAG) S2I = SC2I * CSSR(KFLAG) STR = SC2R STI = SC2I 70 PTR = CSPNR * C1R - CSPNI * C1I PTI = CSPNR * C1I + CSPNI * C1R YR(I) = PTR + CSGNR * C2R - CSGNI * C2I YI(I) = PTI + CSGNR * C2I + CSGNI * C2R CKR = CKR + RZR CKI = CKI + RZI CSPNR = -CSPNR CSPNI = -CSPNI If KFLAG >= 3 Then GoTo 80 PTR = Abs(C1R) PTI = Abs(C1I) C1M = DMAX(PTR, PTI) If C1M <= BSCLE Then GoTo 80 KFLAG = KFLAG + 1 BSCLE = BRY(KFLAG) S1R = S1R * CSR S1I = S1I * CSR S2R = STR S2I = STI S1R = S1R * CSSR(KFLAG) S1I = S1I * CSSR(KFLAG) S2R = S2R * CSSR(KFLAG) S2I = S2I * CSSR(KFLAG) CSR = CSRR(KFLAG) 80 Next I GoTo 100 90 NZ = -1 If NW = -2 Then NZ = -2 100 End Sub 'ZACON 'end of file Cbess3.bas