DECLARE SUB STVH1 (X#, SH1#)
'**************************************************************
'*       Purpose: This program computes Struve function       *
'*                H1(x) using subroutine STVH1                *
'*       Input :  x   --- Argument of H1(x) ( x ò 0 )         *
'*       Output:  SH1 --- H1(x)                               *
'*       Example:                                             *
'*                   x          H1(x)                         *
'*                -----------------------                     *
'*                  0.0       .00000000                       *
'*                  5.0       .80781195                       *
'*                 10.0       .89183249                       *
'*                 15.0       .66048730                       *
'*                 20.0       .47268818                       *
'*                 25.0       .53880362                       *
'* ---------------------------------------------------------- *
'* 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          H1(x)         "
	PRINT "----------------------------"

	CALL STVH1(X, SH1)
	
	PRINT "  "; X; "  "; SH1

	END

SUB STVH1 (X, SH1)
'       =============================================
'       Purpose: Compute Struve function H1(x)
'       Input :  x   --- Argument of H1(x) ( x > 0 )
'       Output:  SH1 --- H1(x)
'       =============================================
	PI = 3.141592653589793#
	R = 1#
	IF (X <= 20!) THEN
	   S = 0!
	   A0 = -2! / PI
	   FOR K = 1 TO 60
	      R = -R * X * X / (4! * K * K - 1#)
	      S = S + R
	      IF (ABS(R) < ABS(S) * .000000000001#) THEN GOTO 15
	   NEXT K
15         SH1 = A0 * S
	ELSE
	   S = 1!
	   KM = INT(.5 * X)
	   IF (X > 50!) THEN KM = 25
	   FOR K = 1 TO KM
	      R = -R * (4! * K * K - 1!) / (X * X)
	      S = S + R
	      IF (ABS(R) < ABS(S) * 1E-12) THEN GOTO 25
	   NEXT K
25 :       T = 4# / X
	   T2 = T * T
	   P1 = ((((4.2414E-06 * T2 - 2.0092E-05) * T2 + 5.80759E-05) * T2 - 2.23203E-04) * T2 + .0029218256#) * T2 + .3989422819#
	   Q1 = T * (((((-3.6594E-06 * T2 + 1.622E-05) * T2 - 3.98708E-05) * T2 + .0001064741#) * T2 - 6.3904E-04) * T2 + .0374008364#)
	   TA1 = X - .75 * PI
	   BY1 = 2! / SQR(X) * (P1 * SIN(TA1) + Q1 * COS(TA1))
	   SH1 = 2! / PI * (1# + S / (X * X)) + BY1
	END IF
END SUB