📄 cbess00_cpp.txt
字号:
/*************************************************************************
* 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. *
* *
* C++ Release By J-P Moreau, Paris (07/01/2005). *
*************************************************************************/
#include <math.h>
#include <basis.h>
#include <vmblock.h>
#include "complex.h"
//Functions defined below.
void ZACAI(REAL, REAL, REAL, int, int, int, REAL *, REAL *, int *, REAL,
REAL, REAL, REAL);
void ZBKNU(REAL, REAL, REAL, int, int, REAL *, REAL *, int *, REAL,
REAL, REAL);
void ZKSCL(REAL, REAL, REAL, int, REAL *, REAL *, int *, REAL *, REAL *,
REAL, REAL, REAL);
//Functions defined in CBess0.cpp.
void ZSERI(REAL, REAL, REAL, int, int, REAL *, REAL *, int *, REAL, REAL, REAL);
void ZASYI(REAL, REAL, REAL, int, int, REAL *, REAL *, int *, REAL, REAL, REAL, REAL);
void ZMLRI(REAL, REAL, REAL, int, int, REAL *, REAL *, int *, REAL);
void ZS1S2(REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, int *, REAL, REAL, int *);
void ZSHCH(REAL, REAL, REAL *, REAL *, REAL *, REAL *);
REAL DGAMLN(REAL, int *);
void ZUCHK(REAL, REAL, int *, REAL, REAL);
void ZAIRY(REAL ZR, REAL ZI, int ID, int KODE, REAL *AIR, REAL *AII,
int *NZ, int *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)
! C=1.0/(PI*SQRT(3.0))
! 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=SQRT(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: e40,e50,e60,e70,e80,e90,e100,e110,e120,e130,e140,e150,e160,e170,
// e180,e190,e200,e210,e260,e270,e280
REAL 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;
int IFLAG, K, K1, K2, MR, NN;
REAL CYI[10], CYR[10];
TTH = 6.66666666666666667e-01;
C1 = 3.55028053887817240e-01;
C2 = 2.58819403792806799e-01;
COEF= 1.83776298473930683e-01;
ZEROR=0.0; ZEROI=0.0; CONER=1.0; CONEI=0.0;
*IERR = 0;
*NZ=0;
if (ID < 0 || ID >1) *IERR=1;
if (KODE < 1 || KODE > 2) *IERR=1;
if (*IERR != 0) return;
AZ = ZABS(ZR,ZI);
TOL = DMAX(D1MACH(4),1e-18);
FID = 1.0*ID;
if (AZ > 1.0) goto e70;
/*----------------------------------------------------------------------
! POWER SERIES FOR CABS(Z).LE.1.
!---------------------------------------------------------------------*/
S1R = CONER;
S1I = CONEI;
S2R = CONER;
S2I = CONEI;
if (AZ < TOL) goto e170;
AA = AZ*AZ;
if (AA < TOL/AZ) goto e40;
TRM1R = CONER;
TRM1I = CONEI;
TRM2R = CONER;
TRM2I = CONEI;
ATRM = 1.0;
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.0 + FID;
BK = 3.0 - FID - FID;
CK = 4.0 - FID;
DK = 3.0 + FID + FID;
D1 = AK*DK;
D2 = BK*CK;
AD = DMIN(D1,D2);
AK = 24.0 + 9.0*FID;
BK = 30.0 - 9.0*FID;
for (K=1; K<=25; K++) {
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) goto e40;
AK += 18.0;
BK += 18.0;
}
e40: if (ID == 1) goto e50;
*AIR = S1R*C1 - C2*(ZR*S2R-ZI*S2I);
*AII = S1I*C1 - C2*(ZR*S2I+ZI*S2R);
if (KODE == 1) return;
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;
return;
e50: *AIR = -S2R*C2;
*AII = -S2I*C2;
if (AZ <= TOL) goto e60;
STR = ZR*S1R - ZI*S1I;
STI = ZR*S1I + ZI*S1R;
CC = C1/(1.0+FID);
*AIR = *AIR + CC*(STR*ZR-STI*ZI);
*AII = *AII + CC*(STR*ZI+STI*ZR);
e60: if (KODE == 1) return;
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;
return;
/*----------------------------------------------------------------------
! CASE FOR CABS(Z).GT.1.0
!---------------------------------------------------------------------*/
e70: FNU = (1.0+FID)/3.0;
/*----------------------------------------------------------------------
! 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 = I1MACH(15);
K2 = I1MACH(16);
R1M5 = D1MACH(5);
K = IMIN(ABS(K1),ABS(K2));
ELIM = 2.303*(K*R1M5-3.0);
K1 = I1MACH(14) - 1;
AA = R1M5*K1;
DIG = DMIN(AA,18.0);
AA = AA*2.303;
ALIM = ELIM + DMAX(-AA,-41.45);
RL = 1.2*DIG + 3.0;
ALAZ = log(AZ);
/*----------------------------------------------------------------------
! TEST FOR PROPER RANGE
!---------------------------------------------------------------------*/
AA=0.5/TOL;
BB=0.5*I1MACH(9);
AA=DMIN(AA,BB);
AA=pow(AA,TTH);
if (AZ > AA) goto e260;
AA=SQRT(AA);
if (AZ > AA) *IERR=3;
ZSQRT(ZR, ZI, &CSQR, &CSQI);
ZTAR = TTH*(ZR*CSQR-ZI*CSQI);
ZTAI = TTH*(ZR*CSQI+ZI*CSQR);
/*----------------------------------------------------------------------
! RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL
!---------------------------------------------------------------------*/
IFLAG = 0;
SFAC = 1.0;
AK = ZTAI;
if (ZR >= 0.0) goto e80;
BK = ZTAR;
CK = -ABS(BK);
ZTAR = CK;
ZTAI = AK;
e80: if (ZI != 0.0) goto e90;
if (ZR > 0.0) goto e90;
ZTAR = 0.0;
ZTAI = AK;
e90: AA = ZTAR;
if (AA >= 0.0 && ZR > 0.0) goto e110;
if (KODE == 2) goto e100;
/*----------------------------------------------------------------------
! OVERFLOW TEST
!---------------------------------------------------------------------*/
if (AA > -ALIM) goto e100;
AA = -AA + 0.25*ALAZ;
IFLAG = 1;
SFAC = TOL;
if (AA > ELIM) goto e270;
/*----------------------------------------------------------------------
! CBKNU AND CACON RETURN EXP(ZTA)*K(FNU,ZTA) ON KODE=2
!---------------------------------------------------------------------*/
e100: MR = 1;
if (ZI < 0.0) MR = -1;
ZACAI(ZTAR, ZTAI, FNU, KODE, MR, 1, CYR, CYI, &NN, RL, TOL, ELIM, ALIM);
if (NN < 0) goto e280;
*NZ = *NZ + NN;
goto e130;
e110: if (KODE == 2) goto e120;
/*----------------------------------------------------------------------
! UNDERFLOW TEST
!---------------------------------------------------------------------*/
if (AA < ALIM) goto e120;
AA = -AA - 0.25*ALAZ;
IFLAG = 2;
SFAC = 1.0/TOL;
if (AA < -ELIM) goto e210;
e120: ZBKNU(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, NZ, TOL, ELIM, ALIM);
e130: S1R = CYR[1]*COEF;
S1I = CYI[1]*COEF;
if (IFLAG != 0) goto e150;
if (ID == 1) goto e140;
*AIR = CSQR*S1R - CSQI*S1I;
*AII = CSQR*S1I + CSQI*S1R;
return;
e140: *AIR = -(ZR*S1R-ZI*S1I);
*AII = -(ZR*S1I+ZI*S1R);
return;
e150: S1R = S1R*SFAC;
S1I = S1I*SFAC;
if (ID == 1) goto e160;
STR = S1R*CSQR - S1I*CSQI;
S1I = S1R*CSQI + S1I*CSQR;
S1R = STR;
*AIR = S1R/SFAC;
*AII = S1I/SFAC;
return;
e160: STR = -(S1R*ZR-S1I*ZI);
S1I = -(S1R*ZI+S1I*ZR);
S1R = STR;
*AIR = S1R/SFAC;
*AII = S1I/SFAC;
return;
e170: AA = 1000*D1MACH(1);
S1R = ZEROR;
S1I = ZEROI;
if (ID == 1) goto e190;
if (AZ <= AA) goto e180;
S1R = C2*ZR;
S1I = C2*ZI;
e180: *AIR = C1 - S1R;
*AII = -S1I;
return;
e190: *AIR = -C2;
*AII = 0.0;
AA = SQRT(AA);
if (AZ <= AA) goto e200;
S1R = 0.5*(ZR*ZR-ZI*ZI);
S1I = ZR*ZI;
e200: *AIR = *AIR + C1*S1R;
*AII = *AII + C1*S1I;
return;
e210: *NZ = 1;
*AIR = ZEROR;
*AII = ZEROI;
return;
e270: *NZ = 0;
*IERR=2;
return;
e280: if (NN == -1) goto e270;
*NZ=0;
*IERR=5;
return;
e260: *IERR=4;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -