📄 cbess00_cpp.txt
字号:
*NZ=0;
} //ZAIRY()
void ZACAI(REAL ZR, REAL ZI, REAL FNU, int KODE, int MR, int N,
REAL *YR, REAL *YI, int *NZ, REAL RL, REAL TOL, REAL ELIM, REAL ALIM) {
/***BEGIN PROLOGUE ZACAI
!***REFER TO ZAIRY
!
! ZACAI APPLIES THE ANALYTIC CONTINUATION FORMULA
!
! K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
! MP=PI*MR*CMPLX(0.0,1.0)
!
! TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
! HALF Z PLANE FOR USE WITH ZAIRY WHERE FNU=1/3 OR 2/3 AND N=1.
! ZACAI IS THE SAME AS ZACON WITH THE PARTS FOR LARGER ORDERS AND
! RECURRENCE REMOVED. A RECURSIVE CALL TO ZACON CAN RESULT IF ZACON
! IS CALLED FROM ZAIRY.
!
!***ROUTINES CALLED ZASYI,ZBKNU,ZMLRI,ZSERI,ZS1S2,D1MACH,ZABS
!***END PROLOGUE ZACAI
! COMPLEX CSGN,CSPN,C1,C2,Y,Z,ZN,CY */
//Labels: e10, e20, e30, e40, e50, e60, e70, e80
REAL ARG, ASCLE, AZ, CSGNR, CSGNI, CSPNR, CSPNI, C1R, C1I, C2R, C2I,
DFNU, FMR, SGN, YY, ZNR, ZNI;
int INU, IUF, NN, NW;
REAL *CYR, *CYI;
void *vmblock = NULL;
// Initialize CYR, CYI
vmblock = vminit();
CYR = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0); //index 0 not used
CYI = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
if (! vmcomplete(vmblock)) {
LogError ("No Memory", 0, __FILE__, __LINE__);
return;
}
*NZ = 0;
ZNR = -ZR;
ZNI = -ZI;
AZ = ZABS(ZR,ZI);
NN = N;
DFNU = FNU + 1.0*(N-1);
if (AZ <= 2.0) goto e10;
if (AZ*AZ*0.25 > DFNU+1.0) goto e20;
/*----------------------------------------------------------------------
! POWER SERIES FOR THE I FUNCTION
!---------------------------------------------------------------------*/
e10: ZSERI(ZNR, ZNI, FNU, KODE, NN, YR, YI, &NW, TOL, ELIM, ALIM);
goto e40;
e20: if (AZ < RL) goto e30;
/*----------------------------------------------------------------------
! ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION
!---------------------------------------------------------------------*/
ZASYI(ZNR, ZNI, FNU, KODE, NN, YR, YI, &NW, RL, TOL, ELIM, ALIM);
if (NW < 0) goto e80;
goto e40;
/*----------------------------------------------------------------------
! MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION
!---------------------------------------------------------------------*/
e30: ZMLRI(ZNR, ZNI, FNU, KODE, NN, YR, YI, &NW, TOL);
if (NW < 0) goto e80;
/*----------------------------------------------------------------------
! ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
!---------------------------------------------------------------------*/
e40: ZBKNU(ZNR, ZNI, FNU, KODE, 1, CYR, CYI, &NW, TOL, ELIM, ALIM);
if (NW != 0) goto e80;
FMR = 1.0*MR;
SGN = -SIGN(PI,FMR);
CSGNR = 0.0;
CSGNI = SGN;
if (KODE == 1) goto e50;
YY = -ZNI;
CSGNR = -CSGNI*SIN(YY);
CSGNI = CSGNI*COS(YY);
/*----------------------------------------------------------------------
! CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
! WHEN FNU IS LARGE
!---------------------------------------------------------------------*/
e50: INU = (int) floor(FNU);
ARG = (FNU-1.0*INU)*SGN;
CSPNR = COS(ARG);
CSPNI = SIN(ARG);
if ((INU % 2) == 0) goto e60;
CSPNR = -CSPNR;
CSPNI = -CSPNI;
e60: C1R = CYR[1];
C1I = CYI[1];
C2R = YR[1];
C2I = YI[1];
if (KODE == 1) goto e70;
IUF = 0;
ASCLE = 1000*D1MACH(1)/TOL;
ZS1S2(&ZNR, &ZNI, &C1R, &C1I, &C2R, &C2I, &NW, ASCLE, ALIM, &IUF);
*NZ = *NZ + NW;
e70: YR[1] = CSPNR*C1R - CSPNI*C1I + CSGNR*C2R - CSGNI*C2I;
YI[1] = CSPNR*C1I + CSPNI*C1R + CSGNR*C2I + CSGNI*C2R;
vmfree(vmblock);
return;
e80: *NZ = -1;
if (NW == -2) *NZ=-2;
vmfree(vmblock);
} //ZACAI()
void ZBKNU(REAL ZR, REAL ZI, REAL FNU, int KODE, int N, REAL *YR, REAL *YI,
int *NZ, REAL TOL, REAL ELIM, REAL ALIM) {
/***BEGIN PROLOGUE ZBKNU
!***REFER TO ZBESI,ZBESK,ZAIRY,ZBESH
!
! ZBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE.
!
!***ROUTINES CALLED DGAMLN,I1MACH,D1MACH,ZKSCL,ZSHCH,ZUCHK,ZABS,ZDIV,
! ZEXP,ZLOG,ZMLT,ZSQRT
!***END PROLOGUE ZBKNU */
//Labels: e10,e30,e40,e50,e60,e70,e80,e90,e100,e110,e120,e130,e140,e160,e170,
// e180,e200,e210,e215,e220,e225,e230,e240,e250,e260,e261,e262,e263,
// e264,e270,e280,e290,e300,e310
REAL AA, AK, ASCLE, A1, A2, BB, BK, CAZ,
CBI, CBR, CCHI, CCHR, CKI, CKR, COEFI, COEFR, CONEI, CONER,
CRSCR, CSCLR, CSHI, CSHR, CSI, CSR, CTWOI, CTWOR,
CZEROI, CZEROR, CZI, CZR, DNU, DNU2, DPI, ETEST, FC, FHS,
FI, FK, FKS, FMUI, FMUR, FPI, FR, G1, G2, HPI, PI0, PR, PTI,
PTR, P1I, P1R, P2I, P2M, P2R, QI, QR, RAK, RCAZ, RTHPI, RZI,
RZR, R1, S, SMUI, SMUR, SPI, STI, STR, S1I, S1R, S2I, S2R, TM,
TTH, T1, T2, ELM, CELMR, ZDR, ZDI, AS, ALAS, HELIM;
int I, IFLAG, INU, K, KFLAG, KK, KMAX, KODED, IDUM, J, IC, INUB, NW;
REAL *CC, *CSSR, *CSRR, *BRY, *CYR, *CYI;
void *vmblock = NULL;
// COMPLEX: Z,Y,A,B,RZ,SMU,FU,FMU,F,FLRZ,CZ,S1,S2,CSH,CCH,
// CK,P,Q,COEF,P1,P2,CBK,PT,CZERO,CONE,CTWO,ST,EZ,CS,DK
// Note: PI is replaced here by PI0 (PI is defined in basis.cpp as PI=3.1415926...).
KMAX=30;
CZEROR=0.0; CZEROI=0.0; CONER=1.0; CONEI=0.0;
CTWOR=2.0; CTWOI=0.0; R1=2.0;
DPI=3.14159265358979324; RTHPI=1.25331413731550025;
SPI=1.90985931710274403; HPI=1.57079632679489662;
FPI=1.89769999331517738; TTH=6.66666666666666666e-01;
// Initialize CC, ..., CYI
vmblock = vminit();
CC = (REAL *) vmalloc(vmblock, VEKTOR, 9, 0); //index 0 not used
CSSR = (REAL *) vmalloc(vmblock, VEKTOR, 4, 0);
CSRR = (REAL *) vmalloc(vmblock, VEKTOR, 4, 0);
BRY = (REAL *) vmalloc(vmblock, VEKTOR, 4, 0);
CYR = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
CYI = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
if (! vmcomplete(vmblock)) {
LogError ("No Memory", 0, __FILE__, __LINE__);
return;
}
CC[1]= 5.77215664901532861e-01; CC[2]=-4.20026350340952355e-02;
CC[3]=-4.21977345555443367e-02; CC[4]= 7.21894324666309954e-03;
CC[5]=-2.15241674114950973e-04; CC[6]=-2.01348547807882387e-05;
CC[7]= 1.13302723198169588e-06; CC[8]= 6.11609510448141582e-09;
CAZ = ZABS(ZR,ZI);
CSCLR = 1.0/TOL;
CRSCR = TOL;
CSSR[1] = CSCLR;
CSSR[2] = 1.0;
CSSR[3] = CRSCR;
CSRR[1] = CRSCR;
CSRR[2] = 1.0;
CSRR[3] = CSCLR;
BRY[1] = 1000*D1MACH(1)/TOL;
BRY[2] = 1.0/BRY[1];
BRY[3] = D1MACH(2);
*NZ = 0;
IFLAG = 0;
KODED = KODE;
RCAZ = 1.0/CAZ;
STR = ZR*RCAZ;
STI = -ZI*RCAZ;
RZR = (STR+STR)*RCAZ;
RZI = (STI+STI)*RCAZ;
INU = (int) floor(FNU+0.5);
DNU = FNU - 1.0*INU;
if (ABS(DNU) == 0.5) goto e110;
DNU2 = 0.0;
if (ABS(DNU) > TOL) DNU2 = DNU*DNU;
if (CAZ > R1) goto e110;
/*----------------------------------------------------------------------
! SERIES FOR CABS(Z).LE.R1
!---------------------------------------------------------------------*/
FC = 1.0;
ZLOG(RZR, RZI, &SMUR, &SMUI, &IDUM);
FMUR = SMUR*DNU;
FMUI = SMUI*DNU;
ZSHCH(FMUR, FMUI, &CSHR, &CSHI, &CCHR, &CCHI);
if (DNU == 0.0) goto e10;
FC = DNU*DPI;
FC = FC/SIN(FC);
SMUR = CSHR/DNU;
SMUI = CSHI/DNU;
e10: A2 = 1.0 + DNU;
/*----------------------------------------------------------------------
! GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU)
!---------------------------------------------------------------------*/
T2 = EXP(-DGAMLN(A2,&IDUM));
T1 = 1.0/(T2*FC);
if (ABS(DNU) > 0.1) goto e40;
/*----------------------------------------------------------------------
! SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
!---------------------------------------------------------------------*/
AK = 1.0;
S = CC[1];
for (K=2; K<9; K++) {
AK = AK*DNU2;
TM = CC[K]*AK;
S = S + TM;
if (ABS(TM) < TOL) goto e30;
}
e30: G1 = -S;
goto e50;
e40: G1 = (T1-T2)/(DNU+DNU);
e50: G2 = (T1+T2)*0.5;
FR = FC*(CCHR*G1+SMUR*G2);
FI = FC*(CCHI*G1+SMUI*G2);
ZEXP(FMUR, FMUI, &STR, &STI);
PR = 0.5*STR/T2;
PI0 = 0.5*STI/T2;
ZDIV(0.5, 0.0, STR, STI, &PTR, &PTI);
QR = PTR/T1;
QI = PTI/T1;
S1R = FR;
S1I = FI;
S2R = PR;
S2I = PI0;
AK = 1.0;
A1 = 1.0;
CKR = CONER;
CKI = CONEI;
BK = 1.0 - DNU2;
if (INU > 0 || N >1) goto e80;
/*----------------------------------------------------------------------
! GENERATE K(FNU,Z), 0.0D0 .LE. FNU .LT. 0.5D0 AND N=1
!---------------------------------------------------------------------*/
if (CAZ < TOL) goto e70;
ZMLT(ZR, ZI, ZR, ZI, &CZR, &CZI);
CZR = 0.25*CZR;
CZI = 0.25*CZI;
T1 = 0.25*CAZ*CAZ;
e60: FR = (FR*AK+PR+QR)/BK;
FI = (FI*AK+PI0+QI)/BK;
STR = 1.0/(AK-DNU);
PR = PR*STR;
PI0 = PI0*STR;
STR = 1.0/(AK+DNU);
QR = QR*STR;
QI = QI*STR;
STR = CKR*CZR - CKI*CZI;
RAK = 1.0/AK;
CKI = (CKR*CZI+CKI*CZR)*RAK;
CKR = STR*RAK;
S1R = CKR*FR - CKI*FI + S1R;
S1I = CKR*FI + CKI*FR + S1I;
A1 = A1*T1*RAK;
BK = BK + AK + AK + 1.0;
AK = AK + 1.0;
if (A1 > TOL) goto e60;
e70: YR[1] = S1R;
YI[1] = S1I;
if (KODED == 1) {
vmfree(vmblock);
return;
}
ZEXP(ZR, ZI, &STR, &STI);
ZMLT(S1R, S1I, STR, STI, &YR[1], &YI[1]);
vmfree(vmblock);
return;
/*----------------------------------------------------------------------
! GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE
!---------------------------------------------------------------------*/
e80: if (CAZ < TOL) goto e100;
ZMLT(ZR, ZI, ZR, ZI, &CZR, &CZI);
CZR = 0.25*CZR;
CZI = 0.25*CZI;
T1 = 0.25*CAZ*CAZ;
e90: FR = (FR*AK+PR+QR)/BK;
FI = (FI*AK+PI0+QI)/BK;
STR = 1.0/(AK-DNU);
PR = PR*STR;
PI0 = PI0*STR;
STR = 1.0/(AK+DNU);
QR = QR*STR;
QI = QI*STR;
STR = CKR*CZR - CKI*CZI;
RAK = 1.0/AK;
CKI = (CKR*CZI+CKI*CZR)*RAK;
CKR = STR*RAK;
S1R = CKR*FR - CKI*FI + S1R;
S1I = CKR*FI + CKI*FR + S1I;
STR = PR - FR*AK;
STI = PI0 - FI*AK;
S2R = CKR*STR - CKI*STI + S2R;
S2I = CKR*STI + CKI*STR + S2I;
A1 = A1*T1*RAK;
BK = BK + AK + AK + 1.0;
AK = AK + 1.0;
if (A1 > TOL) goto e90;
e100: KFLAG = 2;
A1 = FNU + 1.0;
AK = A1*ABS(SMUR);
if (AK > ALIM) KFLAG = 3;
STR = CSSR[KFLAG];
P2R = S2R*STR;
P2I = S2I*STR;
ZMLT(P2R, P2I, RZR, RZI, &S2R, &S2I);
S1R = S1R*STR;
S1I = S1I*STR;
if (KODED == 1) goto e210;
ZEXP(ZR, ZI, &FR, &FI);
ZMLT(S1R, S1I, FR, FI, &S1R, &S1I);
ZMLT(S2R, S2I, FR, FI, &S2R, &S2I);
goto e210;
/*----------------------------------------------------------------------
! IFLAG=0 MEANS NO UNDERFLOW OCCURRED
! IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH
! KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD
! RECURSION
!---------------------------------------------------------------------*/
e110: ZSQRT(ZR, ZI, &STR, &STI);
ZDIV(RTHPI, CZEROI, STR, STI, &COEFR, &COEFI);
KFLAG = 2;
if (KODED == 2) goto e120;
if (ZR > ALIM) goto e290;
STR = EXP(-ZR)*CSSR[KFLAG];
STI = -STR*SIN(ZI);
STR = STR*COS(ZI);
ZMLT(COEFR, COEFI, STR, STI, &COEFR, &COEFI);
e120: if (ABS(DNU) == 0.5) goto e300;
/*----------------------------------------------------------------------
! MILLER ALGORITHM FOR CABS(Z).GT.R1
!---------------------------------------------------------------------*/
AK = COS(DPI*DNU);
AK = ABS(AK);
if (AK == CZEROR) goto e300;
FHS = ABS(0.25-DNU2);
if (FHS == CZEROR) goto e300;
/*----------------------------------------------------------------------
! COMPUTE R2=F(E). IF CABS(Z).GE.R2, USE FORWARD RECURRENCE TO
! DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON
! 12.LE.E.LE.60. E IS COMPUTED FROM 2**(-E)=B**(1-I1MACH(14))=
! TOL WHERE B IS THE BASE OF THE ARITHMETIC.
!---------------------------------------------------------------------*/
T1 = 1.0*(I1MACH(14)-1);
T1 = T1*D1MACH(5)*3.321928094;
T1 = DMAX(T1,12.0);
T1 = DMIN(T1,60.0);
T2 = TTH*T1 - 6.0;
if (ZR != 0.0) goto e130;
T1 = HPI;
goto e140;
e130: T1 = atan(ZI/ZR);
T1 = ABS(T1);
e140: if (T2 > CAZ) goto e170;
/*----------------------------------------------------------------------
! FORWARD RECURRENCE LOOP WHEN CABS(Z).GE.R2
!---------------------------------------------------------------------*/
ETEST = AK/(DPI*CAZ*TOL);
FK = CONER;
if (ETEST < CONER) goto e180;
FKS = CTWOR;
CKR = CAZ + CAZ + CTWOR;
P1R = CZEROR;
P2R = CONER;
for (I=1; I<=KMAX; I++) {
AK = FHS/FKS;
CBR = CKR/(FK+CONER);
PTR = P2R;
P2R = CBR*P2R - P1R*AK;
P1R = PTR;
CKR = CKR + CTWOR;
FKS = FKS + FK + FK + CTWOR;
FHS = FHS + FK + FK;
FK = FK + CONER;
STR = ABS(P2R)*FK;
if (ETEST < STR) goto e160;
}
goto e310;
e160: FK = FK + SPI*T1*SQRT(T2/CAZ);
FHS = ABS(0.25-DNU2);
goto e180;
/*----------------------------------------------------------------------
! COMPUTE BACKWARD INDEX K FOR CABS(Z).LT.R2
!---------------------------------------------------------------------*/
e170: A2 = SQRT(CAZ);
AK = FPI*AK/(TOL*SQRT(A2));
AA = 3.0*T1/(1.0+CAZ);
BB = 14.7*T1/(28.0+CAZ);
AK = (log(AK)+CAZ*COS(AA)/(1.0+0.008*CAZ))/COS(BB);
FK = 0.12125*AK*AK/CAZ + 1.5;
/*----------------------------------------------------------------------
! BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM
!---------------------------------------------------------------------*/
e180: K = (int) floor(FK);
FK = 1.0*K;
FKS = FK*FK;
P1R = CZEROR;
P1I = CZEROI;
P2R = TOL;
P2I = CZEROI;
CSR = P2R;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -