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

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