Attribute VB_Name = "Module1" DefDbl A-H, O-Z DefInt I-N '***************************************************************** '* EVALUATE A Y-BESSEL FUNCTION OF COMPLEX ARGUMENT (SECOND KIND)* '* ------------------------------------------------------------- * '* SAMPLE RUN: * '* (Evaluate Y0 to Y4 for argument Z=(1.0,2.0) ). * '* * '* cyr(0) = 1.367418... * '* cyi(0) = 1.521506... * '* cyr(1) = -1.089469... * '* cyi(1) = 1.314951... * '* cyr(2) = -0.751245... * '* cyi(2) = -0.123950... * '* cyr(3) = 0.290153... * '* cyi(3) = -0.212118... * '* cyr(4) = 0.590344... * '* cyi(4) = -0.826960... * '* NZ = 0 * '* Error code: 0 * '* * '* ------------------------------------------------------------- * '* Ref.: From Numath Library By Tuan Dang Trong in Fortran 77 * '* [BIBLI 18]. * '* * '* Visual Basic Release 1.0 By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '***************************************************************** 'Note: Used files: CBess0,CBess00,CBess1,CBess2,CBess3,Complex. '-------------------------------------------------------------- 'PROGRAM TEST_ZBESY 'Global variables ' NONE Sub ZBESY(ZR, ZI, FNU, KODE, N, cyr(), cyi(), NZ, IERR) '***BEGIN PROLOGUE ZBESY '***DATE WRITTEN 830501 (YYMMDD) (Original Fortran Version). '***REVISION DATE 830501 (YYMMDD) '***CATEGORY NO. B5K '***KEYWORDS Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT, ' BESSEL FUNCTION OF SECOND KIND '***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES '***PURPOSE TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT '***DESCRIPTION ' ' ***A DOUBLE PRECISION ROUTINE*** ' ' ON KODE=1, CBESY COMPUTES AN N MEMBER SEQUENCE OF COMPLEX ' BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL, NONNEGATIVE ' ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE ' -PI < ARG(Z) <= PI. ON KODE=2, CBESY RETURNS THE SCALED ' FUNCTIONS ' ' CY(I)=EXP(-ABS(Y))*Y(FNU+I-1,Z) I = 1,...,N , Y=AIMAG(Z) ' ' WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND ' LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION ' ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS ' (REF. 1). ' ' INPUT ZR,ZI,FNU ARE DOUBLE PRECISION ' ZR,ZI - Z=CMPLX(ZR,ZI), Z <> CMPLX(0.0,0.0), ' -PI < ARG(Z) <= PI ' FNU - ORDER OF INITIAL Y FUNCTION, FNU >= 0.0 ' KODE - A PARAMETER TO INDICATE THE SCALING OPTION ' KODE= 1 RETURNS ' CY(I)=Y(FNU+I-1,Z), I=1,...,N ' = 2 RETURNS ' CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N ' WHERE Y=AIMAG(Z) ' N - NUMBER OF MEMBERS OF THE SEQUENCE, N >= 1 ' CWRKR, - DOUBLE PRECISION WORK VECTORS OF DIMENSION AT ' CWRKI LEAST N ' ' OUTPUT CYR,CYI ARE DOUBLE PRECISION ' CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS ' CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE ' CY(I)=Y(FNU+I-1,Z) OR ' CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)) I=1,...,N ' DEPENDING ON KODE. ' NZ - NZ=0 , A NORMAL RETURN ' NZ > 0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO ' UNDERFLOW (GENERALLY ON KODE=2) ' IERR - ERROR FLAG ' IERR=0, NORMAL RETURN - COMPUTATION COMPLETED ' IERR=1, INPUT ERROR - NO COMPUTATION ' IERR=2, OVERFLOW - NO COMPUTATION, FNU IS ' TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH ' IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE ' BUT LOSSES OF SIGNIFCANCE BY ARGUMENT ' REDUCTION PRODUCE LESS THAN HALF OF MACHINE ' ACCURACY ' IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- ' TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- ' CANCE BY ARGUMENT REDUCTION ' IERR=5, ERROR - NO COMPUTATION, ' ALGORITHM TERMINATION CONDITION NOT MET ' '***LONG DESCRIPTION ' ' THE COMPUTATION IS CARRIED OUT BY THE FORMULA ' 'Y(FNU, Z) = 0.5 * (H(1, FNU, Z) - H(2, FNU, Z)) / I ' ' WHERE I^2 = -1 AND THE HANKEL BESSEL FUNCTIONS H(1,FNU,Z) ' AND H(2,FNU,Z) ARE CALCULATED IN CBESH. ' ' FOR NEGATIVE ORDERS,THE FORMULA ' 'Y(-FNU, Z) = Y(FNU, Z) * Cos(PI * FNU) + J(FNU, Z) * Sin(PI * FNU) ' ' CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD ' INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE ' POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)* ' SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS ' NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A ' LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM ' CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, ' WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF ' ODD INTEGER. HERE, LARGE MEANS FNU > CABS(Z). ' ' IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- ' MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS ' LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. ' CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN ' LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG ' IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS ' DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. ' IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS ' LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS ' MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE ' INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS ' 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 LIMITING IN ' THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT ' TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS ' IN SINGLE 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,1.0E-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 ' BY D. E. AMOS, SAND83-0083, MAY, 1983. ' ' 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 ZBESH,I1MACH,D1MACH '***END PROLOGUE ZBESY ' ' COMPLEX CWRK,CY,C1,C2,EX,HCI,Z 'Labels: 60,70,90,170, 200 ' C1I, C1R, C2I, C2R, ELIM, EXI, EXR, EY, HCII, STI, STR, TAY: double ' I, K, K1, K2, NZ1, NZ2: integer Dim STR As Double Dim CWRKR(10), CWRKI(10) '***FIRST EXECUTABLE STATEMENT ZBESY IERR = 0 NZ = 0 If (ZR = 0#) And (ZI = 0#) Then IERR = 1 If FNU < 0# Then IERR = 1 If (KODE < 1) Or (KODE > 2) Then IERR = 1 If N < 1 Then IERR = 1 If IERR <> 0 Then GoTo 200 'bad parameter(s) HCII = 0.5 ZBESH ZR, ZI, FNU, KODE, 1, N, cyr(), cyi(), NZ1, IERR If (IERR <> 0) And (IERR <> 3) Then GoTo 170 ZBESH ZR, ZI, FNU, KODE, 2, N, CWRKR(), CWRKI(), NZ2, IERR If (IERR <> 0) And (IERR <> 3) Then GoTo 170 NZ = IMIN(NZ1, NZ2) If KODE = 2 Then GoTo 60 For I = 1 To N STR = CWRKR(I) - cyr(I) STI = CWRKI(I) - cyi(I) cyr(I) = -STI * HCII cyi(I) = STR * HCII Next I GoTo 200 60 K1 = Int(XI1MACH(15)) K2 = Int(XI1MACH(16)) K = IMIN(Abs(K1), Abs(K2)) '----------------------------------------------------------------------- ' ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT '----------------------------------------------------------------------- ELIM = 2.303 * (K * D1MACH(5) - 3#) EXR = Cos(ZR) EXI = Sin(ZR) EY = 0# TAY = Abs(ZI + ZI) If TAY < ELIM Then EY = Exp(-TAY) If ZI < 0# Then GoTo 90 C1R = EXR * EY C1I = EXI * EY C2R = EXR C2I = -EXI 70 NZ = 0 For I = 1 To N STR = C1R * cyr(I) - C1I * cyi(I) STI = C1R * cyi(I) + C1I * cyr(I) STR = -STR + C2R * CWRKR(I) - C2I * CWRKI(I) STI = -STI + C2R * CWRKI(I) + C2I * CWRKR(I) cyr(I) = -STI * HCII cyi(I) = STR * HCII If (STR = 0#) And (STI = 0#) And (EY = 0#) Then NZ = NZ + 1 Next I GoTo 200 90 C1R = EXR C1I = EXI C2R = EXR * EY C2I = -EXI * EY GoTo 70 170 NZ = 0 200 End Sub 'ZBESY Sub ZBESH(ZR, ZI, FNU, KODE, M, N, cyr(), cyi(), NZ, IERR) '***BEGIN PROLOGUE ZBESH '***DATE WRITTEN 830501 (YYMMDD) '***REVISION DATE 830501 (YYMMDD) '***CATEGORY NO. B5K '***KEYWORDS H-BESSEL FUNCTIONS,BESSEL FUNCTIONS OF COMPLEX ARGUMENT, ' BESSEL FUNCTIONS OF THIRD KIND,HANKEL FUNCTIONS '***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES '***PURPOSE TO COMPUTE THE H-BESSEL FUNCTIONS OF A COMPLEX ARGUMENT '***DESCRIPTION ' ' ***A DOUBLE PRECISION ROUTINE*** ' ON KODE=1, ZBESH COMPUTES AN N MEMBER SEQUENCE OF COMPLEX ' HANKEL (BESSEL) FUNCTIONS CY(J)=H(M,FNU+J-1,Z) FOR KINDS M=1 ' OR 2, REAL, NONNEGATIVE ORDERS FNU+J-1, J=1,...,N, AND COMPLEX ' Z <> CMPLX(0.0,0.0) IN THE CUT PLANE -PI < ARG(Z) <= PI. ' ON KODE=2, ZBESH RETURNS THE SCALED HANKEL FUNCTIONS ' ' CY(I)=EXP(-MM*Z*I)*H(M,FNU+J-1,Z) MM=3-2*M, I^2=-1. ' ' WHICH REMOVES THE EXPONENTIAL BEHAVIOR IN BOTH THE UPPER AND ' LOWER HALF PLANES. DEFINITIONS AND NOTATION ARE FOUND IN THE ' NBS HANDBOOK OF MATHEMATICAL FUNCTIONS (REF. 1). ' ' INPUT ZR,ZI,FNU ARE DOUBLE PRECISION ' ZR,ZI - Z=CMPLX(ZR,ZI), Z <> CMPLX(0.0,0.0), ' -PI < ARG(Z) <= PI ' FNU - ORDER OF INITIAL H FUNCTION, FNU >= 0.0 ' KODE - A PARAMETER TO INDICATE THE SCALING OPTION ' KODE= 1 RETURNS ' CY(J)=H(M,FNU+J-1,Z), J=1,...,N ' = 2 RETURNS ' CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M)) ' J=1,...,N , I^2=-1 ' M - KIND OF HANKEL FUNCTION, M=1 OR 2 ' N - NUMBER OF MEMBERS IN THE SEQUENCE, N >= 1 ' ' OUTPUT CYR,CYI ARE DOUBLE PRECISION ' CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS ' CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE ' CY(J)=H(M,FNU+J-1,Z) OR ' CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M)) J=1,...,N ' DEPENDING ON KODE, I^2=-1. ' NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW, ' NZ= 0, NORMAL RETURN ' NZ > 0, FIRST NZ COMPONENTS OF CY SET TO ZERO DUE ' TO UNDERFLOW, CY(J)=CMPLX(0.0,0.0) ' J=1,...,NZ WHEN Y > 0.0 AND M=1 OR ' Y < 0.0 AND M=2. FOR THE COMPLMENTARY ' HALF PLANES, NZ STATES ONLY THE NUMBER ' OF UNDERFLOWS. ' IERR - ERROR FLAG ' IERR=0, NORMAL RETURN - COMPUTATION COMPLETED ' IERR=1, INPUT ERROR - NO COMPUTATION ' IERR=2, OVERFLOW - NO COMPUTATION, FNU TOO ' LARGE OR CABS(Z) TOO SMALL OR BOTH ' IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE ' BUT LOSSES OF SIGNIFCANCE BY ARGUMENT ' REDUCTION PRODUCE LESS THAN HALF OF MACHINE ' ACCURACY ' IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- ' TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- ' CANCE BY ARGUMENT REDUCTION ' IERR=5, ERROR - NO COMPUTATION, ' ALGORITHM TERMINATION CONDITION NOT MET ' '***LONG DESCRIPTION ' ' THE COMPUTATION IS CARRIED OUT BY THE RELATION ' ' H(M,FNU,Z)=(1/MP)*EXP(-MP*FNU)*K(FNU,Z*EXP(-MP)) ' MP=MM*HPI*I, MM=3-2*M, HPI=PI/2, I^2=-1 ' ' FOR M=1 OR 2 WHERE THE K BESSEL FUNCTION IS COMPUTED FOR THE ' RIGHT HALF PLANE RE(Z) >= 0. THE K FUNCTION IS CONTINUED ' TO THE LEFT HALF PLANE BY THE RELATION ' ' K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z) ' MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I^2=-1 ' ' WHERE I(FNU,Z) IS THE I BESSEL FUNCTION. ' ' EXPONENTIAL DECAY OF H(M,FNU,Z) OCCURS IN THE UPPER HALF Z ' PLANE FOR M=1 AND THE LOWER HALF Z PLANE FOR M=2. EXPONENTIAL ' GROWTH OCCURS IN THE COMPLEMENTARY HALF PLANES. SCALING ' BY EXP(-MM*Z*I) REMOVES THE EXPONENTIAL BEHAVIOR IN THE ' WHOLE Z PLANE FOR Z TO INFINITY. ' ' FOR NEGATIVE ORDERS,THE FORMULAE ' ' H(1,-FNU,Z) = H(1,FNU,Z)*CEXP( PI*FNU*I) ' H(2,-FNU,Z) = H(2,FNU,Z)*CEXP(-PI*FNU*I) ' I^2=-1 ' ' CAN BE USED. ' ' IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- ' MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS ' LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. ' CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN ' LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG ' IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS ' DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. ' IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS ' LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS ' MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE ' INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS ' 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 ' ARITHMETI' AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION ' ARITHMETI' RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN ' THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT ' TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS ' IN SINGLE 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,1.0D-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 ARITHMETI' 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 ' BY D. E. AMOS, SAND83-0083, MAY, 1983. ' ' 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 ZACON,ZBKNU,ZBUNK,ZUOIK,ZABS,I1MACH,D1MACH '***END PROLOGUE ZBESH ' ' COMPLEX CY,Z,ZN,ZT 'Labels: 60,70,80,90,100,110,120,140,230,240,260, 300 Dim STR As Double ' AA, ALIM, ALN, ARG, AZ, DIG, ELIM, ' FMM, FN, FNUL, HPI, RHPI, RL, R1M5, SGN, STR, TOL, UFL, ' ZNI, ZNR, ZTI, BB: Double ' I, INU, INUH, IR, K, K1, K2, MM, MR, NN, NUF, NW: Integer '***FIRST EXECUTABLE STATEMENT ZBESH HPI = 1.5707963267949 IERR = 0 NZ = 0 If (ZR = 0#) And (ZI = 0#) Then IERR = 1 If FNU < 0# Then IERR = 1 If (M < 1) Or (M > 2) Then IERR = 1 If (KODE < 1) Or (KODE > 2) Then IERR = 1 If N < 1 Then IERR = 1 If IERR <> 0 Then GoTo 300 NN = N '----------------------------------------------------------------------- ' SET PARAMETERS RELATED TO MACHINE CONSTANTS. ' TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1E-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). ' FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU '----------------------------------------------------------------------- TOL = DMAX(D1MACH(4), 1E-18) 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 * 1# * K1 DIG = DMIN(AA, 18#) AA = AA * 2.303 ALIM = ELIM + DMAX(-AA, -41.45) FNUL = 10# + 6# * (DIG - 3#) RL = 1.2 * DIG + 3# FN = FNU + 1# * (NN - 1) MM = 3 - M - M FMM = 1# * MM ZNR = FMM * ZI ZNI = -FMM * ZR '----------------------------------------------------------------------- ' TEST FOR PROPER RANGE '----------------------------------------------------------------------- AZ = ZABS(ZR, ZI) AA = 0.5 / TOL BB = 0.5 * XI1MACH(9) AA = DMIN(AA, BB) If AZ > AA Then GoTo 260 If FN > AA Then GoTo 260 AA = Sqr(AA) If AZ > AA Then IERR = 3 If FN > AA Then IERR = 3 '----------------------------------------------------------------------- ' OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE '----------------------------------------------------------------------- UFL = Exp(-ELIM) If AZ < UFL Then GoTo 230 If FNU > FNUL Then GoTo 90 If FN <= 1# Then GoTo 70 If FN > 2# Then GoTo 60 If AZ > TOL Then GoTo 70 ARG = 0.5 * AZ ALN = -FN * Log(ARG) If ALN > ELIM Then GoTo 230 GoTo 70 60 ZUOIK ZNR, ZNI, FNU, KODE, 2, NN, cyr(), cyi(), NUF, TOL, ELIM, ALIM If NUF < 0 Then GoTo 230 NZ = NZ + NUF NN = NN - NUF '----------------------------------------------------------------------- ' HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK ' IF NUF=NN, THEN CY(I)=CZERO FOR ALL I '----------------------------------------------------------------------- If NN = 0 Then GoTo 140 70 If (ZNR < 0#) Or ((ZNR = 0#) And (ZNI < 0#) And (M = 2)) Then GoTo 80 '----------------------------------------------------------------------- ' RIGHT HALF PLANE COMPUTATION, XN >= 0 AND (XN <> 0 OR YN >= 0 OR M=1) '----------------------------------------------------------------------- ZBKNU ZNR, ZNI, FNU, KODE, NN, cyr(), cyi(), NZ, TOL, ELIM, ALIM GoTo 110 '----------------------------------------------------------------------- ' LEFT HALF PLANE COMPUTATION '----------------------------------------------------------------------- 80 MR = -MM ZACON ZNR, ZNI, FNU, KODE, MR, NN, cyr(), cyi(), NW, RL, FNUL, TOL, ELIM, ALIM If NW < 0 Then GoTo 240 NZ = NW GoTo 110 '----------------------------------------------------------------------- ' UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL '----------------------------------------------------------------------- 90 MR = 0 If (ZNR >= 0#) And ((ZNR <> 0#) Or (ZNI >= 0#) Or (M <> 2)) Then GoTo 100 MR = -MM If (ZNR <> 0#) Or (ZNI >= 0#) Then GoTo 100 ZNR = -ZNR ZNI = -ZNI 100 ZBUNK ZNR, ZNI, FNU, KODE, MR, NN, cyr(), cyi(), NW, TOL, ELIM, ALIM If NW < 0 Then GoTo 240 NZ = NZ + NW '----------------------------------------------------------------------- ' H(M,FNU,Z) = -FMM*(I/HPI)*(ZT^FNU)*K(FNU,-Z*ZT) ' ' ZT=EXP(-FMM*HPI*I) = CMPLX(0.0,-FMM), FMM=3-2*M, M=1,2 '----------------------------------------------------------------------- 110 SGN0 = Sign(HPI, -FMM) '----------------------------------------------------------------------- ' CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE ' WHEN FNU IS LARGE '----------------------------------------------------------------------- INU = Int(FNU) INUH = INU / 2 IR = INU - 2 * INUH ARG = (FNU - 1# * (INU - IR)) * SGN0 RHPI = 1# / SGN0 ZNI = RHPI * Cos(ARG) ZNR = -RHPI * Sin(ARG) If (INUH Mod 2) = 0 Then GoTo 120 ZNR = -ZNR ZNI = -ZNI 120 ZTI = -FMM For I = 1 To NN STR = cyr(I) * ZNR - cyi(I) * ZNI cyi(I) = cyr(I) * ZNI + cyi(I) * ZNR cyr(I) = STR STR = -ZNI * ZTI ZNI = ZNR * ZTI ZNR = STR Next I GoTo 300 140 Form1.Print " 140: ZNR="; ZNR If ZNR < 0# Then GoTo 230 GoTo 300 230 NZ = 0 IERR = 2 GoTo 300 240 If NW = -1 Then GoTo 230 NZ = 0 IERR = 5 GoTo 300 260 NZ = 0 IERR = 4 300 End Sub 'ZBESH Sub Exec_Tzbesy() Dim cyr(10), cyi(10) N = 5 ZR = 1#: ZI = 2# Form1.Cls ZBESY ZR, ZI, 0, 1, N, cyr(), cyi(), NZ, IERR Form1.Print For I = 1 To N Form1.Print " cyr("; I - 1; ") = "; cyr(I) Form1.Print " cyi("; I - 1; ") = "; cyi(I) Next I Form1.Print " NZ = "; NZ Form1.Print " Error code: "; IERR Form1.Print End Sub 'end of file tzbesy.bas