📄 tzbesy_cpp.txt
字号:
/****************************************************************
* EVALUATE A Y-BESSEL FUNCTION OF COMPLEX ARGUMENT (SECOND KIND)*
* ------------------------------------------------------------- *
* SAMPLE RUN: *
* (Evaluate Y0 to Y4 for argument Z=(1.0,2.0) ). *
* *
* zr(0) = 1.367419 *
* zi(0) = 1.521507 *
* zr(1) = -1.089470 *
* zi(1) = 1.314951 *
* zr(2) = -0.751245 *
* zi(2) = -0.123950 *
* zr(3) = 0.290153 *
* zi(3) = -0.212119 *
* zr(4) = 0.590344 *
* zi(4) = -0.826960 *
* NZ = 0 *
* Error code: 0 *
* *
* ------------------------------------------------------------- *
* Ref.: From Numath Library By Tuan Dang Trong in Fortran 77. *
* *
* C++ Release 1.0 By J-P Moreau, Paris *
*****************************************************************
Note: to link with CBess0,CBess00,CBess1,CBess2,CBess3,Complex.
-------------------------------------------------------------- */
#include <basis.h>
#include <vmblock.h>
#include "complex.h"
//Uses Complex, CBess0, CBess00, CBess1, CBess2, CBess3,
Basis_r, Vmblock.
REAL zr, zi;
REAL *cyr, *cyi, *cwr, *cwi;
int i, ierr, n, nz;
void *vmblock = NULL;
//Headers of functions used below
void ZBESH(REAL ZR, REAL ZI, REAL FNU, int KODE, int M, int N,
REAL *CYR, REAL *CYI, int *NZ, int *IERR);
void ZUOIK(REAL, REAL, REAL, int, int, int, REAL *, REAL *,
int *, REAL, REAL, REAL);
void ZBKNU(REAL, REAL, REAL, int, int, REAL *, REAL *, int *, REAL,
REAL, REAL);
void ZACON(REAL, REAL, REAL, int, int, int, REAL *, REAL *, int *,
REAL, REAL, REAL, REAL, REAL);
void ZBUNK(REAL, REAL, REAL, int, int, int, REAL *, REAL *, int *,
REAL, REAL, REAL);
void ZBESY(REAL ZR, REAL ZI, REAL FNU, int KODE, int N, REAL *CYR, REAL *CYI,
int *NZ, REAL *CWRKR, REAL *CWRKI, int *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.LT.ARG(Z).LE.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.NE.CMPLX(0.0D0,0.0D0),
! -PI.LT.ARG(Z).LE.PI
! FNU - ORDER OF INITIAL Y FUNCTION, FNU.GE.0.0D0
! 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.GE.1
! CWRKR, - DOUBLE PRECISION WORK VECTORS OF DIMENSION AT
! CWRKI AT 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.GT.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.GT.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: e60, e70, e90, e170
REAL C1I, C1R, C2I, C2R, ELIM, EXI, EXR, EY, HCII, STI, STR, TAY;
int I, K, K1, K2, NZ1, NZ2;
//***FIRST EXECUTABLE STATEMENT ZBESY
*IERR = 0;
*NZ=0;
if (ZR == 0.0 && ZI == 0.0) *IERR=1;
if (FNU < 0.0) *IERR=1;
if (KODE < 1 || KODE > 2) *IERR=1;
if (N < 1) *IERR=1;
if (*IERR != 0) return; //bad parameter(s)
HCII = 0.5;
ZBESH(ZR, ZI, FNU, KODE, 1, N, CYR, CYI, &NZ1, IERR);
if (*IERR != 0 && *IERR != 3) goto e170;
ZBESH(ZR, ZI, FNU, KODE, 2, N, CWRKR, CWRKI, &NZ2, IERR);
if (*IERR != 0 && *IERR != 3) goto e170;
*NZ = IMIN(NZ1,NZ2);
if (KODE == 2) goto e60;
for (I=1; I<=N; I++) {
STR = CWRKR[I] - CYR[I];
STI = CWRKI[I] - CYI[I];
CYR[I] = -STI*HCII;
CYI[I] = STR*HCII;
}
return;
e60: K1 = I1MACH(15);
K2 = I1MACH(16);
K = IMIN(ABS(K1),ABS(K2));
/*----------------------------------------------------------------------
! ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT
!---------------------------------------------------------------------*/
ELIM = 2.303*(K*D1MACH(5)-3.0);
EXR = COS(ZR);
EXI = SIN(ZR);
EY = 0.0;
TAY = ABS(ZI+ZI);
if (TAY < ELIM) EY = EXP(-TAY);
if (ZI < 0.0) goto e90;
C1R = EXR*EY;
C1I = EXI*EY;
C2R = EXR;
C2I = -EXI;
e70: *NZ = 0;
for (I=1; I<=N; I++) {
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.0 && STI == 0.0 && EY == 0.0) *NZ = *NZ + 1;
}
return;
e90: C1R = EXR;
C1I = EXI;
C2R = EXR*EY;
C2I = -EXI*EY;
goto e70;
e170: *NZ = 0;
} //ZBESY()
void ZBESH(REAL ZR, REAL ZI, REAL FNU, int KODE, int M, int N,
REAL *CYR, REAL *CYI, int *NZ, int *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.NE.CMPLX(0.0,0.0) IN THE CUT PLANE -PI.LT.ARG(Z).LE.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.NE.CMPLX(0.0D0,0.0D0),
! -PT.LT.ARG(Z).LE.PI
! FNU - ORDER OF INITIAL H FUNCTION, FNU.GE.0.0D0
! 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.GE.1
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -