/**************************************************************
*     Purpose: This program computes the spherical Bessel      *
*              functions jn(x) and jn'(x) using subroutine     *
*              SPHJ                                            *
*     Input :  x --- Argument of jn(x)                         *
*              n --- Order of jn(x)  ( n = 0 to 250 )          * 
*     Output:  SJ(n) --- jn(x)                                 *
*              DJ(n) --- jn'(x)                                * 
*     Example:   x =10.0                                       *
*                n          jn(x)              jn'(x)          *
*              --------------------------------------------    *
*                0    -.5440211109D-01    -.7846694180D-01     * 
*                1     .7846694180D-01    -.7009549945D-01     *
*                2     .7794219363D-01     .5508428371D-01     *
*                3    -.3949584498D-01     .9374053162D-01     *
*                4    -.1055892851D+00     .1329879757D-01     *
*                5    -.5553451162D-01    -.7226857814D-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>

int MSTA1(double, int);
int MSTA2(double,int,int);
double ENVJ(int N, double X);



void SPHJ(int N, double X, int *NM, double *SJ, double *DJ) {
/*     =======================================================
!      Purpose: Compute spherical Bessel functions jn(x) and
!               their derivatives
!      Input :  x --- Argument of jn(x)
!               n --- Order of jn(x)  ( n = 0,1,תתת )
!      Output:  SJ(n) --- jn(x)
!               DJ(n) --- jn'(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, SA, SB;
        
	*NM=N;
        if (fabs(X) < 1e-100) {
	   for (K=0; K<=N; K++) {
              SJ[K]=0.0;
              DJ[K]=0.0;
           } 
           SJ[0]=1.0;
           DJ[1]=0.333333333333333;
           return;
        }
        SJ[0]=sin(X)/X;
        SJ[1]=(SJ[0]-cos(X))/X;
        if (N >= 2) {
		   SA=SJ[0];
           SB=SJ[1];
           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)  SJ[K]=F;
              F0=F1;
              F1=F;
           }
           if (fabs(SA) > fabs(SB))  CS=SA/F;
           if (fabs(SA) <= fabs(SB)) CS=SB/F0;
		   for (K=0; K<=*NM; K++)  SJ[K] *= CS;
        }      
        DJ[0]=(cos(X)-sin(X)/X)/X;
		for (K=1; K<=*NM; K++)
          DJ[K]=SJ[K-1]-(K+1.0)*SJ[K]/X;
}

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 SJ[251], DJ[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);
        }
        
	SPHJ(n,X,&nm,SJ,DJ);
        
	printf("\n");
        printf("  n          jn(x)              jn''(x)\n");
        printf("--------------------------------------------\n");
        for (k=0; k<=nm; k+=ns)
          printf(" %d %18.6f %18.6f\n", k, SJ[k], DJ[k]);
}

//end of file msphj.cpp