'**************************************************************
'*      Purpose: This program computes the modified Struve    *
'*               function L1(x) using subroutine STVL1        *
'*      Input :  x   --- Argument of L1(x) ( x >= 0 )         *
'*      Output:  SL1 --- L1(x)                                *
'*      Example:                                              *
'*                    x        L1(x)                          *
'*                -----------------------                     *
'*                   0.0   .00000000D+00                      *
'*                   5.0   .23728216D+02                      *
'*                  10.0   .26703583D+04                      *
'*                  15.0   .32812429D+06                      *
'*                  20.0   .42454973D+08                      *
'*                  30.0   .76853204D+12                      *
'*                  40.0   .14707396D+17                      *
'*                  50.0   .29030786D+21                      *
'* ---------------------------------------------------------- *
'* 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)          *
'**************************************************************
DEFDBL A-H, O-Z
DEFINT I-N

PRINT
INPUT " Please enter x: ", X
PRINT
PRINT "   x        L1(x)"
PRINT "-------------------------"
        
CALL STVL1(X, SL1)
        
PRINT " "; X; " "; SL1

END


SUB STVL1 (X, SL1)
'       ================================================
'       Purpose: Compute modified Struve function L1(x)
'       Input :  x   --- Argument of L1(x) ( x ò 0 )
'       Output:  SL1 --- L1(x)
'       ================================================

PI = 3.141592653589793#
R = 1.0
IF (X <= 20.0) THEN
    S = 0.0
    FOR K = 1 TO 60
        R = R * X * X / (4.0 * K * K - 1.0)
        S = S + R
        IF (ABS(R / S) < 1.0D-12) THEN GOTO 15
    NEXT K
    15 SL1 = 2.0 / PI * S
ELSE
    S = 1.0
    KM = INT(0.50 * X)
    IF (X > 50) THEN KM = 25
    FOR K = 1 TO KM
        R = R * (2.0 * K + 3.0) * (2.0 * K + 1.0) / (X * X)
        S = S + R
        IF (ABS(R / S) < 1.0D-12) THEN GOTO 25
    NEXT K
    25 SL1 = 2.0 / PI * (-1.0 + 1.0 / (X * X) + 3.0 * S / X ^ 4)
    A1 = EXP(X) / SQR(2.0 * PI * X)
    R = 1.0
    BI1 = 1.0
    FOR K = 1 TO 16
        R = -0.125 * R * (4.0 - (2.0 * K - 1.0) ^ 2) / (K * X)
        BI1 = BI1 + R
        IF (ABS(R / BI1) < 1.0D-12) THEN GOTO 35
    NEXT K
    35 SL1 = SL1 + A1 * BI1
END IF
END SUB

'end of file mstvl1.bas