'************************************************************** '* 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". * '* * '* 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 Lv(x)" PRINT " -------------------------" V = VMIN 10 CALL STVLV(V, X, SLV) PRINT " "; V; " "; SLV V = V + DV IF (V <= VMAX) THEN GOTO 10 PRINT " -------------------------" END '******************************************* '* 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 SUB STVLV (V, X, 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 ' ====================================================== PI = 3.141592653589793# IF (X = 0.0) THEN IF (V > -1.0 OR INT(V) - V = 0.5) THEN SLV = 0.0 ELSEIF (V < -1.0) THEN SLV = (-1) ^ (INT(0.5 - V) - 1) * 1.0D+15 ELSEIF (V = -1.0) THEN SLV = 2.0 / PI END IF EXIT SUB END IF IF (X <= 40.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 / S) < 1.0D-12) THEN GOTO 15 NEXT K 15 SLV = (0.5 * X) ^ (V + 1.0) * S ELSE SA = -1.0 / PI * (0.5 * X) ^ (V - 1.0) 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 = U0 + L R = 1.0 BIV = 1.0 FOR K = 1 TO 16 R = -0.125 * R * (4.0 * VT * VT - (2.0 * K - 1.0) ^ 2) / (K * X) BIV = BIV + R IF (ABS(R / BIV) < 1.0D-12) THEN GOTO 30 NEXT K 30 IF (L = 0) THEN BIV0 = BIV NEXT L BF0 = BIV0 BF1 = BIV 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 BIV = BIV0 IF (N > 1) THEN BIV = BF SLV = EXP(X) / SQR(2.0 * PI * X) * BIV + S0 END IF END SUB 'end of file mstvlv.bas