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


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