Attribute VB_Name = "Module6" '************************************************************************** '* 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/15/2005). * '* (www.jpmoreau.fr) * '************************************************************************** 'MODULE CBESS2 DefDbl A-H, O-Z DefInt I-N 'LIST OF SUBROUTINES: ' ZBUNI(ZR, ZI, FNU,KODE, N, YR(), YI(), NZ, NUI, NLAST, FNUL, TOL, ELIM, ALIM) ' ZUNI1(ZR, ZI, FNU, KODE, N, YR(), YI(), NZ, NLAST, FNUL, TOL, ELIM, ALIM) ' ZUNI2(ZR, ZI, FNU, KODE, N, YR(), YI(), NZ, NLAST, FNUL, TOL, ELIM, ALIM) Sub ZBUNI(ZR, ZI, FNU, KODE, N, YR(), YI(), NZ, NUI, NLAST, FNUL, TOL, ELIM, ALIM) '***BEGIN PROLOGUE ZBUNI '***REFER TO ZBESI,ZBESK ' ' ZBUNI COMPUTES THE I BESSEL FUNCTION FOR LARGE CABS(Z) > ' FNUL AND FNU+N-1 < FNUL. THE ORDER IS INCREASED FROM ' FNU+N-1 GREATER THAN FNUL BY ADDING NUI AND COMPUTING ' ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR I(FNU,Z) ' ON IFORM=1 AND THE EXPANSION FOR J(FNU,Z) ON IFORM=2. ' '***ROUTINES CALLED ZUNI1,ZUNI2,ZABS,D1MACH '***END PROLOGUE ZBUNI ' COMPLEX CSCL,CSCR,Y,RZ,ST,S1,S2,Y,Z 'Labels: 10,20,21,25,30,40,50,60,70,80,90, 100 ' AX, AY, CSCLR, CSCRR, DFNU, FNUI, GNU, RAZ, RZI, RZR, STI, STR, ' S1I, S1R, S2I, S2R, ASCLE, C1R, C1I, C1M:Double ' I, IFLAG, IFORM, K, NL, NW:integer Dim BRY(3) Dim STR As Double NZ = 0 AX = Abs(ZR) * 1.7321 AY = Abs(ZI) IFORM = 1 If AY > AX Then IFORM = 2 If NUI = 0 Then GoTo 60 FNUI = 1# * NUI DFNU = FNU + 1# * (N - 1) GNU = DFNU + FNUI If IFORM = 2 Then GoTo 10 '----------------------------------------------------------------------- ' ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN ' -PI/3 <= ARG(Z) <= PI/3 '----------------------------------------------------------------------- ZUNI1 ZR, ZI, GNU, KODE, 2, YR(), YI(), NW, NLAST, FNUL, TOL, ELIM, ALIM GoTo 20 '----------------------------------------------------------------------- ' ASYMPTOTIC EXPANSION FOR J(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 ZUNI2 ZR, ZI, GNU, KODE, 2, YR(), YI(), NW, NLAST, FNUL, TOL, ELIM, ALIM 20 If NW < 0 Then GoTo 50 If NW <> 0 Then GoTo 90 STR = ZABS(YR(1), YI(1)) '---------------------------------------------------------------------- ' SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED '---------------------------------------------------------------------- BRY(1) = 1000 * D1MACH(1) / TOL BRY(2) = 1# / BRY(1) BRY(3) = BRY(2) IFLAG = 2 ASCLE = BRY(2) CSCLR = 1# If STR > BRY(1) Then GoTo 21 IFLAG = 1 ASCLE = BRY(1) CSCLR = 1# / TOL GoTo 25 21 If STR < BRY(2) Then GoTo 25 IFLAG = 3 ASCLE = BRY(3) CSCLR = TOL 25 CSCRR = 1# / CSCLR S1R = YR(2) * CSCLR S1I = YI(2) * CSCLR S2R = YR(1) * CSCLR S2I = YI(1) * CSCLR RAZ = 1# / ZABS(ZR, ZI) STR = ZR * RAZ STI = -ZI * RAZ RZR = (STR + STR) * RAZ RZI = (STI + STI) * RAZ For I = 1 To NUI STR = S2R STI = S2I S2R = (DFNU + FNUI) * (RZR * STR - RZI * STI) + S1R S2I = (DFNU + FNUI) * (RZR * STI + RZI * STR) + S1I S1R = STR S1I = STI FNUI = FNUI - 1# If IFLAG >= 3 Then GoTo 30 STR = S2R * CSCRR STI = S2I * CSCRR C1R = Abs(STR) C1I = Abs(STI) C1M = DMAX(C1R, C1I) If C1M <= ASCLE Then GoTo 30 IFLAG = IFLAG + 1 ASCLE = BRY(IFLAG) S1R = S1R * CSCRR S1I = S1I * CSCRR S2R = STR S2I = STI CSCLR = CSCLR * TOL CSCRR = 1# / CSCLR S1R = S1R * CSCLR S1I = S1I * CSCLR S2R = S2R * CSCLR S2I = S2I * CSCLR 30 Next I YR(N) = S2R * CSCRR YI(N) = S2I * CSCRR If N = 1 Then GoTo 100 NL = N - 1 FNUI = 1# * NL K = NL For I = 1 To NL STR = S2R STI = S2I S2R = (FNU + FNUI) * (RZR * STR - RZI * STI) + S1R S2I = (FNU + FNUI) * (RZR * STI + RZI * STR) + S1I S1R = STR S1I = STI STR = S2R * CSCRR STI = S2I * CSCRR YR(K) = STR YI(K) = STI FNUI = FNUI - 1# K = K - 1 If IFLAG >= 3 Then GoTo 40 C1R = Abs(STR) C1I = Abs(STI) C1M = DMAX(C1R, C1I) If C1M <= ASCLE Then GoTo 40 IFLAG = IFLAG + 1 ASCLE = BRY(IFLAG) S1R = S1R * CSCRR S1I = S1I * CSCRR S2R = STR S2I = STI CSCLR = CSCLR * TOL CSCRR = 1# / CSCLR S1R = S1R * CSCLR S1I = S1I * CSCLR S2R = S2R * CSCLR S2I = S2I * CSCLR 40 Next I GoTo 100 50 NZ = -1 If NW = -2 Then NZ = -2 GoTo 100 60 If IFORM = 2 Then GoTo 70 '----------------------------------------------------------------------- ' ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN ' -PI/3 <= ARG(Z) <= PI/3 '----------------------------------------------------------------------- ZUNI1 ZR, ZI, FNU, KODE, N, YR(), YI(), NW, NLAST, FNUL, TOL, ELIM, ALIM GoTo 80 '----------------------------------------------------------------------- ' ASYMPTOTIC EXPANSION FOR J(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 '----------------------------------------------------------------------- 70 ZUNI2 ZR, ZI, FNU, KODE, N, YR(), YI(), NW, NLAST, FNUL, TOL, ELIM, ALIM 80 If NW < 0 Then GoTo 50 NZ = NW GoTo 100 90 NLAST = N 100 End Sub 'ZBUNI Sub ZUNI1(ZR, ZI, FNU, KODE, N, YR(), YI(), NZ, NLAST, FNUL, TOL, ELIM, ALIM) '***BEGIN PROLOGUE ZUNI1 '***REFER TO ZBESI,ZBESK ' ' ZUNI1 COMPUTES I(FNU,Z) BY MEANS OF THE UNIFORM ASYMPTOTIC ' EXPANSION FOR I(FNU,Z) IN -PI/3 <= ARG Z <= PI/3. ' ' FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC ' EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET. ' NLAST <> 0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER ' FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1 < FNUL. ' Y(I)=CZERO FOR I=NLAST+1,N. ' '***ROUTINES CALLED ZUCHK,ZUNIK,ZUOIK,D1MACH,ZABS '***END PROLOGUE ZUNI1 ' COMPLEX CFN,CONE,CRSC,CSCL,CSR,CSS,CWRK,CZERO,C1,C2,PHI,RZ,SUM,S1, ' S2,Y,Z,ZETA1,ZETA2 'Labels: 10,20,30,40,50,60,70,90,100,110,120,130, 150 ' APHI, ASCLE, CONEI, CONER, CRSC, CSCL, C1R, C2I, C2M, C2R, FN, ' PHII, PHIR, RAST, RS1, RZI, RZR, STI, STR, SUMI, SUMR, S1I, S1R, ' S2I, S2R, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R: Double ' I, IFLAG, INIT, K, M, ND, NN, NUF, NW: Integer Dim BRY(3), CSSR(3), CSRR(3) Dim STR As Double Dim CWRKR(16), CWRKI(16) ZEROR = 0#: ZEROI = 0#: CONER = 1#: CONEI = 0# NZ = 0 ND = N NLAST = 0 '----------------------------------------------------------------------- ' COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG- ' NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE, ' Exp(ALIM) = Exp(ELIM) * TOL '----------------------------------------------------------------------- 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 '----------------------------------------------------------------------- ' CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER '----------------------------------------------------------------------- FN = DMAX(FNU, 1#) INIT = 0 ZUNIK ZR, ZI, FN, 1, 1, TOL, INIT, PHIR, PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR(), CWRKI() If KODE = 1 Then GoTo 10 STR = ZR + ZETA2R STI = ZI + ZETA2I RAST = FN / ZABS(STR, STI) STR = STR * RAST * RAST STI = -STI * RAST * RAST S1R = -ZETA1R + STR S1I = -ZETA1I + STI GoTo 20 10 S1R = -ZETA1R + ZETA2R S1I = -ZETA1I + ZETA2I 20 RS1 = S1R If Abs(RS1) > ELIM Then GoTo 130 30 NN = IMIN(2, ND) For I = 1 To NN FN = FNU + 1# * (ND - I) INIT = 0 ZUNIK ZR, ZI, FN, 1, 0, TOL, INIT, PHIR, PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR(), CWRKI() If KODE = 1 Then GoTo 40 STR = ZR + ZETA2R STI = ZI + ZETA2I RAST = FN / ZABS(STR, STI) STR = STR * RAST * RAST STI = -STI * RAST * RAST S1R = -ZETA1R + STR S1I = -ZETA1I + STI + ZI GoTo 50 40 S1R = -ZETA1R + ZETA2R S1I = -ZETA1I + ZETA2I '----------------------------------------------------------------------- ' TEST FOR UNDERFLOW AND OVERFLOW '----------------------------------------------------------------------- 50 RS1 = S1R If Abs(RS1) > ELIM Then GoTo 110 If I = 1 Then IFLAG = 2 If Abs(RS1) < ALIM Then GoTo 60 '----------------------------------------------------------------------- ' REFINE TEST AND SCALE '----------------------------------------------------------------------- APHI = ZABS(PHIR, PHII) RS1 = RS1 + Log(APHI) If Abs(RS1) > ELIM Then GoTo 110 If I = 1 Then IFLAG = 1 If RS1 < 0# Then GoTo 60 If I = 1 Then IFLAG = 3 '----------------------------------------------------------------------- ' SCALE S1 IF CABS(S1) < ASCLE '----------------------------------------------------------------------- 60 S2R = PHIR * SUMR - PHII * SUMI S2I = PHIR * SUMI + PHII * SUMR STR = Exp(S1R) * CSSR(IFLAG) S1R = STR * Cos(S1I) S1I = STR * Sin(S1I) STR = S2R * S1R - S2I * S1I S2I = S2R * S1I + S2I * S1R S2R = STR If IFLAG <> 1 Then GoTo 70 ZUCHK S2R, S2I, NW, BRY(1), TOL If NW <> 0 Then GoTo 110 70: YR(I) = S2R YI(I) = S2I M = ND - I + 1 YR(M) = S2R * CSRR(IFLAG) YI(M) = S2I * CSRR(IFLAG) Next I If ND <= 2 Then GoTo 100 RAST = 1# / ZABS(ZR, ZI) STR = ZR * RAST STI = -ZI * RAST RZR = (STR + STR) * RAST RZI = (STI + STI) * RAST BRY(2) = 1# / BRY(1) BRY(3) = D1MACH(2) S1R = YR(1) S1I = YI(1) S2R = YR(2) S2I = YI(2) C1R = CSRR(IFLAG) ASCLE = BRY(IFLAG) K = ND - 2 FN = 1# * K For I = 3 To ND C2R = S2R C2I = S2I S2R = S1R + (FNU + FN) * (RZR * C2R - RZI * C2I) S2I = S1I + (FNU + FN) * (RZR * C2I + RZI * C2R) S1R = C2R S1I = C2I C2R = S2R * C1R C2I = S2I * C1R YR(K) = C2R YI(K) = C2I K = K - 1 FN = FN - 1# If IFLAG >= 3 Then GoTo 90 STR = Abs(C2R) STI = Abs(C2I) C2M = DMAX(STR, STI) If C2M <= ASCLE Then GoTo 90 IFLAG = IFLAG + 1 ASCLE = BRY(IFLAG) S1R = S1R * C1R S1I = S1I * C1R S2R = C2R S2I = C2I S1R = S1R * CSSR(IFLAG) S1I = S1I * CSSR(IFLAG) S2R = S2R * CSSR(IFLAG) S2I = S2I * CSSR(IFLAG) C1R = CSRR(IFLAG) 90 Next I 100 GoTo 150 '----------------------------------------------------------------------- ' SET UNDERFLOW AND UPDATE PARAMETERS '----------------------------------------------------------------------- 110 If RS1 > 0# Then GoTo 120 YR(ND) = ZEROR YI(ND) = ZEROI NZ = NZ + 1 ND = ND - 1 If ND = 0 Then GoTo 100 ZUOIK ZR, ZI, FNU, KODE, 1, ND, YR(), YI(), NUF, TOL, ELIM, ALIM If NUF < 0 Then GoTo 120 ND = ND - NUF NZ = NZ + NUF If ND = 0 Then GoTo 100 FN = FNU + 1# * (ND - 1) If FN >= FNUL Then GoTo 30 NLAST = ND GoTo 150 120 NZ = -1 GoTo 150 130 If RS1 > 0# Then GoTo 120 NZ = N For I = 1 To N YR(I) = ZEROR YI(I) = ZEROI Next I 150 End Sub 'ZUNI1 Sub ZUNI2(ZR, ZI, FNU, KODE, N, YR(), YI(), NZ, NLAST, FNUL, TOL, ELIM, ALIM) '***BEGIN PROLOGUE ZUNI2 '***REFER TO ZBESI,ZBESK ' ' ZUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF ' UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I ' OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO. ' ' FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC ' EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET. ' NLAST <> 0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER ' FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1 < FNUL. ' Y(I)=CZERO FOR I=NLAST+1,N. ' '***ROUTINES CALLED ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,ZABS '***END PROLOGUE ZUNI2 ' COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS, ' CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN 'Note: variable IN ==> IN0 (IN is a reserved word in Pascal). 'Labels: 10,20,30,40,50,60,70,80, 100,110,120,130,140,150, 200 ' AARG, AIC, AII, AIR, ANG, APHI, ARGI, ARGR, ASCLE, ASUMI, ASUMR, ' BSUMI, BSUMR, CIDI, CONEI, CONER, CRSC, CSCL, C1R, C2I, C2M, C2R, ' DAII, DAIR, FN, HPI, PHII, PHIR, RAST, RAZ, RS1, RZI, RZR, STI, ' STR, S1I, S1R, S2I, S2R, ZBI, ZBR, ZEROI, ZEROR, ZETA1I, ZETA1R, ' ZETA2I, ZETA2R, ZNI, ZNR: Double ' I, IFLAG, IN0, INU, J, K, NAI, ND, NDAI, NN, NUF, NW, IDUM: Integer Dim BRY(3), CIPR(4), CIPI(4), CSSR(3), CSRR(3) Dim STR As Double ZEROR = 0#: ZEROI = 0#: CONER = 1#: CONEI = 0# CIPR(1) = 1#: CIPI(1) = 0#: CIPR(2) = 0#: CIPI(2) = 1# CIPR(3) = -1#: CIPI(3) = 0#: CIPR(4) = 0#: CIPI(4) = -1# HPI = 1.5707963267949: AIC = 1.26551212348465 NZ = 0 ND = N NLAST = 0 '----------------------------------------------------------------------- ' COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG- ' NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE, ' Exp(ALIM) = Exp(ELIM) * TOL '----------------------------------------------------------------------- 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 '---------------------------------------------------------------------- ' ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI '----------------------------------------------------------------------- ZNR = ZI ZNI = -ZR ZBR = ZR ZBI = ZI CIDI = -CONER INU = Int(FNU) ANG = HPI * (FNU - 1# * INU) C2R = Cos(ANG) C2I = Sin(ANG) IN0 = INU + N - 1 IN0 = (IN0 Mod 4) + 1 STR = C2R * CIPR(IN0) - C2I * CIPI(IN0) C2I = C2R * CIPI(IN0) + C2I * CIPR(IN0) C2R = STR If ZI > 0# Then GoTo 10 ZNR = -ZNR ZBI = -ZBI CIDI = -CIDI C2I = -C2I '----------------------------------------------------------------------- ' CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER '----------------------------------------------------------------------- 10 FN = DMAX(FNU, 1#) ZUNHJ ZNR, ZNI, FN, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI If KODE = 1 Then GoTo 20 STR = ZBR + ZETA2R STI = ZBI + ZETA2I RAST = FN / ZABS(STR, STI) STR = STR * RAST * RAST STI = -STI * RAST * RAST S1R = -ZETA1R + STR S1I = -ZETA1I + STI GoTo 30 20 S1R = -ZETA1R + ZETA2R S1I = -ZETA1I + ZETA2I 30 RS1 = S1R If Abs(RS1) > ELIM Then GoTo 150 40 NN = IMIN(2, ND) For I = 1 To NN FN = FNU + 1# * (ND - I) ZUNHJ ZNR, ZNI, FN, 0, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI If KODE = 1 Then GoTo 50 STR = ZBR + ZETA2R STI = ZBI + ZETA2I RAST = FN / ZABS(STR, STI) STR = STR * RAST * RAST STI = -STI * RAST * RAST S1R = -ZETA1R + STR S1I = -ZETA1I + STI + Abs(ZI) GoTo 60 50 S1R = -ZETA1R + ZETA2R S1I = -ZETA1I + ZETA2I '----------------------------------------------------------------------- ' TEST FOR UNDERFLOW AND OVERFLOW '----------------------------------------------------------------------- 60 RS1 = S1R If Abs(RS1) > ELIM Then GoTo 120 If I = 1 Then IFLAG = 2 If Abs(RS1) < ALIM Then GoTo 70 '----------------------------------------------------------------------- ' REFINE TEST AND SCALE '----------------------------------------------------------------------- APHI = ZABS(PHIR, PHII) AARG = ZABS(ARGR, ARGI) RS1 = RS1 + Log(APHI) - 0.25 * Log(AARG) - AIC If Abs(RS1) > ELIM Then GoTo 120 If I = 1 Then IFLAG = 1 If RS1 < 0# Then GoTo 70 If I = 1 Then IFLAG = 3 '----------------------------------------------------------------------- ' SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR ' EXPONENT EXTREMES '----------------------------------------------------------------------- 70 ZAIRY ARGR, ARGI, 0, 2, AIR, AII, NAI, IDUM ZAIRY ARGR, ARGI, 1, 2, DAIR, DAII, NDAI, IDUM STR = DAIR * BSUMR - DAII * BSUMI STI = DAIR * BSUMI + DAII * BSUMR STR = STR + (AIR * ASUMR - AII * ASUMI) STI = STI + (AIR * ASUMI + AII * ASUMR) S2R = PHIR * STR - PHII * STI S2I = PHIR * STI + PHII * STR STR = Exp(S1R) * CSSR(IFLAG) S1R = STR * Cos(S1I) S1I = STR * Sin(S1I) STR = S2R * S1R - S2I * S1I S2I = S2R * S1I + S2I * S1R S2R = STR If IFLAG <> 1 Then GoTo 80 ZUCHK S2R, S2I, NW, BRY(1), TOL If NW <> 0 Then GoTo 120 80 If ZI <= 0# Then S2I = -S2I STR = S2R * C2R - S2I * C2I S2I = S2R * C2I + S2I * C2R S2R = STR YR(I) = S2R YI(I) = S2I J = ND - I + 1 YR(J) = S2R * CSRR(IFLAG) YI(J) = S2I * CSRR(IFLAG) STR = -C2I * CIDI C2I = C2R * CIDI C2R = STR Next I If ND <= 2 Then GoTo 110 RAZ = 1# / ZABS(ZR, ZI) STR = ZR * RAZ STI = -ZI * RAZ RZR = (STR + STR) * RAZ RZI = (STI + STI) * RAZ BRY(2) = 1# / BRY(1) BRY(3) = D1MACH(2) S1R = YR(1) S1I = YI(1) S2R = YR(2) S2I = YI(2) C1R = CSRR(IFLAG) ASCLE = BRY(IFLAG) K = ND - 2 FN = 1# * K For I = 3 To ND C2R = S2R C2I = S2I S2R = S1R + (FNU + FN) * (RZR * C2R - RZI * C2I) S2I = S1I + (FNU + FN) * (RZR * C2I + RZI * C2R) S1R = C2R S1I = C2I C2R = S2R * C1R C2I = S2I * C1R YR(K) = C2R YI(K) = C2I K = K - 1 FN = FN - 1# If IFLAG >= 3 Then GoTo 100 STR = Abs(C2R) STI = Abs(C2I) C2M = DMAX(STR, STI) If C2M <= ASCLE Then GoTo 100 IFLAG = IFLAG + 1 ASCLE = BRY(IFLAG) S1R = S1R * C1R S1I = S1I * C1R S2R = C2R S2I = C2I S1R = S1R * CSSR(IFLAG) S1I = S1I * CSSR(IFLAG) S2R = S2R * CSSR(IFLAG) S2I = S2I * CSSR(IFLAG) C1R = CSRR(IFLAG) 100 Next I 110 GoTo 200 120 If RS1 > 0# Then GoTo 140 '----------------------------------------------------------------------- ' SET UNDERFLOW AND UPDATE PARAMETERS '----------------------------------------------------------------------- YR(ND) = ZEROR YI(ND) = ZEROI NZ = NZ + 1 ND = ND - 1 If ND = 0 Then GoTo 110 ZUOIK ZR, ZI, FNU, KODE, 1, ND, YR(), YI(), NUF, TOL, ELIM, ALIM If NUF < 0 Then GoTo 140 ND = ND - NUF NZ = NZ + NUF If ND = 0 Then GoTo 110 FN = FNU + 1# * (ND - 1) If FN < FNUL Then GoTo 130 FN = CIDI J = NUF + 1 K = (J Mod 4) + 1 S1R = CIPR(K) S1I = CIPI(K) If FN < 0# Then S1I = -S1I STR = C2R * S1R - C2I * S1I C2I = C2R * S1I + C2I * S1R C2R = STR GoTo 40 130 NLAST = ND GoTo 200 140 NZ = -1 GoTo 200 150: If RS1 < 0# Then GoTo 140 NZ = N For I = 1 To N YR(I) = ZEROR YI(I) = ZEROI Next I 200 End Sub 'ZUNI2 'end of file cbess2.bas