/************************************************************** * Purpose: This program computes the modified spherical * * Bessel functions of the first kind in(x) and * * in'(x) using subroutine SPHI * * ----------------------------------------------------------- * * Input : x --- Argument of in(x) * * n --- Order of in(x) ( 0 to 250 ) * * Output: SI(n) --- in(x) * * DI(n) --- in'(x) * * Example: x = 10.0 * * n in(x) in'(x) * * -------------------------------------------- * * 0 .1101323287D+04 .9911909633D+03 * * 1 .9911909633D+03 .9030850948D+03 * * 2 .8039659985D+03 .7500011637D+03 * * 3 .5892079640D+03 .5682828129D+03 * * 4 .3915204237D+03 .3934477522D+03 * * 5 .2368395827D+03 .2494166741D+03 * * ----------------------------------------------------------- * * REFERENCE: "Fortran Routines for Computation of Special * * Functions, * * jin.ece.uiuc.edu/routines/routines.html". * * * * Visual C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * **************************************************************/ #include #include int MSTA1(double, int); int MSTA2(double,int,int); double ENVJ(int N, double X); void SPHI(int N, double X, int *NM, double *SI, double *DI) { /* ======================================================== ! Purpose: Compute modified spherical Bessel functions ! of the first kind, in(x) and in'(x) ! Input : x --- Argument of in(x) ! n --- Order of in(x) ( n = 0,1,2,... ) ! Output: SI(n) --- in(x) ! DI(n) --- in'(x) ! NM --- Highest order computed ! Routines called: ! MSTA1 and MSTA2 for computing the starting ! point for backward recurrence ! ======================================================== */ int K, M; double CS, F, F0, F1, SI0; *NM=N; if (fabs(X) < 1e-100) { for (K=0; K<=N; K++) { SI[K]=0.0; DI[K]=0.0; } SI[0]=1.0; DI[1]=0.333333333333333; return; } SI[0]=sinh(X)/X; SI[1]=-(sinh(X)/X-cosh(X))/X; SI0=SI[0]; if (N >= 2) { M=MSTA1(X,200); if (M < N) *NM=M; else M=MSTA2(X,N,15); F0=0.0; F1=1.0-100; for (K=M; K>-1; K--) { F=(2.0*K+3.0)*F1/X+F0; if (K <= *NM) SI[K]=F; F0=F1; F1=F; } CS=SI0/F; for (K=0; K<=*NM; K++) SI[K] *= CS; } DI[0]=SI[1]; for (K=1; K<=*NM; K++) DI[K]=SI[K-1]-(K+1.0)/X*SI[K]; } int MSTA1(double X, int MP) { /* =================================================== ! Purpose: Determine the starting point for backward ! recurrence such that the magnitude of ! Jn(x) at that point is about 10^(-MP) ! Input : x --- Argument of Jn(x) ! MP --- Value of magnitude ! Output: MSTA1 --- Starting point ! =================================================== */ double A0,F,F0,F1; int IT,NN,N0,N1; A0=fabs(X); N0=floor(1.1*A0)+1; F0=ENVJ(N0,A0)-MP; N1=N0+5; F1=ENVJ(N1,A0)-MP; for (IT=1; IT<=20; IT++) { NN=N1-(N1-N0)/(1.0-F0/F1); F=ENVJ(NN,A0)-MP; if (abs(NN-N1) < 1) goto e20; N0=N1; F0=F1; N1=NN; F1=F; } e20: return NN; } int MSTA2(double X, int N, int MP) { /* =================================================== ! Purpose: Determine the starting point for backward ! recurrence such that all Jn(x) has MP ! significant digits ! Input : x --- Argument of Jn(x) ! n --- Order of Jn(x) ! MP --- Significant digit ! Output: MSTA2 --- Starting point ! =================================================== */ double A0,EJN,F,F0,F1,HMP,OBJ; int IT,N0,N1,NN; A0=fabs(X); HMP=0.5*MP; EJN=ENVJ(N,A0); if (EJN <= HMP) { OBJ=MP; N0=floor(1.1*A0); } else { OBJ=HMP+EJN; N0=N; } F0=ENVJ(N0,A0)-OBJ; N1=N0+5; F1=ENVJ(N1,A0)-OBJ; for (IT=1; IT<=20; IT++) { NN=N1-(N1-N0)/(1.0-F0/F1); F=ENVJ(NN,A0)-OBJ; if (abs(NN-N1) < 1) goto e20; N0=N1; F0=F1; N1=NN; F1=F; } e20: return NN+10; } double ENVJ(int N, double X) { return (0.5*log(6.28*N)-N*log(1.36*X/N)); } void main() { double SI[251], DI[251]; int k, n, nm, ns; double X; printf("\n Please enter n and x: "); scanf("%d %lf", &n, &X); if (n <= 10) ns=1; else { printf(" Please enter order step Ns: "); scanf("%d", &ns); } SPHI(n,X,&nm,SI,DI); printf("\n"); printf(" n in(x) in''(x)\n"); printf("--------------------------------------------\n"); for (k=0; k<=nm; k+=ns) printf(" %d %18.6f %18.6f\n", k, SI[k], DI[k]); } // end of file msphi.cpp