/*************************************************************
*      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".                                        *
*                                                            *
*                          C++ Release By J-P Moreau, Paris. *
*                                  (www.jpmoreau.fr)         *
*************************************************************/ 
#include <stdio.h>
#include <math.h>

void STVL1(double X, double *SL1) {
/* ================================================
!  Purpose: Compute modified Struve function L1(x)
!  Input :  x   --- Argument of L1(x) ( x ò 0 )
!  Output:  SL1 --- L1(x)
!  ================================================ */
        double A1,BI1,PI,R,S;
		int K,KM;
	    PI=3.141592653589793;
        R=1.0;
        if (X <= 20.0) {
           S=0.0;
           for (K=1; K<61; K++) {
              R=R*X*X/(4.0*K*K-1.0);
              S=S+R;
              if (fabs(R/S) < 1.0e-12)  goto e15;
           }
e15:       *SL1=2.0/PI*S;
        }
        else {
           S=1.0;
           KM=int(0.50*X);
           if (X > 50.0)  KM=25;
           for (K=1; K<=KM; K++) {
              R=R*(2.0*K+3.0)*(2.0*K+1.0)/(X*X);
              S=S+R;
              if (fabs(R/S) < 1.0e-12)  goto e25;
           }
e25:       *SL1=2.0/PI*(-1.0+1.0/(X*X)+3.0*S/(X*X*X*X));
           A1=exp(X)/sqrt(2.0*PI*X);
           R=1.0;
           BI1=1.0;
           for (K=1; K<17; K++) {
              R=-0.125*R*(4.0-(2.0*K-1.0)*(2.0*K-1.0))/(K*X);
              BI1=BI1+R;
              if (fabs(R/BI1) < 1.0e-12)  goto e35;
           }
e35:       *SL1=*SL1+A1*BI1;
        }
}


void main() {
        
	    double X, SL1;
	    printf("\n Please enter x: ");
        scanf("%lf", &X);
        printf("\n   x        L1(x)\n");
        printf("-----------------------\n");
        
		STVL1(X, &SL1);
        
		printf("%5.1f   %e\n\n", X, SL1);

}

// end of file mstvl1.cpp