/************************************************************** !* 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 #include 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