'**************************************************************
'*      Purpose:  This program computes Struve function Hv(x) *
'*                for an arbitrary order using subroutine     *
'*                STVHV                                       *
'*      Input :   v  --- Order of Hv(x)  ( -8.0 to 12.5 )     *
'*                x  --- Argument of Hv(x) ( x > 0 )          *
'*      Output:   HV --- Hv(x)                                *
'*      Example:  x = 10.0                                    *
'*                  v           Hv(x)                         *
'*                -----------------------                     *
'*                  .5     .46402212D+00                      *
'*                 1.5     .14452322D+01                      *
'*                 2.5     .31234632D+01                      *
'*                 3.5     .53730255D+01                      *
'*                 4.5     .72083122D+01                      *
'*                 5.5     .76851132D+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)          *
'**************************************************************
DEFDBL A-H, O-Z
DEFINT I-N

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

V = VMIN
10 CALL STVHV(V, X, HV)

PRINT "  "; V; "  "; HV
V = V + DV
IF V <= VMAX THEN GOTO 10

END


SUB STVHV (V, X, HV)
'       =====================================================
'       Purpose: Compute Struve function Hv(x) with an
'                arbitrary order v
'       Input :  v  --- Order of Hv(x)  ( -8.0 ó v ó 12.5 )
'                x  --- Argument of Hv(x) ( x ò 0 )
'       Output:  HV --- Hv(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
        HV = 0.0
    ELSEIF (V < -1.0) THEN
        HV = (-1) ^ (INT(0.5 - V) - 1) * 1.0E15
    ELSEIF (V = -1.0) THEN
        HV = 2.0 / PI
    END IF
    EXIT SUB
END IF
IF (X <= 20.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) < ABS(S) * 1.0D-12) THEN GOTO 15
    NEXT K
    15 HV = (0.5 * X) ^ (V + 1.0) * S
ELSE
    SA = (0.5 * X) ^ (V - 1.0) / PI
    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 = 4.0 * (U0 + L) ^ 2
        R1 = 1.0
        PU1 = 1.0
        FOR K = 1 TO 12
            R1 = -0.0078125 * R1 * (VT - (4.0 * K - 3.0) ^ 2) * (VT - (4.0 * K - 1.0) ^ 2) / ((2.0 * K - 1.0) * K * X * X)
            PU1 = PU1 + R1
        NEXT K
        QU1 = 1.0
        R2 = 1.0
        FOR K = 1 TO 12
            R2 = -0.0078125 * R2 * (VT - (4.0 * K - 1.0) ^ 2) * (VT - (4.0 * K + 1.0) ^ 2) / ((2.0 * K + 1.0) * K * X * X)
            QU1 = QU1 + R2
        NEXT K
        QU1 = 0.125 * (VT - 1.0) / X * QU1
        IF (L = 0) THEN
            PU0 = PU1
            QU0 = QU1
        END IF
    NEXT L
    T0 = X - (0.5 * U0 + 0.25) * PI
    T1 = X - (0.5 * U0 + 0.75) * PI
    SR = SQR(2.0 / (PI * X))
    BY0 = SR * (PU0 * SIN(T0) + QU0 * COS(T0))
    BY1 = SR * (PU1 * SIN(T1) + QU1 * COS(T1))
    BF0 = BY0
    BF1 = BY1
    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 BYV = BY0
    IF (N = 1) THEN BYV = BY1
    IF (N > 1) THEN BYV = BF
    HV = BYV + S0
END IF
END SUB

'*******************************************
'*           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

'end of file mstvhv.bas