/****************************************************************** * Purpose: This program computes the modified spherical * * Bessel functions kn(x) and kn'(x) using * * subroutine SPHK * * Input : x --- Argument of kn(x) ( x > 0 ) * * n --- Order of kn(x) ( n = 0 to 250 ) * * Output: SK(n) --- kn(x) * * DK(n) --- kn'(x) * * Example: x= 10.0 * * n kn(x) kn'(x) * * -------------------------------------------- * * 0 .7131404291D-05 -.7844544720D-05 * * 1 .7844544720D-05 -.8700313235D-05 * * 2 .9484767707D-05 -.1068997503D-04 * * 3 .1258692857D-04 -.1451953914D-04 * * 4 .1829561771D-04 -.2173473743D-04 * * 5 .2905298451D-04 -.3572740841D-04 * * --------------------------------------------------------------- * * 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 void SPHK(int N, double X, int *NM, double *SK, double *DK) { /* ===================================================== ! Purpose: Compute modified spherical Bessel functions ! of the second kind, kn(x) and kn'(x) ! Input : x --- Argument of kn(x) ( x ò 0 ) ! n --- Order of kn(x) ( n = 0,1,2,... ) ! Output: SK(n) --- kn(x) ! DK(n) --- kn'(x) ! NM --- Highest order computed ! ===================================================== */ const double PI=3.141592653589793; int K; double F, F0, F1; *NM=N; if (X < 1e-60) { for (K=0; K<=N; K++) { SK[K]=1.0e+300; DK[K]=-1.0e+300; } return; } SK[0]=0.5*PI/X*exp(-X); SK[1]=SK[0]*(1.0+1.0/X); F0=SK[0]; F1=SK[1]; for (K=2; K<=N; K++) { F=(2.0*K-1.0)*F1/X+F0; SK[K]=F; if (fabs(F) > 1.0e+300) goto e20; F0=F1; F1=F; } e20: *NM=K-1; DK[0]=-SK[1]; for (K=1; K<=*NM; K++) DK[K]=-SK[K-1]-(K+1.0)/X*SK[K]; } void main() { double SK[251], DK[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); } SPHK(n,X,&nm,SK,DK); printf("\n"); printf(" n kn(x) kn''(x)\n"); printf("--------------------------------------------\n"); for (k=0; k<=nm; k+=ns) printf(" %d %e %e\n", k, SK[k], DK[k]); } //end of file msphk.cpp