DECLARE SUB SPHJ (N%, X!, NM%, SJ!(), DJ!())
DECLARE FUNCTION ENVJ! (N%, X!)
DECLARE FUNCTION MSTA1% (X!, MP%)
DECLARE FUNCTION MSTA2% (X!, N%, MP%)
'****************************************************************
'*     Purpose: This program computes the spherical Bessel      *
'*              functions jn(x) and jn'(x) using subroutine     *
'*              SPHJ                                            *
'*     Input :  x --- Argument of jn(x)  ( x > 0)               *
'*              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".        *
'*                                                              *
'*                          Basic Release By J-P Moreau, Paris. *
'*                                   (www.jpmoreau.fr)          *
'****************************************************************
'PROGRAM MSPHJ
OPTION BASE 0
DEFINT I-N

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

'end of file msphj.bas

FUNCTION ENVJ (N, X)
ENVJ = .5 * LOG(6.28 * N) - N * LOG(1.36 * X / N)
END FUNCTION

FUNCTION MSTA1 (X, MP)

'      ===================================================
'      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
'      ===================================================
A0 = ABS(X)
N0 = INT(1.1 * A0) + 1
F0 = ENVJ(N0, A0) - MP
N1 = N0 + 5
F1 = ENVJ(N1, A0) - MP
FOR IT = 1 TO 20
    NN = N1 - (N1 - N0) / (1! - F0 / F1)
    F = ENVJ(NN, A0) - MP
    IF ABS(NN - N1) < 1 THEN GOTO 30
    N0 = N1
    F0 = F1
    N1 = NN
    F1 = F
NEXT IT
30 MSTA1 = NN
END FUNCTION

FUNCTION MSTA2 (X, N, MP)

'      ===================================================
'      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
'      ===================================================
A0 = ABS(X)
HMP = .5 * MP
EJN = ENVJ(N, A0)
IF (EJN.LE.HMP) THEN
    OBJ = MP
    N0 = INT(1.1 * A0)
ELSE
    OBJ = HMP + EJN
    N0 = N
END IF
F0 = ENVJ(N0, A0) - OBJ
N1 = N0 + 5
F1 = ENVJ(N1, A0) - OBJ
FOR IT = 1 TO 20
    NN = N1 - (N1 - N0) / (1! - F0 / F1)
    F = ENVJ(NN, A0) - OBJ
    IF ABS(NN - N1) < 1 THEN GOTO 10
    N0 = N1
    F0 = F1
    N1 = NN
    F1 = F
NEXT IT
10 MSTA2 = NN + 10
END FUNCTION

SUB SPHJ (N, X, NM, SJ(), DJ())

'      =======================================================
'      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
'      =======================================================
NM = N
IF ABS(X) < 0! THEN
    FOR K = 0 TO N
        SJ(K) = 0!
        DJ(K) = 0!
    NEXT K
    SJ(0) = 1!
    DJ(1) = .3333333333333333#
    GOTO 20
END IF
SJ(0) = SIN(X) / X
SJ(1) = (SJ(0) - COS(X)) / X
IF N >= 2 THEN
    SA = SJ(0)
    SB = SJ(1)
    M = MSTA1(X, 200)
    IF M < N THEN
        NM = M
    ELSE
        M = MSTA2(X, N, 15)
    END IF
    F0 = 0!
    F1 = 1! - 100
    FOR K = M TO 0 STEP -1
        F = (2! * K + 3!) * F1 / X - F0
        IF K <= NM THEN SJ(K) = F
        F0 = F1
        F1 = F
    NEXT K
    IF ABS(SA) > ABS(SB) THEN CS = SA / F
    IF ABS(SA) <= ABS(SB) THEN CS = SB / F0
    FOR K = 0 TO NM
        SJ(K) = CS * SJ(K)
    NEXT K
END IF
DJ(0) = (COS(X) - SIN(X) / X) / X
FOR K = 1 TO NM
    DJ(K) = SJ(K - 1) - (K + 1!) * SJ(K) / X
NEXT K
 
20 END SUB