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

DIM SK(250), DK(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 SPHK(N, X, NM, SK(), DK())
        
PRINT
PRINT "  n             kn(x)         kn'(x)"
PRINT "--------------------------------------------"
FOR K = 0 TO NM STEP NS
    PRINT K; TAB(12); SK(K); TAB(28); DK(K)
NEXT K
END

SUB SPHK (N, X, NM, SK(), 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
'   =====================================================
    PI = 3.141592653589793#
    NM = N
    IF X < 0! THEN
        FOR K = 0 TO N
            SK(K) = 1.0E+300
            DK(K) = -1.0E+300
        NEXT K
        GOTO 20
    END IF
    SK(0) = .5 * PI / X * EXP(-X)
    SK(1) = SK(0) * (1! + 1! / X)
    F0 = SK(0)
    F1 = SK(1)
    FOR K = 2 TO N
        F = (2! * K - 1!) * F1 / X + F0
        SK(K) = F
        If ABS(F) > 1.0E+300 Then GoTo 10
        F0 = F1
        F1 = F
    NEXT K
10  NM = K - 1
    DK(0) = -SK(1)
    FOR K = 1 TO NM
        DK(K) = -SK(K - 1) - (K + 1!) / X * SK(K)
    NEXT K
20  End Sub

'end of file msphk.bas