📄 cbess3_cpp.txt
字号:
/*************************************************************************
* 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. *
* *
* C++ Release By J-P Moreau, Paris (07/10/2005). *
*************************************************************************/
#include <basis.h>
#include <vmblock.h>
#include "complex.h"
//Headers of functions used below
void ZACON(REAL ZR, REAL ZI, REAL FNU, int KODE, int MR, int N, REAL *YR,
REAL *YI, int *NZ, REAL RL, REAL FNUL, REAL TOL, REAL ELIM, REAL ALIM);
void ZUNK1(REAL ZR, REAL ZI, REAL FNU, int KODE, int MR, int N, REAL *YR,
REAL *YI, int *NZ, REAL TOL, REAL ELIM, REAL ALIM);
void ZUNK2(REAL ZR, REAL ZI, REAL FNU, int KODE, int MR, int N, REAL *YR,
REAL *YI, int *NZ, REAL TOL, REAL ELIM, REAL ALIM);
void ZUNIK(REAL, REAL, REAL, int, int, REAL, int, REAL *, REAL *, REAL *, REAL *,
REAL *, REAL *, REAL *, REAL *, REAL *, REAL *);
void ZUCHK(REAL, REAL, int *, REAL, REAL);
void ZS1S2(REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, int *, REAL, REAL, int *);
void ZUNHJ(REAL, REAL, REAL, int, REAL, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *,
REAL *, REAL *, REAL *, REAL *, REAL *, REAL *);
void ZAIRY(REAL, REAL, int, int, REAL *, REAL *, int *, int *);
void ZBINU(REAL, REAL, REAL, int, int, REAL *, REAL *, int *, REAL, REAL, REAL,
REAL, REAL);
void ZBKNU(REAL, REAL, REAL, int, int, REAL *, REAL *, int *, REAL,
REAL, REAL);
void ZBUNK(REAL ZR, REAL ZI, REAL FNU, int KODE, int MR, int N, REAL *YR, REAL *YI,
int *NZ, REAL TOL, REAL ELIM, REAL 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 */
//Label: e10
REAL AX, AY;
*NZ = 0;
AX = ABS(ZR)*1.7321;
AY = ABS(ZI);
if (AY > AX) goto e10;
/*----------------------------------------------------------------------
! 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);
return;
/*----------------------------------------------------------------------
! ASYMPTOTIC EXPANSION FOR H(2,FNU,Z*EXP(M*HPI)) FOR LARGE FNU
! APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
! AND HPI=PI/2
!---------------------------------------------------------------------*/
e10: ZUNK2(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM, ALIM);
} //ZBUNK()
void ZUNK1(REAL ZR, REAL ZI, REAL FNU, int KODE, int MR, int N, REAL *YR,
REAL *YI, int *NZ, REAL TOL, REAL ELIM, REAL 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: e10,e20,e30,e40,e50,e60,e70,e75,e80,e90,e95,e100,e120,e160,e170,
// e172,e175,e180,e200,e210,e220,e230,e250,e255,e260,e270,e275,e280,
// e290,e300
REAL 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;
int I, IB, IFLAG, IFN, II, IL, INU, IUF, K, KDFLG, KFLAG, KK, NW, INITD,
IC, IPARD, J, M;
REAL *BRY, *SUMR, *SUMI, *ZETA1R, *ZETA1I, *ZETA2R, *ZETA2I, *CYR, *CYI;
int *INIT; //size 0..3
REAL **CWRKR, **CWRKI; //size 0..16, 0..3
REAL *CSSR, *CSRR, *PHIR, *PHII;
REAL *TMPR, *TMPI; //size 0..16
void *vmblock = NULL;
//*** First executable statement ZUNK1
//initialize pointers to vectors
vmblock = vminit();
BRY = (REAL *) vmalloc(vmblock, VEKTOR, 4, 0); //index 0 not used
SUMR = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
SUMI = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
ZETA1R = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
ZETA1I = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
ZETA2R = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
ZETA2I = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
CYR = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
CYI = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
INIT = (int *) vmalloc(vmblock, VEKTOR, 4, 0);
CWRKR = (REAL **) vmalloc(vmblock, MATRIX, 17, 4);
CWRKI = (REAL **) vmalloc(vmblock, MATRIX, 17, 4);
TMPR = (REAL *) vmalloc(vmblock, VEKTOR, 17, 0);
TMPI = (REAL *) vmalloc(vmblock, VEKTOR, 17, 0);
CSSR = (REAL *) vmalloc(vmblock, VEKTOR, 4, 0);
CSRR = (REAL *) vmalloc(vmblock, VEKTOR, 4, 0);
PHIR = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
PHII = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
if (! vmcomplete(vmblock)) {
LogError ("No Memory", 0, __FILE__, __LINE__);
return;
}
printf(" ZUNK1:allocations Ok.\n"); getchar();
ZEROR=0.0; ZEROI=0.0; CONER=1.0; CONEI=0.0;
KDFLG = 1;
*NZ = 0;
/*----------------------------------------------------------------------
! EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
! THE UNDERFLOW LIMIT
!---------------------------------------------------------------------*/
CSCL = 1.0/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.0/BRY[1];
BRY[3] = D1MACH(2);
ZRR = ZR;
ZRI = ZI;
if (ZR >= 0.0) goto e10;
ZRR = -ZR;
ZRI = -ZI;
e10: J = 2;
for (I=1; I<=N; I++) {
/*----------------------------------------------------------------------
! J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
!---------------------------------------------------------------------*/
J = 3 - J;
FN = FNU + 1.0*(I-1);
INIT[J] = 0;
for (II=1; II<=16; II++) {
TMPR[II]=CWRKR[II][J];
TMPI[II]=CWRKI[II][J];
}
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; II<=16; II++) {
CWRKR[II][J]=TMPR[II];
CWRKI[II][J]=TMPI[II];
}
if (KODE == 1) goto e20;
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 e30;
e20: S1R = ZETA1R[J] - ZETA2R[J];
S1I = ZETA1I[J] - ZETA2I[J];
e30: RS1 = S1R;
/*----------------------------------------------------------------------
! TEST FOR UNDERFLOW AND OVERFLOW
!---------------------------------------------------------------------*/
if (ABS(RS1) > ELIM) goto e60;
if (KDFLG == 1) KFLAG = 2;
if (ABS(RS1) < ALIM) goto e40;
/*----------------------------------------------------------------------
! REFINE TEST AND SCALE
!---------------------------------------------------------------------*/
APHI = ZABS(PHIR[J],PHII[J]);
RS1 = RS1 + log(APHI);
if (ABS(RS1) > ELIM) goto e60;
if (KDFLG == 1) KFLAG = 1;
if (RS1 < 0.0) goto e40;
if (KDFLG == 1) KFLAG = 3;
/*----------------------------------------------------------------------
! SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
! EXPONENT EXTREMES
!---------------------------------------------------------------------*/
e40: 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) goto e50;
ZUCHK(S2R, S2I, &NW, BRY[1], TOL);
if (NW != 0) goto e60;
e50: CYR[KDFLG] = S2R;
CYI[KDFLG] = S2I;
YR[I] = S2R*CSRR[KFLAG];
YI[I] = S2I*CSRR[KFLAG];
if (KDFLG == 2) goto e75;
KDFLG = 2;
goto e70;
e60: if (RS1 > 0.0) goto e300;
/*----------------------------------------------------------------------
! FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
!---------------------------------------------------------------------*/
if (ZR < 0.0) goto e300;
KDFLG = 1;
YR[I]=ZEROR;
YI[I]=ZEROI;
NZ=NZ+1;
if (I == 1) goto e70;
if (YR[I-1] == ZEROR && YI[I-1] == ZEROI) goto e70;
YR[I-1]=ZEROR;
YI[I-1]=ZEROI;
*NZ=*NZ+1;
e70: ;} // I loop
I = N;
e75: RAZR = 1.0/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) goto e160;
/*----------------------------------------------------------------------
! TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
! ON UNDERFLOW.
!---------------------------------------------------------------------*/
FN = FNU + 1.0*(N-1);
IPARD = 1;
if (MR != 0) IPARD = 0;
INITD = 0;
for (II=1; II<=16; II++) {
TMPR[II]=CWRKR[II][3];
TMPI[II]=CWRKI[II][3];
}
ZUNIK(ZRR, ZRI, FN, 2, IPARD, TOL, INITD, &PHIDR, &PHIDI, &ZET1DR,
&ZET1DI, &ZET2DR, &ZET2DI, &SUMDR, &SUMDI, TMPR, TMPI);
for (II=1; II<=16; II++) {
CWRKR[II][3]=TMPR[II];
CWRKI[II][3]=TMPI[II];
}
if (KODE == 1) goto e80;
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 e90;
e80: S1R = ZET1DR - ZET2DR;
S1I = ZET1DI - ZET2DI;
e90: RS1 = S1R;
if (ABS(RS1) > ELIM) goto e95;
if (ABS(RS1) < ALIM) goto e100;
/*----------------------------------------------------------------------
! REFINE ESTIMATE AND TEST
!---------------------------------------------------------------------*/
APHI = ZABS(PHIDR,PHIDI);
RS1 = RS1+log(APHI);
if (ABS(RS1) < ELIM) goto e100;
e95: if (ABS(RS1) > 0.0) goto e300;
/*----------------------------------------------------------------------
! FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
!---------------------------------------------------------------------*/
if (ZR < 0.0) goto e300;
*NZ = N;
for (I=1; I<=N; I++) {
YR[I] = ZEROR;
YI[I] = ZEROI;
}
vmfree(vmblock);
return;
/*----------------------------------------------------------------------
! FORWARD RECUR FOR REMAINDER OF THE SEQUENCE
!---------------------------------------------------------------------*/
e100: S1R = CYR[1];
S1I = CYI[1];
S2R = CYR[2];
S2I = CYI[2];
C1R = CSRR[KFLAG];
ASCLE = BRY[KFLAG];
for (I=IB; I<=N; I++) {
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) goto e120;
STR = ABS(C2R);
STI = ABS(C2I);
C2M = DMAX(STR,STI);
if (C2M <= ASCLE) goto e120;
KFLAG++;
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];
e120:;} // I loop
e160: if (MR == 0) {
vmfree(vmblock);
return;
}
/*----------------------------------------------------------------------
! ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
!---------------------------------------------------------------------*/
*NZ = 0;
FMR = 1.0*MR;
SGN = -SIGN(PI,FMR);
/*----------------------------------------------------------------------
! CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
!---------------------------------------------------------------------*/
CSGNI = SGN;
INU = (int) floor(FNU);
FNF = FNU - 1.0*INU;
IFN = INU + N - 1;
ANG = FNF*SGN;
CSPNR = COS(ANG);
CSPNI = SIN(ANG);
if ((IFN % 2) == 0) goto e170;
CSPNR = -CSPNR;
CSPNI = -CSPNI;
e170: ASC = BRY[1];
IUF = 0;
KK = N;
KDFLG = 1;
IB = IB - 1;
IC = IB - 1;
for (K=1; K<=N; K++) {
FN = FNU + 1.0*(KK-1);
/*----------------------------------------------------------------------
! LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
! FUNCTION ABOVE
!---------------------------------------------------------------------*/
M=3;
if (N > 2) goto e175;
e172: 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 e180;
e175: if (KK == N && IB < N) goto e180;
if (KK == IB || KK == IC) goto e172;
INITD = 0;
e180: for (II=1; II<=16; II++) {
TMPR[II]=CWRKR[II][M];
TMPI[II]=CWRKI[II][M];
}
ZUNIK(ZRR, ZRI, FN, 1, 0, TOL, INITD, &PHIDR, &PHIDI,
&ZET1DR, &ZET1DI, &ZET2DR, &ZET2DI, &SUMDR, &SUMDI,
TMPR,TMPI);
for (II=1; II<=16; II++) {
CWRKR[II][M]=TMPR[II];
CWRKI[II][M]=TMPI[II];
}
if (KODE == 1) goto e200;
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 e210;
e200: S1R = -ZET1DR + ZET2DR;
S1I = -ZET1DI + ZET2DI;
/*----------------------------------------------------------------------
! TEST FOR UNDERFLOW AND OVERFLOW
!---------------------------------------------------------------------*/
e210: RS1 = S1R;
if (ABS(RS1) > ELIM) goto e260;
if (KDFLG == 1) IFLAG = 2;
if (ABS(RS1) < ALIM) goto e220;
/*----------------------------------------------------------------------
! REFINE TEST AND SCALE
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -