!**************************************************************
!*      Purpose: This program computes modified Struve        *
!*               function L0(x) using subroutine STVL0        *
!*      Input :  x   --- Argument of L0(x) ( x ò 0 )          *
!*      Output:  SL0 --- L0(x)                                *
!*      Example:                                              *
!*                  x        L0(x)                            *
!*              ------------------------                      *
!*                 0.0   .00000000D+00                        *
!*                 5.0   .27105917D+02                        *
!*                10.0   .28156522D+04                        *
!*                15.0   .33964933D+06                        *
!*                20.0   .43558283D+08                        *
!*                30.0   .78167230D+12                        *
!*                40.0   .14894775D+17                        *
!*                50.0   .29325538D+21                        *
!* ---------------------------------------------------------- *
!* 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 MSTVL0
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        WRITE(*,*)'Please enter x '
        READ(*,*) X
		print *,' '
        WRITE(*,*)'   x        L0(x)'
        WRITE(*,*)'-----------------------'
        CALL STVL0(X,SL0)
        WRITE(*,10)X,SL0
10      FORMAT(1X,F5.1,D16.8/)
        END


        SUBROUTINE STVL0(X,SL0)
!       ================================================
!       Purpose: Compute modified Struve function L0(x)
!       Input :  x   --- Argument of L0(x) ( x ò 0 )
!       Output:  SL0 --- L0(x)
!       ================================================
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        PI=3.141592653589793D0
        S=1.0D0
        R=1.0D0
        IF (X.LE.20.0D0) THEN
           A0=2.0D0*X/PI
           DO 10 K=1,60
              R=R*(X/(2.0D0*K+1.0D0))**2
              S=S+R
              IF (DABS(R/S).LT.1.0D-12) GO TO 15
10         CONTINUE
15         SL0=A0*S
        ELSE
           KM=INT(.5*(X+1.0))
           IF (X.GE.50.0) KM=25
           DO 20 K=1,KM
              R=R*((2.0D0*K-1.0D0)/X)**2
              S=S+R
              IF (DABS(R/S).LT.1.0D-12) GO TO 25
20         CONTINUE
25         A1=DEXP(X)/DSQRT(2.0D0*PI*X)
           R=1.0D0
           BI0=1.0D0
           DO 30 K=1,16
              R=0.125D0*R*(2.0D0*K-1.0D0)**2/(K*X)
              BI0=BI0+R
              IF (DABS(R/BI0).LT.1.0D-12) GO TO 35
30         CONTINUE
35         BI0=A1*BI0
           SL0=-2.0D0/(PI*X)*S+BI0
        ENDIF
        RETURN
        END

!end of file mstvl0.f90