!**************************************************************
!*       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".                                        *
!*                                                            *
!*                          F90 Release By J-P Moreau, Paris. *
!*                                  (www.jpmoreau.fr)         *
!************************************************************** 
        PROGRAM MSTVH1
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        WRITE(*,*)'Please enter x '
        READ(*,*)X
	PRINT *,' '
        WRITE(*,*)'   x       H1(x)'
        WRITE(*,*)'-----------------------'
        CALL STVH1(X,SH1)
        WRITE(*,10)X,SH1
        PRINT *,' ' 
10      FORMAT(1X,F5.1,E16.8)
        END


        SUBROUTINE STVH1(X,SH1)
!       =============================================
!       Purpose: Compute Struve function H1(x)
!       Input :  x   --- Argument of H1(x) ( x ò 0 )
!       Output:  SH1 --- H1(x)
!       =============================================
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        PI=3.141592653589793D0
        R=1.0D0
        IF (X.LE.20.0D0) THEN
           S=0.0D0
           A0=-2.0D0/PI
           DO 10 K=1,60
              R=-R*X*X/(4.0D0*K*K-1.0D0)
              S=S+R
              IF (DABS(R).LT.DABS(S)*1.0D-12) GO TO 15
10         CONTINUE
15         SH1=A0*S
        ELSE
           S=1.0D0
           KM=INT(.5*X)
           IF (X.GT.50.D0) KM=25
           DO 20 K=1,KM
              R=-R*(4.0D0*K*K-1.0D0)/(X*X)
              S=S+R
              IF (DABS(R).LT.DABS(S)*1.0D-12) GO TO 25
20         CONTINUE
25         T=4.0D0/X
           T2=T*T
           P1=((((.42414D-5*T2-.20092D-4)*T2+.580759D-4)*T2   &
              -.223203D-3)*T2+.29218256D-2)*T2+.3989422819D0
           Q1=T*(((((-.36594D-5*T2+.1622D-4)*T2-.398708D-4)*  &
              T2+.1064741D-3)*T2-.63904D-3)*T2+.0374008364D0)
           TA1=X-.75D0*PI
           BY1=2.0D0/DSQRT(X)*(P1*DSIN(TA1)+Q1*DCOS(TA1))
           SH1=2.0/PI*(1.0D0+S/(X*X))+BY1
        ENDIF
        RETURN
        END

!end of file mstvh1.f90