/**************************************************************
!*      Purpose:  This program computes Struve function Hv(x) *
!*                for an arbitrary order using subroutine     *
!*                STVHV                                       * 
!*      Input :   v  --- Order of Hv(x)  ( -8.0 to 12.5 )     *
!*                x  --- Argument of Hv(x) ( x > 0 )          *
!*      Output:   HV --- Hv(x)                                *
!*      Example:  x = 10.0                                    *
!*                  v           Hv(x)                         *
!*                -----------------------                     *
!*                  .5     .46402212D+00                      * 
!*                 1.5     .14452322D+01                      *                      
!*                 2.5     .31234632D+01                      *
!*                 3.5     .53730255D+01                      *
!*                 4.5     .72083122D+01                      *
!*                 5.5     .76851132D+01                      *
!* ---------------------------------------------------------- *
!* 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 STVHV(double V, double X, double *HV) {
/*      =====================================================
!       Purpose: Compute Struve function Hv(x) with an
!                arbitrary order v
!       Input :  v  --- Order of Hv(x)  ( -8.0 ó v ó 12.5 )
!                x  --- Argument of Hv(x) ( x ò 0 )
!       Output:  HV --- Hv(x)
!       Routine called: GAMMA to compute the gamma function
!       ===================================================== */
        double GA,GB,PI,PU0,PU1,QU0,QU1,R1,R2,S,S0,SA,T0,T1,U,U0,V0,VA,VB,VT;
		double BF,BF0,BF1,BY0,BY1,BYV,SR;
        int K, L, N;

		PI=3.141592653589793;
        if (X == 0.0) {
           if (V > -1.0 || int(V)-V == 0.5)
              *HV=0.0;
           else if (V < -1.0)
		      *HV=pow(-1, int(0.5-V)-1)*1.0e+300;
           else if (V == -1.0)
              *HV=2.0/PI;
           return;
        }
        if (X <= 20.0) {
           V0=V+1.5;
           GA = GAMMA(V0);
           S=2.0/(sqrt(PI)*GA);
           R1=1.0;
           for (K=1; K<101; K++) {
              VA=K+1.5;
              GA = GAMMA(VA);
              VB=V+K+1.5;
              GB = GAMMA(VB);
              R1=-R1*pow(0.5*X,2);
              R2=R1/(GA*GB);
              S=S+R2;
              if (fabs(R2) < fabs(S)*1.0e-12) goto e15;
		   }
e15:       *HV=pow(0.5*X,V+1.0)*S;
        }
        else {
           SA=pow(0.5*X, V-1.0)/PI;
           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=S+R1*GA/GB;
           }
           S0=SA*S;
           U=fabs(V);
           N=int(U);
           U0=U-N;
           for (L=0; L<2; L++) {
              VT=4.0*pow(U0+L,2);
              R1=1.0;
              PU1=1.0;
              for (K=1; K<13; K++) {
                 R1=-0.0078125*R1*(VT-pow(4.0*K-3.0,2))*(VT-pow(4.0*K-1.0,2))/((2.0*K-1.0)*K*X*X);
                 PU1=PU1+R1;
			  }
              QU1=1.0;
              R2=1.0;
              for (K=1; K<13; K++) {
                 R2=-0.0078125*R2*(VT-pow(4.0*K-1.0,2))*(VT-pow(4.0*K+1.0,2))/((2.0*K+1.0)*K*X*X);
                 QU1=QU1+R2;
			  }
              QU1=0.125*(VT-1.0)/X*QU1;
              if (L == 0) {
                 PU0=PU1;
                 QU0=QU1;
              }
		   }
           T0=X-(0.5*U0+0.25)*PI;
           T1=X-(0.5*U0+0.75)*PI;
           SR=sqrt(2.0/(PI*X));
           BY0=SR*(PU0*sin(T0)+QU0*cos(T0));
           BY1=SR*(PU1*sin(T1)+QU1*cos(T1));
           BF0=BY0;
           BF1=BY1;
           for (K=2; K<=N; K++) {
              BF=2.0*(K-1.0+U0)/X*BF1-BF0;
              BF0=BF1;
              BF1=BF;
           }
           if (N == 0) BYV=BY0;
           if (N == 1) BYV=BY1;
           if (N > 1)  BYV=BF;
           *HV=BYV+S0;
        }
}

       
void main() {

        double HV, 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         Hv(x)  \n");
        printf(" ---------------------\n");
        
        V=VMIN;                                  

e10:    STVHV(V,X,&HV);
        
        printf(" %5.1f      %f\n", V, HV);
	V += DV;
	if (V<=VMAX) goto e10;

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

}

// end of file mstvhv.cpp