'******************************************************************** '* Purpose: This program computes a sequence of characteristic * '* values of Mathieu functions using subroutine CVA1 * '* 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: CV(I) --- Characteristic values; I = 1,2,3,... * '* For KD=1, CV(1), CV(2), CV(3),..., correspond to * '* For KD=2, CV(1), CV(2), CV(3),..., correspond to * '* the characteristic values of cem for m = 1,3,5,.. * '* For KD=3, CV(1), CV(2), CV(3),..., correspond to * '* the characteristic values of sem for m = 1,3,5,.. * '* For KD=4, CV(1), CV(2), CV(3),..., correspond to * '* the characteristic values of sem for m = 0,2,4,.. * '* * '* Example: Mmax = 12, q = 25.00 * '* * '* Characteristic values of Mathieu functions * '* * '* m a b * '* ------------------------------------------ * '* 0 -40.256779547 * '* 1 -21.314899691 -40.256778985 * '* 2 -3.522164727 -21.314860622 * '* 3 12.964079444 -3.520941527 * '* 4 27.805240581 12.986489953 * '* 5 40.050190986 28.062765899 * '* 6 48.975786716 41.801071292 * '* 7 57.534689001 55.002957151 * '* 8 69.524065166 69.057988351 * '* 9 85.076999882 85.023356505 * '* 10 103.230204804 103.225680042 * '* 11 123.643012376 123.642713667 * '* 12 146.207690643 146.207674647 * '* ---------------------------------------------------------------- * '* 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 DIM CV1(200) DIM CV2(200) DIM CVE(200) DIM CVS1(200) PRINT INPUT " Please enter Mmax,q: ", MMAX, Q PRINT CALL CVA1(1, MMAX, Q, CV1()) CALL CVA1(2, MMAX, Q, CV2()) FOR J = 1 TO MMAX / 2 + 1 CVE(2 * J - 1) = CV1(J) CVE(2 * J) = CV2(J) NEXT J CALL CVA1(3, MMAX, Q, CV1()) CALL CVA1(4, MMAX, Q, CV2()) FOR J = 1 TO MMAX / 2 + 1 CVS1(2 * J) = CV1(J) CVS1(2 * J + 1) = CV2(J) NEXT J PRINT " Characteristic values of Mathieu functions" PRINT PRINT " m a b " PRINT "------------------------------------------" FOR J = 0 TO MMAX IF (J = 0) THEN PRINT " "; J; " "; CVE(J + 1) END IF IF (J <> 0) THEN PRINT " "; J; " "; CVE(J + 1); " "; CVS1(J + 1) END IF NEXT J END SUB CVA1 (KD, M, Q, CV()) ' ============================================================ ' Purpose: Compute a sequence of characteristic values of ' Mathieu functions ' Input : M --- Maximum 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: CV(I) --- Characteristic values; I = 1,2,3,... ' For KD=1, CV(1), CV(2), CV(3),..., correspond to ' the characteristic values of cem for m = 0,2,4,... ' For KD=2, CV(1), CV(2), CV(3),..., correspond to ' the characteristic values of cem for m = 1,3,5,... ' For KD=3, CV(1), CV(2), CV(3),..., correspond to ' the characteristic values of sem for m = 1,3,5,... ' For KD=4, CV(1), CV(2), CV(3),..., correspond to ' the characteristic values of sem for m = 0,2,4,... ' ============================================================ DIM G(200), H(200), D(500), E(500), F(500) EPS = 1.0D-14 ICM = INT(M / 2) + 1 IF (KD = 4) THEN ICM = M / 2 IF (Q = 0.0) THEN IF (KD.EQ.1) THEN FOR IC = 1 TO ICM CV(IC) = 4.0 * (IC - 1.0) ^ 2 NEXT IC ELSEIF (KD <> 4) THEN FOR IC = 1 TO ICM CV(IC) = (2.0 * IC - 1.0) ^ 2 NEXT IC ELSE FOR IC = 1 TO ICM CV(IC) = 4.0D0 * IC * IC NEXT IC END IF ELSE NM = INT(10 + 1.5 * M + 0.5 * Q) E(1) = 0.0 F(1) = 0.0 IF (KD = 1) THEN D(1) = 0.0 FOR I = 2 TO NM D(I) = 4.0 * (I - 1.0) ^ 2 E(I) = Q F(I) = Q * Q NEXT I E(2) = SQR(2.0) * Q F(2) = 2.0 * Q * Q ELSEIF (KD <> 4) THEN D(1) = 1.0 + (-1) ^ KD * Q FOR I = 2 TO NM D(I) = (2.0 * I - 1.0) ^ 2 E(I) = Q F(I) = Q * Q NEXT I ELSE D(1) = 4.0 FOR I = 2 TO NM D(I) = 4.0 * I * I E(I) = Q F(I) = Q * Q NEXT I END IF XA = D(NM) + ABS(E(NM)) XB = D(NM) - ABS(E(NM)) NM1 = NM - 1 FOR I = 1 TO NM1 T = ABS(E(I)) + ABS(E(I + 1)) T1 = D(I) + T IF (XA < T1) THEN XA = T1 T1 = D(I) - T IF (T1 < XB) THEN XB = T1 NEXT I FOR I = 1 TO ICM G(I) = XA H(I) = XB NEXT I FOR K = 1 TO ICM FOR K1 = K TO ICM IF (G(K1) < G(K)) THEN G(K) = G(K1) GOTO 55 END IF NEXT K1 55 IF (K <> 1 AND H(K) < H(K - 1)) THEN H(K) = H(K - 1) 60 X1 = (G(K) + H(K)) / 2.0 CV(K) = X1 IF (ABS((G(K) - H(K)) / X1) < EPS) THEN GOTO 70 J = 0 S = 1.0 FOR I = 1 TO NM IF (S = 0.0) THEN S = S + 1.0D-30 T = F(I) / S S = D(I) - T - X1 IF (S < 0.0) THEN J = J + 1 NEXT I IF (J < K) THEN H(K) = X1 ELSE G(K) = X1 IF (J >= ICM) THEN G(ICM) = X1 ELSE IF (H(J + 1) < X1) THEN H(J + 1) = X1 IF (X1 < G(J)) THEN G(J) = X1 END IF END IF GOTO 60 70 CV(K) = X1 NEXT K END IF END SUB ' end of file mcva1.bas