DECLARE SUB CVQM (M%, Q#, A0#) DECLARE SUB CVQL (KD%, M%, Q#, A0#) DECLARE SUB CV0 (KD%, M%, Q#, A0#) DECLARE SUB CVF (KD%, M%, Q#, A#, MJ%, F#) DECLARE SUB REFINE (KD%, M%, Q#, A#, IFLAG%) DECLARE SUB CVA2 (KD%, M%, Q#, A#) DECLARE SUB FCOEF (KD%, M%, Q#, A#, FC#()) DECLARE SUB MTU0 (KF%, M%, Q#, X#, CSF#, CSD#) '*********************************************************************** '* Calculate Mathieu Functions and their First Derivatives * '* ------------------------------------------------------------------- * '* EXPLANATION: * '* * '* Purpose: This program computes Mathieu functions cem(x,q), * '* sem(x,q) and their derivatives using subroutine * '* MTU0 ( q = 0 ) * '* Input : KF --- Function code * '* KF=1 for computing cem(x,q) and cem'(x,q) * '* KF=2 for computing sem(x,q) and sem'(x,q) * '* m --- Order of Mathieu functions * '* q --- Parameter of Mathieu functions * '* x --- Argument of Mathieu functions (in degrees) * '* Output: CSF --- cem(x,q) or sem(x,q) * '* CSD --- cem'x,q) or sem'x,q) * '* Example: x = 40 * '* m q cem(x,q) cem'(x,q) sem(x,q) sem'(x,q) * '* -------------------------------------------------------- * '* 0 5.0 .3025683 .9470247 * '* 1 5.0 .7669652 1.2873097 .2988052 .9606824 * '* 2 5.0 .9102723 -.3463855 .7549264 1.4743128 * '* 5 5.0 -.9810931 -.6328576 .1694850 -4.8676455 * '* 0 25.0 .0515371 .3823737 * '* 1 25.0 .2074402 1.2646301 .0515365 .3823777 * '* 2 25.0 -.5297051 -2.4292679 .2074275 1.2646996 * '* 5 25.0 .7507159 -3.9047012 1.1881232 .3258081 * '* -------------------------------------------------------- * '* * '* ------------------------------------------------------------------- * '* SAMPLE RUNS: * '* * '* Please enter KF, m, q and x (in degrees): 1,5,25,40 * '* KF = 1, m = 5, q = 25.0, x = 40.0 * '* * '* x(degs) cem(x,q) cem'(x,q) * '* ------------------------------------- * '* 40.0 0.750715922 -3.904701223 * '* * '* Please enter KF, m, q and x (in degrees): 2,5,25,40 * '* KF = 2, m = 5, q = 25.0, x = 40.0 * '* * '* x(degs) sem(x,q) sem'(x,q) * '* ------------------------------------- * '* 40.0 1.188123243 0.325808111 * '* * '* ------------------------------------------------------------------- * '* REFERENCE: "Fortran Routines for Computation of Special Functions, * '* jin.ece.uiuc.edu/routines/routines.html". * '* * '* Quick Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '*********************************************************************** ' PROGRAM MMTU0 DEFDBL A-H, O-Z DEFINT I-N ' KF,M: Integer ' CSF,CSD,Q,X : Double F\$ = "##.#########" CLS PRINT INPUT " Please enter KF, m, q and x (in degrees): ", KF, M, Q, X PRINT " KF = "; KF; ", M = "; M; ", Q = "; Q; ", X = "; X PRINT IF KF = 1 THEN PRINT " x(degs.) cem(x,q) cem'(x,q)" IF KF = 2 THEN PRINT " x(degs.) sem(x,q) sem'(x,q)" PRINT " --------------------------------------" CALL MTU0(KF, M, Q, X, CSF, CSD) PRINT " "; X; " "; PRINT USING F\$; CSF; PRINT " "; PRINT USING F\$; CSD PRINT END 'of main program SUB CV0 (KD, M, Q, A0) ' ====================================================== ' Purpose: Compute the initial characteristic value of ' Mathieu functions for m = 12 or q = 300 or ' Q = M * M ' Input : m --- Order of Mathieu functions ' q --- Parameter of Mathieu functions ' Output: A0 --- Characteristic value ' Routines called: ' (1) CVQM for computing initial characteristic ' value for q = 3*m ' (2) CVQL for computing initial characteristic ' value for q = m*m ' ====================================================== Q2 = Q * Q IF M = 0 THEN IF Q <= 1# THEN A0 = (((.0036392# * Q2 - .0125868#) * Q2 + .0546875#) * Q2 - .5) * Q2 ELSEIF Q <= 10# THEN A0 = ((3.999267E-03 * Q - 9.638957E-02) * Q - .88297) * Q + .5542818# ELSE CALL CVQL(KD, M, Q, A0) END IF ELSEIF M = 1 THEN IF (Q <= 1#) AND (KD = 2) THEN A0 = (((-.000651 * Q - .015625) * Q - .125) * Q + 1#) * Q + 1# ELSEIF (Q <= 1#) AND (KD = 3) THEN A0 = (((-.000651 * Q + .015625) * Q - .125) * Q - 1#) * Q + 1# ELSEIF (Q <= 10#) AND (KD = 2) THEN A0 = (((-4.94603E-04 * Q + .0192917) * Q - .3089229#) * Q + 1.33372) * Q + .811752 ELSEIF (Q <= 10#) AND (KD = 3) THEN A0 = ((1.971096E-03 * Q - 5.482465E-02) * Q - 1.152218) * Q + 1.10427 ELSE CALL CVQL(KD, M, Q, A0) END IF ELSEIF M = 2 THEN IF (Q <= 1#) AND (KD = 1) THEN A0 = (((-.0036391# * Q2 + .0125888#) * Q2 - .0551939#) * Q2 + .416667) * Q2 + 4# ELSEIF (Q <= 1#) AND (KD = 4) THEN A0 = (.0003617# * Q2 - .0833333#) * Q2 + 4# ELSEIF (Q <= 15#) AND (KD = 1) THEN A0 = (((3.200972E-04 * Q - 8.667445E-03) * Q - 1.829032E-04) * Q + .9919999#) * Q + 3.3290504# ELSEIF (Q <= 10#) AND (KD = 4) THEN A0 = ((2.38446E-03 * Q - .08725329#) * Q - 4.732542E-03) * Q + 4.00909 ELSE CALL CVQL(KD, M, Q, A0) END IF ELSEIF M = 3 THEN IF (Q <= 1#) AND (KD = 2) THEN A0 = ((.0006348 * Q + .015625) * Q + .0625) * Q2 + 9# ELSEIF (Q <= 1#) AND (KD = 3) THEN A0 = ((.0006348 * Q - .015625) * Q + .0625) * Q2 + 9# ELSEIF (Q <= 20#) AND (KD = 2) THEN A0 = (((3.035731E-04 * Q - 1.453021E-02) * Q + .19069602#) * Q - .1039356#) * Q + 8.944927399999999# ELSEIF (Q <= 15#) AND (KD = 3) THEN A0 = ((9.369364E-05 * Q - .03569325#) * Q + .2689874#) * Q + 8.771735 ELSE CALL CVQL(KD, M, Q, A0) END IF ELSEIF M = 4 THEN IF (Q <= 1#) AND (KD = 1) THEN A0 = ((-.0000021 * Q2 + .0005012) * Q2 + .0333333#) * Q2 + 16# ELSEIF (Q <= 1#) AND (KD = 4) THEN A0 = ((.0000037 * Q2 - .0003669) * Q2 + .0333333#) * Q2 + 16# ELSEIF (Q <= 25#) AND (KD = 1) THEN A0 = (((1.076676E-04 * Q - .0079684875#) * Q + .17344854#) * Q - .5924058#) * Q + 16.620847# ELSEIF (Q <= 20#) AND (KD = 4) THEN A0 = ((-7.08719E-04 * Q + .0038216144#) * Q + .1907493#) * Q + 15.744 ELSE CALL CVQL(KD, M, Q, A0) END IF ELSEIF M = 5 THEN IF (Q <= 1#) AND (KD = 2) THEN A0 = ((.0000068 * Q + .0000142) * Q2 + .0208333#) * Q2 + 25# ELSEIF (Q <= 1#) AND (KD = 3) THEN A0 = ((-.0000068 * Q + .0000142) * Q2 + .0208333#) * Q2 + 25# ELSEIF (Q <= 35#) AND (KD = 2) THEN A0 = (((2.238231E-05 * Q - 2.983416E-03) * Q + .10706975#) * Q - .600205) * Q + 25.93515 ELSEIF (Q <= 25#) AND (KD = 3) THEN A0 = ((-7.425364E-04 * Q + .0218225) * Q + .0416399) * Q + 24.897 ELSE CALL CVQL(KD, M, Q, A0) END IF ELSEIF (M = 6) THEN IF (Q <= 1#) THEN A0 = (.0000004 * Q2 + .0142857#) * Q2 + 36# ELSEIF (Q <= 40#) AND (KD = 1) THEN A0 = (((-1.66846E-05 * Q + 4.80263E-04) * Q + .0253998) * Q - .181233) * Q + 36.423 ELSEIF (Q <= 35#) AND (KD = 4) THEN A0 = ((-4.57146E-04 * Q + .0216609) * Q - 2.349616E-02) * Q + 35.99251 ELSE CALL CVQL(KD, M, Q, A0) END IF ELSEIF (M = 7) THEN IF (Q <= 10#) THEN CALL CVQM(M, Q, A0) ELSEIF (Q <= 50#) AND (KD = 2) THEN A0 = (((-1.411114E-05 * Q + 9.730514E-04) * Q - 3.097887E-03) * Q + 3.533597E-02) * Q + 49.0547 ELSEIF (Q <= 40#) AND (KD = 3) THEN A0 = ((-3.043872E-04 * Q + .0205511) * Q - .0916292) * Q + 49.19035 ELSE CALL CVQL(KD, M, Q, A0) END IF ELSEIF M >= 8 THEN IF Q <= 3# * M THEN CALL CVQM(M, Q, A0) ELSEIF Q > M * M THEN CALL CVQL(KD, M, Q, A0) ELSE IF (M = 8) AND (KD = 1) THEN A0 = (((8.634308E-06 * Q - 2.100289E-03) * Q + .169072) * Q - 4.64336) * Q + 109.4211 ELSEIF (M = 8) AND (KD = 4) THEN A0 = ((-6.7842E-05 * Q + .0022057) * Q + .48296) * Q + 56.59 ELSEIF (M = 9) AND (KD = 2) THEN A0 = (((2.906435E-06 * Q - 1.019893E-03) * Q + .1101965#) * Q - 3.821851) * Q + 127.6098 ELSEIF (M = 9) AND (KD = 3) THEN A0 = ((-9.577289E-05 * Q + .01043839#) * Q + .06588934#) * Q + 78.0198 ELSEIF (M = 10) AND (KD = 1) THEN A0 = (((5.44927E-07 * Q - 3.926119E-04) * Q + .0612099#) * Q - 2.600805) * Q + 138.1923 ELSEIF (M = 10) AND (KD = 4) THEN A0 = ((-7.660143E-05 * Q + .01132506#) * Q - 9.746022999999999D-02) * Q + 99.29494 ELSEIF (M = 11) AND (KD = 2) THEN A0 = (((-5.67615E-07 * Q + 7.152722E-06) * Q + .01920291#) * Q - 1.081583) * Q + 140.88 ELSEIF (M = 11) AND (KD = 3) THEN A0 = ((-6.310551E-05 * Q + .0119247#) * Q - .2681195#) * Q + 123.667 ELSEIF (M = 12) AND (KD = 1) THEN A0 = (((-2.38351E-07 * Q - 2.90139E-05) * Q + .02023088#) * Q - 1.289) * Q + 171.2723 ELSEIF (M = 12) AND (KD = 4) THEN A0 = (((3.08902E-07 * Q - 1.577869E-04) * Q + .0247911#) * Q - 1.05454) * Q + 161.471 END IF END IF END IF END SUB SUB CVA2 (KD, M, Q, A) ' ======================================================= ' Purpose: Calculate a specific characteristic value of ' Mathieu functions ' Input : m --- Order of Mathieu functions ' q --- Parameter of Mathieu functions ' KD --- Case code ' KD=1 for cem(x,q) ( m = 0,2,4,...) ' KD=2 for cem(x,q) ( m = 1,3,5,...) ' KD=3 for sem(x,q) ( m = 1,3,5,...) ' KD=4 for sem(x,q) ( m = 2,4,6,...) ' Output: A --- Characteristic value ' Routines called: ' (1) REFINE for finding accurate characteristic ' values using an iteration method ' (2) CV0 for finding initial characteristic ' values using polynomial approximation ' (3) CVQM for computing initial characteristic ' values for q ó 3*m ' (3) CVQL for computing initial characteristic ' values for q ò m*m ' ======================================================= ' Labels: 5, 16 IF (M <= 12) OR (Q <= 3# * M) OR (Q > M * M) THEN CALL CV0(KD, M, Q, A) IF Q <> 0# THEN CALL REFINE(KD, M, Q, A, 1) ELSE NDIV = 10 DELTA = (M - 3#) * M / NDIV IF Q - 3# * M <= M * M - Q THEN 5 NN = INT((Q - 3# * M) / DELTA) + 1 DELTA = (Q - 3# * M) / NN Q1 = 2# * M CALL CVQM(M, Q1, A1) Q2 = 3# * M CALL CVQM(M, Q2, A2) QQ = 3# * M FOR I = 1 TO NN QQ = QQ + DELTA A = (A1 * Q2 - A2 * Q1 + (A2 - A1) * QQ) / (Q2 - Q1) IFLAG = 1 IF I = NN THEN IFLAG = -1 CALL REFINE(KD, M, QQ, A, IFLAG) Q1 = Q2 Q2 = QQ A1 = A2 A2 = A NEXT I IF IFLAG = -10 THEN NDIV = NDIV * 2 DELTA = (M - 3#) * M / NDIV GOTO 5 END IF ELSE 16 NN = INT((M * M - Q) / DELTA) + 1 DELTA = (M * M - Q) / NN Q1 = M * (M - 1#) CALL CVQL(KD, M, Q1, A1) Q2 = M * M CALL CVQL(KD, M, Q2, A2) QQ = M * M FOR I = 1 TO NN QQ = QQ - DELTA A = (A1 * Q2 - A2 * Q1 + (A2 - A1) * QQ) / (Q2 - Q1) IFLAG = 1 IF I = NN THEN IFLAG = -1 CALL REFINE(KD, M, QQ, A, IFLAG) Q1 = Q2 Q2 = QQ A1 = A2 A2 = A NEXT I IF IFLAG = -10 THEN NDIV = NDIV * 2 DELTA = (M - 3#) * M / NDIV GOTO 16 END IF END IF END IF END SUB SUB CVF (KD, M, Q, A, MJ, F) ' ====================================================== ' Purpose: Compute the value of F for characteristic ' equation of Mathieu functions ' Input : m --- Order of Mathieu functions ' q --- Parameter of Mathieu functions ' A --- Characteristic value ' Output: F --- Value of F for characteristic equation ' ====================================================== B = A IC = M / 2 L = 0 L0 = 0 J0 = 2 JF = IC IF KD = 1 THEN L0 = 2 IF KD = 1 THEN J0 = 3 IF (KD = 2) OR (KD = 3) THEN L = 1 IF KD = 4 THEN JF = IC - 1 T1 = 0# FOR J = MJ TO IC + 1 STEP -1 T1 = -Q * Q / ((2# * J + L) ^ 2 - B + T1) NEXT J IF M <= 2 THEN T2 = 0# IF (KD = 1) AND (M = 0) THEN T1 = T1 + T1 IF (KD = 1) AND (M = 2) THEN T1 = -2# * Q * Q / (4# - B + T1) - 4# IF (KD = 2) AND (M = 1) THEN T1 = T1 + Q IF (KD = 3) AND (M = 1) THEN T1 = T1 - Q ELSE IF KD = 1 THEN T0 = 4# - B + 2# * Q * Q / B IF KD = 2 THEN T0 = 1# - B + Q IF KD = 3 THEN T0 = 1# - B - Q IF KD = 4 THEN T0 = 4# - B T2 = -Q * Q / T0 FOR J = J0 TO JF T2 = -Q * Q / ((2# * J - L - L0) ^ 2 - B + T2) NEXT J END IF F = (2# * IC + L) ^ 2 + T1 + T2 - B END SUB SUB CVQL (KD, M, Q, A0) ' ======================================================== ' Purpose: Compute the characteristic value of Mathieu ' functions for q ò 3m ' Input : m --- Order of Mathieu functions ' q --- Parameter of Mathieu functions ' Output: A0 --- Initial characteristic value ' ======================================================== IF (KD = 1) OR (KD = 2) THEN W = 2# * M + 1# IF (KD = 3) OR (KD = 4) THEN W = 2# * M - 1# W2 = W * W W3 = W * W2 W4 = W2 * W2 W6 = W2 * W4 D1 = 5# + 34# / W2 + 9# / W4 D2 = (33# + 410# / W2 + 405# / W4) / W D3 = (63# + 1260# / W2 + 2943# / W4 + 486# / W6) / W2 D4 = (527# + 15617# / W2 + 69001# / W4 + 41607# / W6) / W3 C1 = 128# P2 = Q / W4 P1 = SQR(P2) CV1 = -2# * Q + 2# * W * SQR(Q) - (W2 + 1#) / 8# CV2 = (W + 3# / W) + D1 / (32# * P1) + D2 / (8# * C1 * P2) CV2 = CV2 + D3 / (64# * C1 * P1 * P2) + D4 / (16# * C1 * C1 * P2 * P2) A0 = CV1 - CV2 / (C1 * P1) END SUB SUB CVQM (M, Q, A0) ' ====================================================== ' Purpose: Compute the characteristic value of Mathieu ' functions for q ó m*m ' Input : m --- Order of Mathieu functions ' q --- Parameter of Mathieu functions ' Output: A0 --- Initial characteristic value ' ====================================================== HM1 = .5 * Q / (M * M - 1#) HM3 = .25 * HM1 * HM1 * HM1 / (M * M - 4#) HM5 = HM1 * HM3 * Q / ((M * M - 1#) * (M * M - 9#)) A0 = M * M + Q * (HM1 + (5# * M * M + 7#) * HM3 + (9# * M * M * M * M + 58# * M * M + 29#) * HM5) END SUB SUB FCOEF (KD, M, Q, A, FC()) ' ====================================================== ' Purpose: Compute expansion coefficients for Mathieu ' functions and modified Mathieu functions ' Input : m --- Order of Mathieu functions ' q --- Parameter of Mathieu functions ' KD --- Case code ' KD=1 for cem(x,q) ( m = 0,2,4,...) ' KD=2 for cem(x,q) ( m = 1,3,5,...) ' KD=3 for sem(x,q) ( m = 1,3,5,...) ' KD=4 for sem(x,q) ( m = 2,4,6,...) ' A --- Characteristic value of Mathieu ' functions for given m and q ' Output: FC(k) --- Expansion coefficients of Mathieu ' functions ( k= 1,2,...,KM ) ' FC(1),FC(2),FC(3),... correspond to ' A0,A2,A4,... for KD=1 case, A1,A3, ' A5,... for KD=2 case, B1,B3,B5,... ' for KD=3 case and B2,B4,B6,... for ' KD=4 case ' ====================================================== ' Labels: 45, 70, 85 IF Q <= 1# THEN QM = 7.5 + 56.1 * SQR(Q) - 134.7 * Q + 90.7 * SQR(Q) * Q ELSE QM = 17# + 3.1 * SQR(Q) - .126 * Q + .0037 * SQR(Q) * Q END IF KM = INT(QM + .5 * M) IF Q = 0# THEN FOR K = 1 TO KM FC(K) = 0# NEXT K IF KD = 1 THEN FC((M + 2) / 2) = 1# IF M = 0 THEN FC(1) = 1# / SQR(2#) ELSEIF KD = 4 THEN FC(M / 2) = 1# ELSE FC((M + 1) / 2) = 1# END IF RETURN END IF KB = 0 S = 0# F = 1D-100 U = 0# FC(KM) = 0# IF KD = 1 THEN FOR K = KM TO 3 STEP -1 V = U U = F F = (A - 4# * K * K) * U / Q - V IF ABS(F) < ABS(FC(K + 1)) THEN KB = K FC(1) = 1D-100 SP = 0# F3 = FC(K + 1) FC(2) = A / Q * FC(1) FC(3) = (A - 4#) * FC(2) / Q - 2# * FC(1) U = FC(2) F1 = FC(3) FOR I = 3 TO KB V = U U = F1 F1 = (A - 4# * (I - 1#) ^ 2) * U / Q - V FC(I + 1) = F1 IF I = KB THEN F2 = F1 IF I <> KB THEN SP = SP + F1 * F1 NEXT I SP = SP + 2# * FC(1) ^ 2 + FC(2) ^ 2 + FC(3) ^ 2 SS = S + SP * (F3 / F2) ^ 2 S0 = SQR(1# / SS) FOR J = 1 TO KM IF J <= KB + 1 THEN FC(J) = S0 * FC(J) * F3 / F2 ELSE FC(J) = S0 * FC(J) END IF NEXT J GOTO 85 ELSE FC(K) = F S = S + F * F END IF NEXT K FC(2) = Q * FC(3) / (A - 4# - 2# * Q * Q / A) FC(1) = Q / A * FC(2) S = S + 2# * FC(1) ^ 2 + FC(2) ^ 2 S0 = SQR(1# / S) FOR K = 1 TO KM FC(K) = S0 * FC(K) NEXT K ELSEIF KD = 2 OR KD = 3 THEN FOR K = KM TO 3 STEP -1 V = U U = F F = (A - (2# * K - 1#) ^ 2) * U / Q - V IF ABS(F) >= ABS(FC(K)) THEN FC(K - 1) = F S = S + F * F ELSE KB = K F3 = FC(K) GOTO 45 END IF NEXT K TMP = (-1) ^ KD FC(1) = Q / (A - 1# - TMP * Q) * FC(2) S = S + FC(1) * FC(1) S0 = SQR(1# / S) FOR K = 1 TO KM FC(K) = S0 * FC(K) NEXT K GOTO 85 45 FC(1) = 1D-100 TMP = (-1) ^ KD FC(2) = (A - 1# - TMP * Q) / Q * FC(1) SP = 0# U = FC(1) F1 = FC(2) FOR I = 2 TO KB - 1 V = U U = F1 F1 = (A - (2# * I - 1#) ^ 2) * U / Q - V IF I <> KB - 1 THEN FC(I + 1) = F1 SP = SP + F1 * F1 ELSE F2 = F1 END IF NEXT I SP = SP + FC(1) ^ 2 + FC(2) ^ 2 SS = S + SP * (F3 / F2) ^ 2 S0 = 1# / SQR(SS) FOR J = 1 TO KM IF J < KB THEN FC(J) = S0 * FC(J) * F3 / F2 IF J >= KB THEN FC(J) = S0 * FC(J) NEXT J ELSEIF KD = 4 THEN FOR K = KM TO 3 STEP -1 V = U U = F F = (A - 4# * K * K) * U / Q - V IF ABS(F) >= ABS(FC(K)) THEN FC(K - 1) = F S = S + F * F ELSE KB = K F3 = FC(K) GOTO 70 END IF NEXT K FC(1) = Q / (A - 4#) * FC(2) S = S + FC(1) * FC(1) S0 = SQR(1# / S) FOR K = 1 TO KM FC(K) = S0 * FC(K) NEXT K GOTO 85 70 FC(1) = 1D-100 FC(2) = (A - 4#) / Q * FC(1) SP = 0# U = FC(1) F1 = FC(2) FOR I = 2 TO KB - 1 V = U U = F1 F1 = (A - 4# * I * I) * U / Q - V IF I <> KB - 1 THEN FC(I + 1) = F1 SP = SP + F1 * F1 ELSE F2 = F1 END IF NEXT I SP = SP + FC(1) ^ 2 + FC(2) ^ 2 SS = S + SP * (F3 / F2) ^ 2 S0 = 1# / SQR(SS) FOR J = 1 TO KM IF J < KB THEN FC(J) = S0 * FC(J) * F3 / F2 IF J >= KB THEN FC(J) = S0 * FC(J) NEXT J END IF 85 IF FC(1) < 0# THEN FOR J = 1 TO KM FC(J) = -FC(J) NEXT J END IF END SUB SUB MTU0 (KF, M, Q, X, CSF, CSD) ' =============================================================== ' Purpose: Compute Mathieu functions cem(x,q) and sem(x,q) ' and their derivatives ( q = 0 ) ' Input : KF --- Function code ' KF=1 for computing cem(x,q) and cem'(x,q) ' KF=2 for computing sem(x,q) and sem'(x,q) ' m --- Order of Mathieu functions ' q --- Parameter of Mathieu functions ' x --- Argument of Mathieu functions (in degrees) ' Output: CSF --- cem(x,q) or sem(x,q) ' CSD --- cem'x,q) or sem'x,q) ' Routines called: ' (1) CVA2 for computing the characteristic values ' (2) FCOEF for computing the expansion coefficients ' =============================================================== ' Label: 15 DIM FG(251) EPS = .00000000000001# IF KF = 1 AND M = 2 * INT(M / 2) THEN KD = 1 IF KF = 1 AND M <> 2 * INT(M / 2) THEN KD = 2 IF KF = 2 AND M <> 2 * INT(M / 2) THEN KD = 3 IF KF = 2 AND M = 2 * INT(M / 2) THEN KD = 4 CALL CVA2(KD, M, Q, A) IF Q <= 1# THEN QM = 7.5 + 56.1 * SQR(Q) - 134.7 * Q + 90.7 * SQR(Q) * Q ELSE QM = 17# + 3.1 * SQR(Q) - .126 * Q + .0037 * SQR(Q) * Q END IF KM = INT(QM + .5 * M) CALL FCOEF(KD, M, Q, A, FG()) IC = INT(M / 2) + 1 RD = .0174532925199433# XR = X * RD CSF = 0# FOR K = 1 TO KM IF KD = 1 THEN CSF = CSF + FG(K) * COS((2 * K - 2) * XR) ELSEIF KD = 2 THEN CSF = CSF + FG(K) * COS((2 * K - 1) * XR) ELSEIF KD = 3 THEN CSF = CSF + FG(K) * SIN((2 * K - 1) * XR) ELSEIF KD = 4 THEN CSF = CSF + FG ^ (K) * SIN(2 * K * XR) END IF IF K >= IC AND ABS(FG(K)) < ABS(CSF) * EPS THEN GOTO 15 NEXT K 15 CSD = 0# FOR K = 1 TO KM IF KD = 1 THEN CSD = CSD - (2 * K - 2) * FG(K) * SIN((2 * K - 2) * XR) ELSEIF KD = 2 THEN CSD = CSD - (2 * K - 1) * FG(K) * SIN((2 * K - 1) * XR) ELSEIF KD = 3 THEN CSD = CSD + (2 * K - 1) * FG(K) * COS((2 * K - 1) * XR) ELSEIF KD = 4 THEN CSD = CSD + 2# * K * FG(K) * COS(2 * K * XR) END IF IF K >= IC AND ABS(FG(K)) < ABS(CSD) * EPS THEN EXIT SUB NEXT K END SUB SUB REFINE (KD, M, Q, A, IFLAG) ' ======================================================== ' Purpose: calculate the accurate characteristic value ' by the secant method ' Input : m --- Order of Mathieu functions ' q --- Parameter of Mathieu functions ' A --- Initial characteristic value ' Output: A --- Refineed characteristic value ' Routine called: CVF for computing the value of F for ' characteristic equation ' ======================================================== ' Labels: 6, 20 EPS = .00000000000001# MJ = 10 + M CA = A DELTA = 0# X0 = A CALL CVF(KD, M, Q, X0, MJ, F0) X1 = 1.002 * A CALL CVF(KD, M, Q, X1, MJ, F1) 6 FOR IT = 1 TO 100 MJ = MJ + 1 X = X1 - (X1 - X0) / (1# - F0 / F1) CALL CVF(KD, M, Q, X, MJ, F) IF (ABS(1# - X1 / X) < EPS) OR (F = 0#) THEN GOTO 20 X0 = X1 F0 = F1 X1 = X F1 = F NEXT IT 20 A = X IF DELTA > .05 THEN A = CA IF IFLAG < 0 THEN IFLAG = -10 EXIT SUB END IF IF ABS((A - CA) / CA) > .05 THEN X0 = CA DELTA = DELTA + .005 CALL CVF(KD, M, Q, X0, MJ, F0) X1 = (1# + DELTA) * CA CALL CVF(KD, M, Q, X1, MJ, F1) GOTO 6 END IF END SUB