/**************************************************************
!*      Purpose:  This program computes the modified Struve   *
!*                function Lv(x) for an arbitrary order v     * 
!*                using subroutine STVLV                      *
!*      Input :   v   --- Order of Lv(x)  ( |v| ó 20 )        *
!*                x   --- Argument of Lv(x) ( x ò 0 )         *
!*      Output:   SLV --- Lv(x)                               *
!*      Example:  x = 10.0                                    * 
!*                  v          Lv(x)                          *
!*                ------------------------                    *
!*                 0.5     .27785323D+04                      * 
!*                 1.5     .24996698D+04                      *
!*                 2.5     .20254774D+04                      *
!*                 3.5     .14816746D+04                      *
!*                 4.5     .98173460D+03                      *
!*                 5.5     .59154277D+03                      *
!*                ------------------------                    *
!*                                                            *
!* ---------------------------------------------------------- *
!* 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>

double GAMMA(double X) {

/*  ==================================================
!   Purpose: Compute gamma function gamma(x)
!   Input :  x  --- Argument of gamma(x)
!                   ( x is not equal to 0,-1,-2)
!   Output:  gamma(x)
!   ================================================== */
        double GA, GR,PI, R, Z;
        double G[27];
        int K, M, M1;

        PI=3.141592653589793;
        if (X==int(X))
	   if (X>0.0) {
              GA=1.0;
              M1=(int) (X-1);
              for (K=2; K<=M1; K++)  GA=GA*K;
	   }
           else
              GA=1e+100;
	else {
	   if (fabs(X)>1.0) {
              Z=fabs(X);
              M=int(Z);
              R=1.0;
              for (K=1; K<=M; K++)  R=R*(Z-K);
              Z=Z-M;
           }
           else
              Z=X;

           G[1]=1.0; G[2]=0.5772156649015329;
	       G[3]=-0.6558780715202538; G[4]=-0.420026350340952e-1;
	       G[5]=0.1665386113822915; G[6]=-0.421977345555443e-1;
	       G[7]=-0.96219715278770e-2; G[8]=0.72189432466630e-2;
	       G[9]=-0.11651675918591e-2; G[10]=-0.2152416741149e-3;
           G[11]=0.1280502823882e-3; G[12]=-0.201348547807e-4;
	       G[13]=-0.12504934821e-5; G[14]=0.11330272320e-5;
	       G[15]=-0.2056338417e-6; G[16]=0.61160950e-8;
	       G[17]=0.50020075e-8; G[18]=-0.11812746e-8;
	       G[19]=0.1043427e-9; G[20]=0.77823e-11;
	       G[21]=-0.36968e-11; G[22]=0.51e-12;
	       G[23]=-0.206e-13; G[24]=-0.54e-14;
	       G[25]=0.14e-14; G[26]=0.1e-15;

           GR=G[26];
           for (K=25; K>0; K--) GR=GR*Z+G[K];
           GA=1.0/(GR*Z);
           if (fabs(X)>1.0) {
              GA=GA*R;
              if (X < 0.0)  GA=-PI/(X*GA*sin(PI*X));
           }
        }
        return GA;
}


void STVLV(double V, double X, double *SLV) {
/*      ======================================================
!       Purpose:  Compute modified Struve function Lv(x) with
!                 an arbitrary order v
!       Input :   v   --- Order of Lv(x)  ( |v| ó 20 )
!                 x   --- Argument of Lv(x) ( x ò 0 )
!       Output:   SLV --- Lv(x)
!       Routine called: GAMMA to compute the gamma function
!       ====================================================== */
        double BIV, BIV0, GA, GB, PI, R, R1, R2, S, S0, SA, U, U0, V0, VA, VB, VT;
		double BF, BF0, BF1;
		int K, L, N;
	    PI=3.141592653589793;
        if (X == 0.0) {
           if (V > -1.0 || int(V)-V == 0.5)
              *SLV=0.0;
           else if (V < -1.0)
              *SLV=pow(-1, int(0.5-V)-1)*1.0e+300;
           else if (V == -1.0)
              *SLV=2.0/PI;
		   return;
        }
        if (X <= 40.0) {
           V0=V+1.5;
           GA = GAMMA(V0);
           S=2.0/(sqrt(PI)*GA);
           R1=1.0;
           for (K=1; K<=100; K++) {
              VA=K+1.5;
              GA = GAMMA(VA);
              VB=V+K+1.5;
              GB = GAMMA(VB);
              R1=R1*(0.5*X)*(0.5*X);
              R2=R1/(GA*GB);
              S=S+R2;
              if (fabs(R2/S) < 1.0e-12) goto e15;
	   }
e15:       *SLV=pow(0.5*X,V+1.0)*S;
        }
        else {
           SA=-1.0/PI*pow(0.5*X, V-1.0);
           V0=V+0.5;
           GA = GAMMA(V0);
           S=-sqrt(PI)/GA;
           R1=-1.0;
           for (K=1; K<13; K++) {
              VA=K+0.5;
              GA = GAMMA(VA);
              VB=-K+V+0.5;
              GB = GAMMA(VB);
              R1=-R1/pow(0.5*X,2);
              S += R1*GA/GB;
		   }
           S0=SA*S;
           U=fabs(V);
           N=int(U);
           U0=U-N;
           for (L=0; L<2; L++) {
              VT=U0+L;
              R=1.0;
              BIV=1.0;
              for (K=1; K<17; K++) {
                 R=-0.125*R*(4.0*VT*VT-pow(2.0*K-1.0,2))/(K*X);
                 BIV += R;
                 if (fabs(R/BIV) < 1.0e-12) goto e30;
	      }
e30:          if (L == 0)  BIV0=BIV;
	   }
           BF0=BIV0;
           BF1=BIV;
           for (K=2; K<=N; K++) {
              BF=-2.0*(K-1.0+U0)/X*BF1+BF0;
              BF0=BF1;
              BF1=BF;
           }
           if (N == 0) BIV=BIV0;
           if (N > 1)  BIV=BF;
           *SLV=exp(X)/sqrt(2.0*PI*X)*BIV+S0;
        }
}


void main() {

        double SLV, V, VMIN, VMAX, DV, X;

        printf("\n Please enter vmin, vmax, dv and x: ");
        scanf(" %lf %lf %lf %lf", &VMIN,&VMAX,&DV,&X);

        printf("\n");
        printf("   v              Lv(x)     \n");
        printf(" ---------------------------\n");
        
        V=VMIN;                                  

e10:    STVLV(V,X,&SLV);
        
        printf(" %5f %18.8f\n", V, SLV);
	V += DV;
	if (V<=VMAX) goto e10;

	printf(" ---------------------------\n\n");

}

// end of file mstvlv.cpp