/************************************************************* * CALCULATE THE Kth ZERO OF THE DERIVATIVE OF BESSEL * * FUNCTION OF ORDER N, J(N,X) * * ---------------------------------------------------------- * * SAMPLE RUN: * * (Calculate the 10th zero of the derivative of J(1,X) ). * * * * X= 30.60192297 * * * * BESSJP(X)= -1.743397e-016 * * * * ---------------------------------------------------------- * * Reference: From Numath Library By Tuan Dang Trong in * * Fortran 77 [BIBLI 18]. * * * * C++ Release By J-P Moreau, Paris * * (www.jpmoreau.fr) * *************************************************************/ #include #include #define TRUE 1 double RES; int K, N; double Sign(double a, double b) { if (b<0) return (-fabs(a)); else return (fabs(a)); } double BESSJ (int N, double X); double BESSJ0(double X); double BESSJ1(double X); double BESSJP(int N, double X); double ZEROJP(int N, int K) { /*------------------------------------------------------------------- ! CALCULATE THE Kth ZERO OF THE DERIVATIVE OF BESSEL FUNCTION ! OF ORDER N, J(N,X) !-------------------------------------------------------------------- ! CALLING MODE: ! RES = ZEROJP(N,K) ! ! INPUTS: ! N ORDER OF BESSEL FUNCTION J (INTEGER >= 0) I*4 ! K RANK OF ZERO (INTEGER > 0) I*4 ! OUTPUT: ! ZEROJP R*8 ! REFERENCE: ! ABRAMOWITZ M. & STEGUN IRENE A. ! HANDBOOK OF MATHEMATICAL FUNCTIONS !-------------------------------------------------------------------*/ // Labels: e10, e20, e25 double B0,B1,B2,B3,B5,B7,T0,T1,T3,T5,T7,FN,FK; double C1,C2,C3,C4,F1,F2,P,DP,P0,P1,Q0,Q1,TOL; double PI, ZJP; int IER,IT,NEV,NITMX; bool IMPROV; PI = 4.0*atan(1.0); TOL=1e-7; NITMX=15; C1=0.8086165; C2=0.072490; C3=0.05097; C4=0.0094; IMPROV=TRUE; FN = 1.0*N; FK = 1.0*K; if (K > 1) goto e10; // if N = 0 ET K = 1 if (N == 0) return 0.0; // TCHEBYCHEV'S SERIES FOR K <= N else { F1 = pow(FN,1.0/3.0); F2 = F1*F1*FN; ZJP = FN+C1*F1+(C2/F1)-(C3/FN)+(C4/F2); goto e20; } // MAC MAHON'S SERIES FOR K >> N e10: B0 = (FK+0.5*FN-0.75)*PI; B1 = 8.0*B0; B2 = B1*B1; B3 = 3.0*B1*B2; B5 = 5.0*B3*B2; B7 = 7.0*B5*B2; T0 = 4.0*FN*FN; T1 = T0 + 3.0; T3 = 4.0*((7.0*T0+82.0)*T0-9.0); T5 = 32.0*(((83.0*T0+2075.0)*T0-3039.0)*T0+3537.0); T7 = 64.0*((((6949.0*T0+296492.0)*T0-1248002.0)*T0 +7414380.0)*T0-5853627.0); ZJP = B0-(T1/B1)-(T3/B3)-(T5/B5)-(T7/B7); e20: if (IMPROV) { // IMPROVE SOLUTION BY SECANT METHOD WHEN K > N // AND IMPROV = TRUE P0 = 0.9*ZJP; P1 = ZJP; IER = 0; NEV = 2; Q0 = BESSJP(N,P0); Q1 = BESSJP(N,P1); for (IT=1; IT<=NITMX; IT++) { P = P1-Q1*(P1-P0)/(Q1-Q0); DP = P-P1; if (IT == 1) goto e25; if (fabs(DP) < TOL) return P; e25: NEV = NEV+1; P0 = P1; Q0 = Q1; P1 = P ; Q1 = BESSJP(N,P1); } IER = 1; printf(" ** ZEROJP ** NITMX EXCEEDED.\n"); return ZJP; } return ZJP; } double BESSJP(int N, double X) { /* ---------------------------------------------------------------------- ! NAME : BESSJP ! DATE : 06/01/1982 ! IV : 1 ! IE : 1 ! AUTHOR: DANG TRONG TUAN ! ...................................................................... ! ! FIRST DERIVATIVE OF FIRST KIND BESSEL FUNCTION OF ORDER N, FOR REAL X ! ! MODULE BESSJP ! ...................................................................... ! ! THIS SUBROUTINE CALCULATES THE FIRST DERIVATIVE OF FIRST KIND BESSEL ! FUNCTION OF ORDER N, FOR REAL X. ! ! ...................................................................... ! ! I VARIABLE DIMENSION/TYPE DESCRIPTION (INPUTS) ! N I*4 ORDER OF FUNCTION ! X R*8 ABSCISSA OF FUNCTION BESSJP(N,X) ! ! O VARIABLE,DIMENSION/TYPE DESCRIPTION (OUTPUT) ! ! BESSJP R*8 FUNCTION EVALUATION AT X !....................................................................... ! CALLED SUBROUTINE ! ! BESSJ FIRST KIND BESSEL FUNCTION ! ! ----------------------------------------------------------------------*/ if (N == 0) return (-BESSJ(1,X)); else if (X == 0.0) return 0.0; else return (BESSJ(N-1,X)-(1.0*N/X)*BESSJ(N,X)); } double BESSJ(int N, double X) { /*---------------------------------------------------------------------- ! THIS SUBROUTINE CALCULATES THE 1ST KIND BESSEL FUNCTION OF ORDER ! N, INTEGER FOR ANY REAL X. THE CLASSICAL RECURRENCE FORMULA IS USED ! HERE, WHEN X IS > N, FOR X < N, THE MILLER'S ALGORITHM ALLOWS ! AVOIDING OVERFLOWS. ! REFERENCE: ! C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS, ! MATHEMATICAL TABLES, VOL.5, 1962. !---------------------------------------------------------------------*/ int IACC=40; double BIGNO = 1e10, BIGNI = 1e-10; double TOX,BJM,BJ,BJP,BSJ,SUM; int J,JSUM,M; if (N == 0) return BESSJ0(X); if (N == 1) return BESSJ1(X); if (X == 0.0) return 0.0; TOX = 2.0/X; if (X > 1.0*N) { BJM = BESSJ0(X); BJ = BESSJ1(X); for (J=1; J0; J--) { BJM = J*TOX*BJ-BJP; BJP = BJ; BJ = BJM; if (fabs(BJ) > BIGNO) { BJ = BJ*BIGNI; BJP = BJP*BIGNI; BSJ = BSJ*BIGNI; SUM = SUM*BIGNI; } if (JSUM != 0) SUM = SUM + BJ; JSUM = 1 - JSUM; if (J == N) BSJ = BJP; } SUM = 2.0*SUM-BJ; return (BSJ/SUM); } } double BESSJ0(double X) { /*--------------------------------------------------------------------- ! THIS SUBROUTINE CALCULATES THE 1ST KIND BESSEL FUNCTION OF ORDER ! ZERO FOR ANY REAL X. THE CHEBYSHEV'S POLYNOMIAL SERIES IS USED ! FOR 0