/************************************************************** * Purpose: This program computes the spherical Bessel * * functions yn(x) and yn'(x) using subroutine * * SPHY * * Input : x --- Argument of yn(x) ( x > 0 ) * * n --- Order of yn(x) ( n = 0 to 250 ) * * Output: SY(n) --- yn(x) * * DY(n) --- yn'(x) * * Example: x = 10.0 * * n yn(x) yn'(x) * * -------------------------------------------- * * 0 .8390715291D-01 -.6279282638D-01 * * 1 .6279282638D-01 .7134858763D-01 * * 2 -.6506930499D-01 .8231361788D-01 * * 3 -.9532747888D-01 -.2693831344D-01 * * 4 -.1659930220D-02 -.9449751377D-01 * * 5 .9383354168D-01 -.5796005523D-01 * * ----------------------------------------------------------- * * 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 SPHY(int N, double X, int *NM, double *SY, double *DY) { /* ====================================================== ! Purpose: Compute spherical Bessel functions yn(x) and ! their derivatives ! Input : x --- Argument of yn(x) ( x ע 0 ) ! n --- Order of yn(x) ( n = 0,1,תתת ) ! Output: SY(n) --- yn(x) ! DY(n) --- yn'(x) ! NM --- Highest order computed ! ====================================================== */ int K; double F, F0, F1; *NM=N; if (X < 1e-60) { for (K=0; K<=N; K++) { SY[K]=-1.0e+300; DY[K]=1e+300; } return; } SY[0]=-cos(X)/X; SY[1]=(SY[0]-sin(X))/X; F0=SY[0]; F1=SY[1]; for (K=2; K<=N; K++) { F=(2.0*K-1.0)*F1/X-F0; SY[K]=F; if (fabs(F) >= 1e+300) goto e20; F0=F1; F1=F; } e20: *NM=K-1; DY[0]=(sin(X)+cos(X)/X)/X; for (K=1; K<=*NM; K++) DY[K]=SY[K-1]-(K+1.0)*SY[K]/X; } void main() { double SY[251], DY[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); } SPHY(n,X,&nm,SY,DY); printf("\n"); printf(" n yn(x) yn''(x)\n"); printf("--------------------------------------------\n"); for (k=0; k<=nm; k+=ns) printf(" %d %e %e\n", k, SY[k], DY[k]); } //end of file msphy.cpp