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