'***************************************************************
'*     Purpose: This program computes the modified spherical   *
'*              Bessel functions of the first kind in(x) and   *
'*              in'(x) using subroutine SPHI                   *
'* ----------------------------------------------------------- *
'*     Input :  x --- Argument of in(x)  (x > 0 )              *
'*              n --- Order of in(x) ( 0 to 250 )              *
'*     Output:  SI(n) --- in(x)                                *
'*              DI(n) --- in'(x)                               *
'*     Example: x = 10.0                                       *
'*                n          in(x)               in'(x)        *
'*              --------------------------------------------   *
'*                0     .1101323287D+04     .9911909633D+03    *
'*                1     .9911909633D+03     .9030850948D+03    *
'*                2     .8039659985D+03     .7500011637D+03    *
'*                3     .5892079640D+03     .5682828129D+03    *
'*                4     .3915204237D+03     .3934477522D+03    *
'*                5     .2368395827D+03     .2494166741D+03    *
'* ----------------------------------------------------------- *
'* 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)            *
'***************************************************************
OPTION BASE 0
DEFINT I-N

DIM SI(250), DI(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 SPHI(N, X, NM, SI(), DI())
        
PRINT
PRINT "  n             in(x)         in'(x)"
PRINT "--------------------------------------------"
FOR K = 0 TO NM STEP NS
    PRINT K, SI(K), DI(K)
NEXT K
END


SUB SPHI (N, X, NM, SI(), DI())

'      Purpose: Compute modified spherical Bessel functions
'               of the first kind, in(x) and in'(x)
'      Input :  x --- Argument of in(x)
'               n --- Order of in(x) ( n = 0,1,2,... )
'      Output:  SI(n) --- in(x)
'               DI(n) --- in'(x)
'               NM --- Highest order computed
'      Routines called:
'               MSTA1 and MSTA2 for computing the starting
'               point for backward recurrence
'      ========================================================

NM = N
IF ABS(X) < 1E-100 THEN
    FOR K = 0 TO N
        SI(K) = 0.0
        DI(K) = 0.0
    NEXT K
    SI(0) = 1.0
    DI(1) = 0.333333333333333
    GOTO 20
END IF
SI(0) = SINH(X) / X
SI(1) = -(SINH(X) / X - COSH(X)) / X
SI0 = SI(0)
IF N >= 2 THEN
    M = MSTA1(X, 200)
    IF M < N THEN
        NM = M
    ELSE
        M = MSTA2(X, N, 15)
    END IF
    F0 = 0.0
    F1 = 1.0 - 100
    FOR K = M TO 0 STEP -1
        F = (2.0 * K + 3.0) * F1 / X + F0
        IF K <= NM THEN SI(K) = F
        F0 = F1
        F1 = F
    NEXT K
    CS = SI0 / F
    FOR K = 0 TO NM
        SI(K) = CS * SI(K)
    NEXT K
END IF
DI(0) = SI(1)
FOR K = 1 TO NM
    DI(K) = SI(K - 1) - (K + 1.0) / X * SI(K)
NEXT K

20 END SUB


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.0 - F0 / F1)
    F = ENVJ(NN, A0) - MP
    IF ABS(NN - N1) < 1 THEN GOTO 20
    N0 = N1
    F0 = F1
    N1 = NN
    F1 = F
NEXT IT
20 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 = 0.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.0 - F0 / F1)
    F = ENVJ(NN, A0) - OBJ
    IF ABS(NN - N1) < 1 THEN GOTO 20
    N0 = N1
    F0 = F1
    N1 = NN
    F1 = F
NEXT IT
20 MSTA2 = NN + 10
END FUNCTION

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

FUNCTION COSH (X)
COSH = (EXP(X) + EXP(-X)) / 2
END FUNCTION

FUNCTION SINH (X)
SINH = (EXP(X) - EXP(-X)) / 2
END FUNCTION

'end of file msphi.bas