{***************************************************************
*     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".        *
*                                                              *
*                         Pascal Release By J-P Moreau, Paris. *
*                                 (www.jpmoreau.fr)            *
***************************************************************}
Program MSPHJ;

Uses WinCrt;

Type VEC = Array[0..250] of double;

Const MM = 2.30258509;

Var
     SJ, DJ: VEC;
     k, n, nm, ns: Integer;
     X: double;

Function ENVJ(N:integer; X:double): double;
Begin
  ENVJ := (0.5*Ln(6.28*N)-N*Ln(1.36*X/N))/MM
End;

Function MSTA1(X: double; MP: integer): Integer;

{      ===================================================
!      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   
!      =================================================== }
Label 20;
Var
        A0,F,F0,F1:Double; IT,NN,N0,N1:Integer;
Begin
        A0:=abs(X);
        N0:=Round(1.1*A0)+1;
        F0:=ENVJ(N0,A0)-MP;
        N1:=N0+5;
        F1:=ENVJ(N1,A0)-MP;
        for IT:=1 to 20 do
        begin             
           NN:=Round(N1-(N1-N0)/(1.0-F0/F1));                  
           F:=ENVJ(NN,A0)-MP;
           if abs(NN-N1) < 1 then goto 20;
           N0:=N1;
           F0:=F1;
           N1:=NN;
           F1:=F
        end;
20:	MSTA1 := NN
End;

Function MSTA2(X:double; N, MP: integer): Integer;

{      ===================================================
!      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
!      =================================================== }
Label 20;
Var
        A0,EJN,F,F0,F1,HMP,OBJ:Double; IT,N0,N1,NN: Integer;
Begin
        A0:=abs(X);
        HMP:=0.5*MP;
        EJN:=ENVJ(N,A0);
        if EJN <= HMP then
        begin
           OBJ:=MP;
           N0:=Round(1.1*A0)
        end
        else
        begin
           OBJ:=HMP+EJN;
           N0:=N
        end;
        F0:=ENVJ(N0,A0)-OBJ;
        N1:=N0+5;
        F1:=ENVJ(N1,A0)-OBJ;
        for IT:=1 to 20 do
        begin
           NN:=Round(N1-(N1-N0)/(1.0-F0/F1));
           F:=ENVJ(NN,A0)-OBJ;
           if abs(NN-N1) < 1 then goto 20;
           N0:=N1;
           F0:=F1;
           N1:=NN;
           F1:=F
        end;
20:     MSTA2 := NN+10
End;

Procedure SPHJ(N:integer; X:double; Var NM:integer; Var SJ, DJ: VEC);
{      =======================================================
!      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
!      ======================================================= }
Label   10;
Var
        K, M: Integer; CS, F, F0, F1, SA, SB: Double;
Begin
        NM:=N;
        if abs(X) < 1e-100 then
        begin
	   for K:=0 to N do
           begin
             SJ[K]:=0.0;
             DJ[K]:=0.0
           end; 
           SJ[0]:=1.0;
           DJ[1]:=0.333333333333333;
           goto 10;
        end;
        SJ[0]:=sin(X)/X;
        SJ[1]:=(SJ[0]-cos(X))/X;
        if N >= 2 then
        begin
	   SA:=SJ[0];
           SB:=SJ[1];
           M:=MSTA1(X,200);
           if M < N then
              NM:=M
           else
              M:=MSTA2(X,N,15);
           F0:=0.0;
           F1:=1.0-100;
           for K:=M Downto 0 do
           begin
              F:=(2.0*K+3.0)*F1/X-F0;
              if K <= NM then SJ[K]:=F;
              F0:=F1;
              F1:=F
           end;
           if abs(SA) > abs(SB) then CS:=SA/F;
           if abs(SA) <= abs(SB) then CS:=SB/F0;
	   for K:=0 to NM do  SJ[K] := SJ[K]*CS
        end;      
        DJ[0]:=(cos(X)-sin(X)/X)/X;
	for K:=1 to NM do
          DJ[K]:=SJ[K-1]-(K+1.0)*SJ[K]/X;
10:End;


{main program}
BEGIN
        InitWinCrt;
        writeln;
        write(' Please enter n and x: ');
        readln(n, X);
        
	if n <= 10 then
           ns:=1
        else
        begin
           write(' Please enter order step Ns: ');
           readln(ns)
        end;
        
	SPHJ(n,X,nm,SJ,DJ);
        
	writeln;
        writeln('  n          jn(x)             jn''(x)');
        writeln('--------------------------------------------');
        k:=0;
        Repeat 
          writeln(k:3, SJ[k]:18:6, DJ[k]:18:6);
          k:=k+ns
        until k = nm+1;
        ReadKey;
        DoneWinCrt
END.

{end of file msphj.pas}