DECLARE SUB SPHK (N%, X!, NM%, SK!(), DK!())
'***************************************************************
'*      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".       *
'*                                                             *
'*                         Basic Release By J-P Moreau, Paris. *
'*                                 (www.jpmoreau.fr)           *
'***************************************************************
'PROGRAM MSPHY
OPTION BASE 0
DEFINT I-N

DIM SY(250), DY(250)
PRINT
INPUT " Please enter n and x: ", N, X
IF N <= 10 THEN
    NS = 1
ELSE
    INPUT " Please enter order step Ns: ", NS
END IF
        
CALL SPHY(N, X, NM, SY(), DY())
        
PRINT
PRINT "  n             yn(x)         yn'(x)"
PRINT "--------------------------------------------"
FOR K = 0 TO NM STEP NS
    PRINT K; TAB(12); SY(K); TAB(28); DY(K)
NEXT K
END


Sub SPHY(N, X, NM, SY(), DY())

    '   ======================================================
    '    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
    '   ======================================================
    NM = N
    If X < 1.0E-60 Then
        For K = 0 To N
            SY(K) = -1.0E+300
            DY(K) = 1.0E+300
        Next K
        GoTo 20
    End If
    SY(0) = -DCOS(X) / X
    SY(1) = (SY(0) - DSIN(X)) / X
    F0 = SY(0)
    F1 = SY(1)
    For K = 2 To N
        F = (2.0 * K - 1.0) * F1 / X - F0
        SY(K) = F
      IF ABS(F) >= 1E+300) Then GoTo 10              
            F0 = F1
            F1 = F
    NEXT K
10: NM = K - 1
    DY(0) = (SIN(X) + COS(X) / X) / X
    For K = 1 To NM
      DY(K)=SY(K-1)-(K+1.0D0)*SY(K)/X
    Next K
20 End Sub