!************************************************************** !* Purpose: This program computes the modified Struve * !* function L1(x) using subroutine STVL1 * !* Input : x --- Argument of L1(x) ( x >= 0 ) * !* Output: SL1 --- L1(x) * !* Example: * !* x L1(x) * !* ----------------------- * !* 0.0 .00000000D+00 * !* 5.0 .23728216D+02 * !* 10.0 .26703583D+04 * !* 15.0 .32812429D+06 * !* 20.0 .42454973D+08 * !* 30.0 .76853204D+12 * !* 40.0 .14707396D+17 * !* 50.0 .29030786D+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 MSTVL1 IMPLICIT DOUBLE PRECISION (A-H,O-Z) WRITE(*,*)'Please enter x ' READ(*,*)X WRITE(*,*)' x L1(x)' WRITE(*,*)'-----------------------' CALL STVL1(X,SL1) WRITE(*,10)X,SL1 10 FORMAT(1X,F5.1,D16.8) END SUBROUTINE STVL1(X,SL1) ! ================================================ ! Purpose: Compute modified Struve function L1(x) ! Input : x --- Argument of L1(x) ( x ò 0 ) ! Output: SL1 --- L1(x) ! ================================================ IMPLICIT DOUBLE PRECISION (A-H,O-Z) PI=3.141592653589793D0 R=1.0D0 IF (X.LE.20.0D0) THEN S=0.0D0 DO 10 K=1,60 R=R*X*X/(4.0D0*K*K-1.0D0) S=S+R IF (DABS(R/S).LT.1.0D-12) GO TO 15 10 CONTINUE 15 SL1=2.0D0/PI*S ELSE S=1.0D0 KM=INT(.50*X) IF (X.GT.50) KM=25 DO 20 K=1,KM R=R*(2.0D0*K+3.0D0)*(2.0D0*K+1.0D0)/(X*X) S=S+R IF (DABS(R/S).LT.1.0D-12) GO TO 25 20 CONTINUE 25 SL1=2.0D0/PI*(-1.0D0+1.0D0/(X*X)+3.0D0*S/X**4) A1=DEXP(X)/DSQRT(2.0D0*PI*X) R=1.0D0 BI1=1.0D0 DO 30 K=1,16 R=-0.125D0*R*(4.0D0-(2.0D0*K-1.0D0)**2)/(K*X) BI1=BI1+R IF (DABS(R/BI1).LT.1.0D-12) GO TO 35 30 CONTINUE 35 SL1=SL1+A1*BI1 ENDIF RETURN END !end of file mstvl1.f90