/**************************************************************
*      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 <stdio.h>
#include <math.h>


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