'************************************************************** '* 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". * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '************************************************************** DEFDBL A-H, O-Z DEFINT I-N PRINT INPUT " Please enter vmin, vmax, dv and x: ", VMIN, VMAX, DV, X PRINT PRINT " v Hv(x) " PRINT " ---------------------------" V = VMIN 10 CALL STVHV(V, X, HV) PRINT " "; V; " "; HV V = V + DV IF V <= VMAX THEN GOTO 10 END SUB STVHV (V, X, 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 ' ===================================================== PI = 3.141592653589793# IF (X = 0.0) THEN IF (V > -1.0 OR INT(V) - V = 0.5) THEN HV = 0.0 ELSEIF (V < -1.0) THEN HV = (-1) ^ (INT(0.5 - V) - 1) * 1.0E15 ELSEIF (V = -1.0) THEN HV = 2.0 / PI END IF EXIT SUB END IF IF (X <= 20.0) THEN V0 = V + 1.5 GA = Gamma(V0) S = 2.0 / (SQR(PI) * GA) R1 = 1.0 FOR K = 1 TO 100 VA = K + 1.5 GA = Gamma(VA) VB = V + K + 1.5 GB = Gamma(VB) R1 = -R1 * (0.5 * X) ^ 2 R2 = R1 / (GA * GB) S = S + R2 IF (ABS(R2) < ABS(S) * 1.0D-12) THEN GOTO 15 NEXT K 15 HV = (0.5 * X) ^ (V + 1.0) * S ELSE SA = (0.5 * X) ^ (V - 1.0) / PI V0 = V + 0.5 GA = Gamma(V0) S = SQR(PI) / GA R1 = 1.0 FOR K = 1 TO 12 VA = K + 0.5 GA = Gamma(VA) VB = -K + V + 0.5 GB = Gamma(VB) R1 = R1 / (0.5 * X) ^ 2 S = S + R1 * GA / GB NEXT K S0 = SA * S U = ABS(V) N = INT(U) U0 = U - N FOR L = 0 TO 1 VT = 4.0 * (U0 + L) ^ 2 R1 = 1.0 PU1 = 1.0 FOR K = 1 TO 12 R1 = -0.0078125 * R1 * (VT - (4.0 * K - 3.0) ^ 2) * (VT - (4.0 * K - 1.0) ^ 2) / ((2.0 * K - 1.0) * K * X * X) PU1 = PU1 + R1 NEXT K QU1 = 1.0 R2 = 1.0 FOR K = 1 TO 12 R2 = -0.0078125 * R2 * (VT - (4.0 * K - 1.0) ^ 2) * (VT - (4.0 * K + 1.0) ^ 2) / ((2.0 * K + 1.0) * K * X * X) QU1 = QU1 + R2 NEXT K QU1 = 0.125 * (VT - 1.0) / X * QU1 IF (L = 0) THEN PU0 = PU1 QU0 = QU1 END IF NEXT L T0 = X - (0.5 * U0 + 0.25) * PI T1 = X - (0.5 * U0 + 0.75) * PI SR = SQR(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 TO N BF = 2.0 * (K - 1.0 + U0) / X * BF1 - BF0 BF0 = BF1 BF1 = BF NEXT K IF (N = 0) THEN BYV = BY0 IF (N = 1) THEN BYV = BY1 IF (N > 1) THEN BYV = BF HV = BYV + S0 END IF END SUB '******************************************* '* FUNCTION GAMMA(X) * '* --------------------------------------- * '* Returns the value of Gamma(x) in double * '* precision as EXP(LN(GAMMA(X))) for X>0. * '* (input: x output: y) * '******************************************* FUNCTION Gamma (x) DIM cof(6) cof(1) = 76.18009173#: cof(2) = -86.50532033# cof(3) = 24.01409822#: cof(4) = -1.231739516# cof(5) = 0.00120858003#: cof(6) = -0.00000536382# stp = 2.50662827465# half = 0.5#: one = 1#: fpf = 5.5# xx = x - one tmp = xx + fpf tmp = (xx + half) * LOG(tmp) - tmp ser = one FOR j = 1 TO 6 xx = xx + one ser = ser + cof(j) / xx NEXT j Gamma = EXP(tmp + LOG(stp * ser)) END FUNCTION 'end of file mstvhv.bas