'********************************************************************
'*      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