/********************************************************************** * 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.000000, x = 40.000000 * * * * 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.000000, x = 40.000000 * * * * 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". * * * * C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * **********************************************************************/ #include #include #define NMAX 252 //index zero not used int KF,M; double CSF,CSD,Q,X; //Headers of functions used below void CVA2(int KD, int M, double Q, double *A); void FCOEF(int KD, int M, double Q, double A, double *FC); void CV0(int KD, int M, double Q, double *A0); void REFINE(int KD, int M, double Q, double *A, int IFLAG); void CVQM(int M, double Q, double *A0); void CVQL(int KD, int M, double Q, double *A0); void CVF(int KD, int M, double Q, double A, int MJ, double *F); double Sqr(double x) { return x*x; } //main subroutine void MTU0(int KF, int M, double Q, double X, double *CSF, double *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 =============================================================== */ // Labels: e15, e25 double FG[NMAX]; double A,EPS,QM,RD,XR; int IC,K,KD,KM; EPS=1e-14; if (KF == 1 && M == 2*(M / 2)) KD=1; if (KF == 1 && M != 2*(M / 2)) KD=2; if (KF == 2 && M != 2*(M / 2)) KD=3; if (KF == 2 && M == 2*(M / 2)) KD=4; CVA2(KD,M,Q,&A); if (Q <= 1.0) QM=7.5+56.1*sqrt(Q)-134.7*Q+90.7*sqrt(Q)*Q; else QM=17.0+3.1*sqrt(Q)-0.126*Q+0.0037*sqrt(Q)*Q; KM=(int) floor(QM+0.5*M); FCOEF(KD,M,Q,A,FG); IC=int(M/2)+1; RD=1.74532925199433e-2; XR=X*RD; *CSF=0.0; for (K=1; K<=KM; K++) { if (KD == 1) *CSF = *CSF + FG[K]*cos((2*K-2)*XR); else if (KD == 2) *CSF = *CSF + FG[K]*cos((2*K-1)*XR); else if (KD == 3) *CSF = *CSF + FG[K]*sin((2*K-1)*XR); else if (KD == 4) *CSF = *CSF + FG[K]*sin(2*K*XR); if (K >= IC && fabs(FG[K]) < fabs(*CSF)*EPS) goto e15; } e15: *CSD=0.0; for (K=1; K<=KM; K++) { if (KD == 1) *CSD = *CSD - (2*K-2)*FG[K]*sin((2*K-2)*XR); else if (KD == 2) *CSD = *CSD - (2*K-1)*FG[K]*sin((2*K-1)*XR); else if (KD == 3) *CSD = *CSD + (2*K-1)*FG[K]*cos((2*K-1)*XR); else if (KD == 4) *CSD = *CSD + 2.0*K*FG[K]*cos(2*K*XR); if (K >= IC && fabs(FG[K]) < fabs(*CSD)*EPS) return; } } //MTU0 void FCOEF(int KD, int M, double Q, double A, double *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: e45, e70, e85 double F,F1,F2,F3,QM,S,S0,SP,SS,TMP,U,V; int I,J,K,KB,KM; if (Q <= 1.0) QM=7.5+56.1*sqrt(Q)-134.7*Q+90.7*sqrt(Q)*Q; else QM=17.0+3.1*sqrt(Q)-0.126*Q+0.0037*sqrt(Q)*Q; KM=(int) floor(QM+0.5*M); if (Q == 0.0) { for (K=1; K<=KM; K++) FC[K]=0.0; if (KD == 1) { FC[(M+2) / 2] = 1.0; if (M == 0) FC[1]=1.0/sqrt(2.0); } else if (KD == 4) FC[M/2]=1.0; else FC[(M+1)/2]=1.0; return; } KB=0; S=0.0; F=1.0e-100; U=0.0; FC[KM]=0.0; if (KD == 1) { for (K=KM; K>2; K--) { V=U; U=F; F=(A-4.0*K*K)*U/Q-V; if (fabs(F) < fabs(FC[K+1])) { KB=K; FC[1]=1.0e-100; SP=0.0; F3=FC[K+1]; FC[2]=A/Q*FC[1]; FC[3]=(A-4.0)*FC[2]/Q-2.0*FC[1]; U=FC[2]; F1=FC[3]; for (I=3 ; I<=KB; I++) { V=U; U=F1; F1=(A-4.0*Sqr(I-1.0))*U/Q-V; FC[I+1]=F1; if (I == KB) F2=F1; if (I != KB) SP += F1*F1; } SP += 2.0*Sqr(FC[1])+Sqr(FC[2])+Sqr(FC[3]); SS=S+SP*Sqr(F3/F2); S0=sqrt(1.0/SS); for (J=1; J<=KM; J++) if (J <= KB+1) FC[J]=S0*FC[J]*F3/F2; else FC[J]=S0*FC[J]; goto e85; } else { FC[K]=F; S += F*F; } } //of K loop FC[2]=Q*FC[3]/(A-4.0-2.0*Q*Q/A); FC[1]=Q/A*FC[2]; S += 2.0*Sqr(FC[1])+Sqr(FC[2]); S0=sqrt(1.0/S); for (K=1; K<=KM; K++) FC[K]=S0*FC[K]; } else if (KD == 2 || KD == 3) { for (K=KM; K>2; K--) { V=U; U=F; F=(A-Sqr(2.0*K-1.0))*U/Q-V; if (fabs(F) >= fabs(FC[K])) { FC[K-1]=F; S += F*F; } else { KB=K; F3=FC[K]; goto e45; } } //TMP=(-1)^KD if (KD==2*(KD/2)) TMP=1.0; else TMP=-1.0; FC[1]=Q/(A-1.0-TMP*Q)*FC[2]; S=S+FC[1]*FC[1]; S0=sqrt(1.0/S); for (K=1; K<=KM; K++) FC[K]=S0*FC[K]; goto e85; e45: FC[1]=1.0e-100; //TMP=(-1)^KD if (KD==2*(KD/2)) TMP=1.0; else TMP=-1.0; FC[2]=(A-1.0-TMP*Q)/Q*FC[1]; SP=0.0; U=FC[1]; F1=FC[2]; for (I=2; I= KB) FC[J]=S0*FC[J]; } } else if (KD == 4) { for (K=KM; K>2; K--) { V=U; U=F; F=(A-4.0*K*K)*U/Q-V; if (fabs(F) >= fabs(FC[K])) { FC[K-1]=F; S += F*F; } else { KB=K; F3=FC[K]; goto e70; } } FC[1]=Q/(A-4.0)*FC[2]; S += FC[1]*FC[1]; S0=sqrt(1.0/S); for (K=1; K<=KM; K++) FC[K]=S0*FC[K]; goto e85; e70: FC[1]=1.0e-100; FC[2]=(A-4.0)/Q*FC[1]; SP=0.0; U=FC[1]; F1=FC[2]; for (I=2; I= KB) FC[J]=S0*FC[J]; } } e85: if (FC[1] < 0.0) for (J=1; J<=KM; J++) FC[J]=-FC[J]; } void CVA2(int KD, int M, double Q, double *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: e5, e15 int I,IFLAG,NDIV,NN; double A1,A2,DELTA,Q1,Q2,QQ; if (M <= 12 || Q <= 3.0*M || Q > M*M) { CV0(KD,M,Q,A); if (Q != 0.0) REFINE(KD,M,Q,A,1); } else { NDIV=10; DELTA=(M-3.0)*M/NDIV; if (Q-3.0*M <= M*M-Q) { e5: NN=(int) floor((Q-3.0*M)/DELTA)+1; DELTA=(Q-3.0*M)/NN; Q1=2.0*M; CVQM(M,Q1,&A1); Q2=3.0*M; CVQM(M,Q2,&A2); QQ=3.0*M; for (I=1; I<=NN; I++) { QQ=QQ+DELTA; *A=(A1*Q2-A2*Q1+(A2-A1)*QQ)/(Q2-Q1); IFLAG=1; if (I == NN) IFLAG=-1; REFINE(KD,M,QQ,A,IFLAG); Q1=Q2; Q2=QQ; A1=A2; A2=*A; } if (IFLAG == -10) { NDIV=NDIV*2; DELTA=(M-3.0)*M/NDIV; goto e5; } } else { e15: NN=(int) floor((M*M-Q)/DELTA)+1; DELTA=(M*M-Q)/NN; Q1=M*(M-1.0); CVQL(KD,M,Q1,&A1); Q2=M*M; CVQL(KD,M,Q2,&A2); QQ=M*M; for (I=1; I<=NN; I++) { QQ=QQ-DELTA; *A=(A1*Q2-A2*Q1+(A2-A1)*QQ)/(Q2-Q1); IFLAG=1; if (I == NN) IFLAG=-1; REFINE(KD,M,QQ,A,IFLAG); Q1=Q2; Q2=QQ; A1=A2; A2=*A; } if (IFLAG == -10) { NDIV *= 2; DELTA=(M-3.0)*M/NDIV; goto e15; } } } } //CVA2() void REFINE(int KD, int M, double Q, double *A, int 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: e5, e15 double CA,DELTA,F,F0,F1,EPS,X,X0,X1; int IT,MJ; EPS=1E-14; MJ=10+M; CA=*A; DELTA=0.0; X0=*A; CVF(KD,M,Q,X0,MJ,&F0); X1=1.002*(*A); CVF(KD,M,Q,X1,MJ,&F1); e5: for (IT=1; IT<=100; IT++) { MJ=MJ+1; X=X1-(X1-X0)/(1.0-F0/F1); CVF(KD,M,Q,X,MJ,&F); if (fabs(1.0-X1/X) < EPS || F == 0.0) goto e15; X0=X1; F0=F1; X1=X; F1=F; } e15: *A=X; if (DELTA > 0.05) { *A=CA; if (IFLAG < 0) IFLAG=-10; return; } if (fabs((*A-CA)/CA) > 0.05) { X0=CA; DELTA=DELTA+0.005; CVF(KD,M,Q,X0,MJ,&F0); X1=(1.0+DELTA)*CA; CVF(KD,M,Q,X1,MJ,&F1); goto e5; } } void CVF(int KD, int M, double Q, double A, int MJ, double *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 ====================================================== */ double B,T0,T1,T2; int IC,J,J0,JF,L,L0; B=A; IC=M/2; L=0; L0=0; J0=2; JF=IC; if (KD == 1) L0=2; if (KD == 1) J0=3; if (KD == 2 || KD == 3) L=1; if (KD == 4) JF=IC-1; T1=0.0; for (J=MJ; J>IC; J--) T1=-Q*Q/(Sqr(2.0*J+L)-B+T1); if (M <= 2) { T2=0.0; if (KD == 1 && M == 0) T1=T1+T1; if (KD == 1 && M == 2) T1=-2.0*Q*Q/(4.0-B+T1)-4.0; if (KD == 2 && M == 1) T1=T1+Q; if (KD == 3 && M == 1) T1=T1-Q; } else { if (KD == 1) T0=4.0-B+2.0*Q*Q/B; if (KD == 2) T0=1.0-B+Q; if (KD == 3) T0=1.0-B-Q; if (KD == 4) T0=4.0-B; T2=-Q*Q/T0; for (J=J0; J<=JF; J++) T2=-Q*Q/(Sqr(2.0*J-L-L0)-B+T2); } *F=Sqr(2.0*IC+L) + T1 + T2 - B; } //CVF() void CV0(int KD, int M, double Q, double *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 ==================================================== */ double Q2; Q2=Q*Q; if (M == 0) { if (Q <= 1.0) *A0=(((0.0036392*Q2-0.0125868)*Q2+0.0546875)*Q2-0.5)*Q2; else if (Q <= 10.0) *A0=((3.999267E-3*Q-9.638957E-2)*Q-0.88297)*Q +0.5542818; else CVQL(KD,M,Q,A0); } else if (M == 1) { if (Q <= 1.0 && KD == 2) *A0=(((-6.51e-4*Q-0.015625)*Q-0.125)*Q+1.0)*Q+1.0; else if (Q <= 1.0 && KD == 3) *A0=(((-6.51e-4*Q+0.015625)*Q-0.125)*Q-1.0)*Q+1.0; else if (Q <= 10.0 && KD == 2) *A0=(((-4.94603e-4*Q+1.92917E-2)*Q-0.3089229)*Q+1.33372)*Q+0.811752; else if (Q <= 10.0 && KD == 3) *A0=((1.971096e-3*Q-5.482465E-2)*Q-1.152218)*Q+1.10427; else CVQL(KD,M,Q,A0); } else if (M == 2) { if (Q <= 1.0 && KD == 1) *A0=(((-0.0036391*Q2+0.0125888)*Q2-0.0551939)*Q2+0.416667)*Q2+4.0; else if (Q <= 1.0 && KD == 4) *A0=(0.0003617*Q2-0.0833333)*Q2+4.0; else if (Q <= 15.0 && KD == 1) *A0=(((3.200972e-4*Q-8.667445E-3)*Q-1.829032e-4)*Q+0.9919999)*Q+3.3290504; else if (Q <= 10.0 && KD == 4) *A0=((2.38446e-3*Q-0.08725329)*Q-4.732542e-3)*Q+4.00909; else CVQL(KD,M,Q,A0); } else if (M == 3) { if (Q <= 1.0 && KD == 2) *A0=((6.348e-4*Q+0.015625)*Q+0.0625)*Q2+9.0; else if (Q <= 1.0 && KD == 3) *A0=((6.348e-4*Q-0.015625)*Q+0.0625)*Q2+9.0; else if (Q <= 20.0 && KD == 2) *A0=(((3.035731e-4*Q-1.453021e-2)*Q+0.19069602)*Q-0.1039356)*Q+8.9449274; else if (Q <= 15.0 && KD == 3) *A0=((9.369364e-5*Q-0.03569325)*Q+0.2689874)*Q+8.771735; else CVQL(KD,M,Q,A0); } else if (M == 4) { if (Q <= 1.0 && KD == 1) *A0=((-2.1e-6*Q2+5.012e-4)*Q2+0.0333333)*Q2+16.0; else if (Q <= 1.0 && KD == 4) *A0=((3.7e-6*Q2-3.669e-4)*Q2+0.0333333)*Q2+16.0; else if (Q <= 25.0 && KD == 1) *A0=(((1.076676e-4*Q-7.9684875e-3)*Q+0.17344854)*Q-0.5924058)*Q+16.620847; else if (Q <= 20.0 && KD == 4) *A0=((-7.08719e-4*Q+3.8216144e-3)*Q+0.1907493)*Q+15.744; else CVQL(KD,M,Q,A0); } else if (M == 5) { if (Q <= 1.0 && KD == 2) *A0=((6.8e-6*Q+1.42e-5)*Q2+0.0208333)*Q2+25.0; else if (Q <= 1.0 && KD == 3) *A0=((-6.8e-6*Q+1.42e-5)*Q2+0.0208333)*Q2+25.0; else if (Q <= 35.0 && KD == 2) *A0=(((2.238231e-5*Q-2.983416e-3)*Q+0.10706975)*Q-0.600205)*Q+25.93515; else if (Q <= 25.0 && KD == 3) *A0=((-7.425364e-4*Q+2.18225e-2)*Q+4.16399e-2)*Q+24.897; else CVQL(KD,M,Q,A0); } else if (M == 6) { if (Q <= 1.0) *A0=(0.4e-6*Q2+0.0142857)*Q2+36.0; else if (Q <= 40.0 && KD == 1) *A0=(((-1.66846e-5*Q+4.80263e-4)*Q+2.53998e-2)*Q-0.181233)*Q+36.423; else if (Q <= 35.0 && KD == 4) *A0=((-4.57146e-4*Q+2.16609e-2)*Q-2.349616e-2)*Q+35.99251; else CVQL(KD,M,Q,A0); } else if (M == 7) { if (Q <= 10.0) CVQM(M,Q,A0); else if (Q <= 50.0 && KD == 2) *A0=(((-1.411114e-5*Q+9.730514e-4)*Q-3.097887e-3)*Q+3.533597e-2)*Q+49.0547; else if (Q <= 40.0 && KD == 3) *A0=((-3.043872e-4*Q+2.05511e-2)*Q-9.16292e-2)*Q+49.19035; else CVQL(KD,M,Q,A0); } else if (M >= 8) { if (Q <= 3.*M) CVQM(M,Q,A0); else if (Q > M*M) CVQL(KD,M,Q,A0); else { if (M == 8 && KD == 1) *A0=(((8.634308e-6*Q-2.100289e-3)*Q+0.169072)*Q-4.64336)*Q+109.4211; else if (M == 8 && KD == 4) *A0=((-6.7842e-5*Q+2.2057e-3)*Q+0.48296)*Q+56.59; else if (M == 9 && KD == 2) *A0=(((2.906435e-6*Q-1.019893e-3)*Q+0.1101965)*Q-3.821851)*Q+127.6098; else if (M == 9 && KD == 3) *A0=((-9.577289e-5*Q+0.01043839)*Q+0.06588934)*Q+78.0198; else if (M == 10 && KD == 1) *A0=(((5.44927e-7*Q-3.926119e-4)*Q+0.0612099)*Q-2.600805)*Q+138.1923; else if (M == 10 && KD == 4) *A0=((-7.660143e-5*Q+0.01132506)*Q-0.09746023)*Q+99.29494; else if (M == 11 && KD == 2) *A0=(((-5.67615e-7*Q+7.152722e-6)*Q+0.01920291)*Q-1.081583)*Q+140.88; else if (M == 11 && KD == 3) *A0=((-6.310551e-5*Q+0.0119247)*Q-0.2681195)*Q+123.667; else if (M == 12 && KD == 1) *A0=(((-2.38351e-7*Q-2.90139e-5)*Q+0.02023088)*Q-1.289)*Q+171.2723; else if (M == 12 && KD == 4) *A0=(((3.08902e-7*Q-1.577869e-4)*Q+0.0247911)*Q-1.05454)*Q+161.471; } } } //CV0() void CVQL(int KD, int M, double Q, double *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 ======================================================== */ double C1,CV1,CV2,D1,D2,D3,D4,P1,P2,W,W2,W3,W4,W6; if (KD == 1 || KD == 2) W=2.0*M+1.0; if (KD == 3 || KD == 4) W=2.0*M-1.0; W2=W*W; W3=W*W2; W4=W2*W2; W6=W2*W4; D1=5.0+34.0/W2+9.0/W4; D2=(33.0+410.0/W2+405.0/W4)/W; D3=(63.0+1260.0/W2+2943.0/W4+486.0/W6)/W2; D4=(527.0+15617.0/W2+69001.0/W4+41607.0/W6)/W3; C1=128.0; P2=Q/W4; P1=sqrt(P2); CV1=-2.0*Q+2.0*W*sqrt(Q)-(W2+1.0)/8.0; CV2=(W+3.0/W)+D1/(32.0*P1)+D2/(8.0*C1*P2); CV2=CV2+D3/(64.0*C1*P1*P2)+D4/(16.0*C1*C1*P2*P2); *A0=CV1-CV2/(C1*P1); } void CVQM(int M, double Q, double *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 ===================================================== */ double HM1,HM3,HM5; HM1=0.5*Q/(M*M-1.0); HM3=0.25*HM1*HM1*HM1/(M*M-4.0); HM5=HM1*HM3*Q/((M*M-1.0)*(M*M-9.0)); *A0=M*M+Q*(HM1+(5.0*M*M+7.0)*HM3+(9.0*M*M*M*M+58.0*M*M+29.0)*HM5); } void main() { printf("\n Please enter KF, m, q and x (in degrees): "); scanf("%d %d %lf %lf", &KF, &M, &Q, &X); printf(" KF = %d M = %d Q = %f X = %f\n", KF,M,Q,X); printf("\n"); if (KF == 1) printf(" x(degs.) cem(x,q) cem'(x,q)\n"); if (KF == 2) printf(" x(degs.) sem(x,q) sem'(x,q)\n"); printf(" --------------------------------------\n"); MTU0(KF,M,Q,X,&CSF,&CSD); printf(" %5.1f%16.9f%16.9f\n\n", X, CSF, CSD); } // end of file mmtu0.cpp