'**************************************************************
'*      Purpose:  This program computes the modified Struve   *
'*                function Lv(x) for an arbitrary order v     *
'*                using subroutine STVLV                      *
'*      Input :   v   --- Order of Lv(x)  ( |v| ó 20 )        *
'*                x   --- Argument of Lv(x) ( x ò 0 )         *
'*      Output:   SLV --- Lv(x)                               *
'*      Example:  x = 10.0                                    *
'*                  v          Lv(x)                          *
'*                ------------------------                    *
'*                 0.5     .27785323D+04                      *
'*                 1.5     .24996698D+04                      *
'*                 2.5     .20254774D+04                      *
'*                 3.5     .14816746D+04                      *
'*                 4.5     .98173460D+03                      *
'*                 5.5     .59154277D+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)          *
'**************************************************************
DEFDBL A-H, O-Z
DEFINT I-N

PRINT
INPUT " Please enter vmin, vmax, dv and x: ", VMIN, VMAX, DV, X
        
PRINT
PRINT "   v          Lv(x)"
PRINT " -------------------------"
        
V = VMIN

10 CALL STVLV(V, X, SLV)
PRINT " "; V; "  "; SLV
V = V + DV
IF (V <= VMAX) THEN GOTO 10

PRINT " -------------------------"

END

'*******************************************
'*           FUNCTION  GAMMA(X)            *
'* --------------------------------------- *
'* Returns the value of Gamma(x) in double *
'* precision as EXP(LN(GAMMA(X))) for X>0. *
'* (input: x  output: y)                   *
'*******************************************
FUNCTION Gamma (x)
DIM cof(6)
cof(1) = 76.18009173#: cof(2) = -86.50532033#
cof(3) = 24.01409822#: cof(4) = -1.231739516#
cof(5) = 0.00120858003#: cof(6) = -0.00000536382#
stp = 2.50662827465#
half = 0.5#: one = 1#: fpf = 5.5#
xx = x - one
tmp = xx + fpf
tmp = (xx + half) * LOG(tmp) - tmp
ser = one
FOR j = 1 TO 6
    xx = xx + one
    ser = ser + cof(j) / xx
NEXT j
Gamma = EXP(tmp + LOG(stp * ser))
END FUNCTION


SUB STVLV (V, X, SLV)
'   ======================================================
'   Purpose:  Compute modified Struve function Lv(x) with
'             an arbitrary order v
'   Input :   v   --- Order of Lv(x)  ( |v| ó 20 )
'             x   --- Argument of Lv(x) ( x ò 0 )
'   Output:   SLV --- Lv(x)
'   Routine called: GAMMA to compute the gamma function
'   ======================================================
PI = 3.141592653589793#
IF (X = 0.0) THEN
    IF (V > -1.0 OR INT(V) - V = 0.5) THEN
        SLV = 0.0
    ELSEIF (V < -1.0) THEN
        SLV = (-1) ^ (INT(0.5 - V) - 1) * 1.0D+15
    ELSEIF (V = -1.0) THEN
        SLV = 2.0 / PI
    END IF
    EXIT SUB
END IF
IF (X <= 40.0) THEN
    V0 = V + 1.5
    GA = Gamma(V0)
    S = 2.0 / (SQR(PI) * GA)
    R1 = 1.0
    FOR K = 1 TO 100
        VA = K + 1.5
        GA = Gamma(VA)
        VB = V + K + 1.5
        GB = Gamma(VB)
        R1 = R1 * (0.5 * X) ^ 2
        R2 = R1 / (GA * GB)
        S = S + R2
        IF (ABS(R2 / S) < 1.0D-12) THEN GOTO 15
    NEXT K
    15 SLV = (0.5 * X) ^ (V + 1.0) * S
ELSE
    SA = -1.0 / PI * (0.5 * X) ^ (V - 1.0)
    V0 = V + 0.5
    GA = Gamma(V0)
    S = -SQR(PI) / GA
    R1 = -1.0
    FOR K = 1 TO 12
        VA = K + 0.5
        GA = Gamma(VA)
        VB = -K + V + 0.5
        GB = Gamma(VB)
        R1 = -R1 / (0.5 * X) ^ 2
        S = S + R1 * GA / GB
    NEXT K
    S0 = SA * S
    U = ABS(V)
    N = INT(U)
    U0 = U - N
    FOR L = 0 TO 1
        VT = U0 + L
        R = 1.0
        BIV = 1.0
        FOR K = 1 TO 16
            R = -0.125 * R * (4.0 * VT * VT - (2.0 * K - 1.0) ^ 2) / (K * X)
            BIV = BIV + R
            IF (ABS(R / BIV) < 1.0D-12) THEN GOTO 30
        NEXT K
        30 IF (L = 0) THEN BIV0 = BIV
    NEXT L
    BF0 = BIV0
    BF1 = BIV
    FOR K = 2 TO N
        BF = -2.0 * (K - 1.0 + U0) / X * BF1 + BF0
        BF0 = BF1
        BF1 = BF
    NEXT K
    IF (N = 0) THEN BIV = BIV0
    IF (N > 1) THEN BIV = BF
    SLV = EXP(X) / SQR(2.0 * PI * X) * BIV + S0
END IF

END SUB


'end of file mstvlv.bas