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

Uses WinCrt;

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

Var
     SY, DY: VEC;
     k, n, nm, ns: Integer;
     X: double;


Procedure SPHY(N:integer; X:double; Var NM:integer; Var SY, DY: VEC);
{       ======================================================
!       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
!       ====================================================== }
Label 10, 20;
Var
        K:Integer; F, F0, F1:Double;
Begin
	NM:=N;
        if X < 1E-60 then
        begin
           for K:=0 to N do
           begin
             SY[K]:=-1E+60;
             DY[K]:=1E+60
           end;
           goto 10
        end;
        SY[0]:=-cos(X)/X;
        SY[1]:=(SY[0]-sin(X))/X;
        F0:=SY[0];
        F1:=SY[1];
        for K:=2 to N do
        begin
           F:=(2.0*K-1.0)*F1/X-F0;
           SY[K]:=F;
           if abs(F) >= 1E+60 then goto 20;              
           F0:=F1;
           F1:=F
        end;
20:     NM:=K-1;
        DY[0]:=(sin(X)+cos(X)/X)/X;
        for K:=1 to NM do
           DY[K]:=SY[K-1]-(K+1.0)*SY[K]/X;
10: End;


{main program}
BEGIN
        InitWinCrt;
        writeln;
        write(' Please enter n and x: ');
        readln(n, X);

        Inc(n);

	if n <= 10 then
           ns:=1
        else
        begin
           write(' Please enter order step Ns: ');
           readln(ns)
        end;
        
	SPHY(n,X,nm,SY,DY);
        
	writeln;
        writeln('  n          yn(x)             yn''(x)');
        writeln('--------------------------------------------');
        k:=0;
        Repeat 
          writeln(k:3, SY[k]:18:6, DY[k]:18:6);
          k:=k+ns
        until k = nm+1;
        ReadKey;
        DoneWinCrt
END.

{end of file msphy.pas}