/******************************************************************** !* 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" * !* * !* C++ Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !*******************************************************************/ #include #include double CV1[200], CV2[200], CVE[200], CVS[200]; int J, KD, MMAX; double Q; void CVA1(int KD, int M, double Q, double *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,... ! ============================================================ */ //Labels: 55, 60, 70 double G[200], H[200]; double D[500], E[500], F[500]; double EPS, S,T,T1, X1,XA,XB; int I,IC,ICM, J,K,K1, NM,NM1; EPS=1e-14; ICM=int(M/2)+1; if (KD == 4) ICM=(M/2); if (Q == 0.0) if (KD == 1) for (IC=1; IC<=ICM; IC++) CV[IC]=4.0*(IC-1.0)*(IC-1.0); else if (KD != 4) for (IC=1; IC<=ICM; IC++) CV[IC]=(2.0*IC-1.0)*(2.0*IC-1.0); else for (IC=1; IC<=ICM; IC++) CV[IC]=4.0*IC*IC; else { NM=int(10+1.5*M+0.5*Q); E[1]=0.0; F[1]=0.0; if (KD == 1) { D[1]=0.0; for (I=2; I<=NM; I++) { D[I]=4.0*(I-1)*(I-1); E[I]=Q; F[I]=Q*Q; } E[2]=sqrt(2.0)*Q; F[2]=2.0*Q*Q; } else if (KD != 4) { D[1]=1.0+pow(-1,KD)*Q; for (I=2; I<=NM; I++) { D[I]=(2.0*I-1.0)*(2.0*I-1.0); E[I]=Q; F[I]=Q*Q; } } else { D[1]=4.0; for (I=2; I<=NM; I++) { D[I]=4.0*I*I; E[I]=Q; F[I]=Q*Q; } } XA=D[NM]+fabs(E[NM]); XB=D[NM]-fabs(E[NM]); NM1=NM-1; for (I=1; I<=NM1; I++) { T=fabs(E[I])+fabs(E[I+1]); T1=D[I]+T; if (XA < T1) XA=T1; T1=D[I]-T; if (T1 < XB) XB=T1; } for (I=1; I<=ICM; I++) { G[I]=XA; H[I]=XB; } for (K=1; K<=ICM; K++) { for (K1=K; K1<=ICM; K1++) { if (G[K1] < G[K]) { G[K]=G[K1]; goto e55; } } e55: if (K != 1 && H[K] < H[K-1]) H[K]=H[K-1]; e60: X1=(G[K]+H[K])/2.0; CV[K]=X1; if (fabs((G[K]-H[K])/X1) < EPS) goto e70; J=0; S=1.0; for (I=1; I<=NM; I++) { if (S == 0.0) S += 1e-30; T=F[I]/S; S=D[I]-T-X1; if (S < 0.0) J++; } if (J < K) H[K]=X1; else { G[K]=X1; if (J >= ICM) G[ICM]=X1; else { if (H[J+1] < X1) H[J+1]=X1; if (X1 < G[J]) G[J]=X1; } } goto e60; e70: CV[K]=X1; } } } void main() { printf("\n Please enter Mmax,q: "); scanf("%d %lf", &MMAX, &Q); printf("\n"); CVA1(1,MMAX,Q,CV1); CVA1(2,MMAX,Q,CV2); for (J=1; J<=MMAX/2+1; J++) { CVE[2*J-1]=CV1[J]; CVE[2*J]=CV2[J]; } CVA1(3,MMAX,Q,CV1); CVA1(4,MMAX,Q,CV2); for (J=1; J<=MMAX/2+1; J++) { CVS[2*J]=CV1[J]; CVS[2*J+1]=CV2[J]; } printf(" Characteristic values of Mathieu functions\n"); printf("\n"); printf(" m a b \n"); printf("------------------------------------\n"); for (J=0; J<=MMAX; J++) { if (J == 0) printf(" %2d %f\n", J,CVE[J+1]); if (J != 0) printf(" %2d %f %f\n", J,CVE[J+1],CVS[J+1]); } } // end of file mcva1a.cpp