/************************************************************************* * Calculate the Asociated Legendre Functions and their First Derivatives * * for a Complex Argument * * ---------------------------------------------------------------------- * * EXPLANATION: * * * * ============================================================ * * Purpose: This program computes the associated Legendre * * functions Pmn(z) and their derivatives Pmn'(z) for * * a complex argument using subroutine CLPMN * * Input : x --- Real part of z * * y --- Imaginary part of z * * m --- Order of Pmn(z), m = 0,1,2,...,n * * n --- Degree of Pmn(z), n = 0,1,2,...,N * * Output: CPM(m,n) --- Pmn(z) * * CPD(m,n) --- Pmn'(z) * * Examples: * * n = 5, x = 0.5, y = 0.2 * * * * m Re[Pmn(z)] Im[Pmn(z)] Re[Pmn'(z)] Im[Pmn'(z)] * * ------------------------------------------------------------- * * 0 .252594D+00 -.530293D+00 -.347606D+01 -.194250D+01 * * 1 .333071D+01 .135206D+01 .117643D+02 -.144329D+02 * * 2 -.102769D+02 .125759D+02 .765713D+02 .598500D+02 * * 3 -.572879D+02 -.522744D+02 -.343414D+03 .147389D+03 * * 4 .335711D+03 -.389151D+02 -.226328D+03 -.737100D+03 * * 5 -.461125D+03 .329122D+03 .187180D+04 .160494D+02 * * * * n = 5, x = 2.5, y = 1.0 * * * * m Re[Pmn(z)] Im[Pmn(z)] Re[Pmn'(z)] Im[Pmn'(z)] * * ------------------------------------------------------------- * * 0 -.429395D+03 .900336D+03 -.350391D+02 .193594D+04 * * 1 -.216303D+04 .446358D+04 -.208935D+03 .964685D+04 * * 2 -.883477D+04 .174005D+05 -.123703D+04 .381938D+05 * * 3 -.273211D+05 .499684D+05 -.568080D+04 .112614D+06 * * 4 -.565523D+05 .938503D+05 -.167147D+05 .219713D+06 * * 5 -.584268D+05 .863328D+05 -.233002D+05 .212595D+06 * * ============================================================ * * * * ---------------------------------------------------------------------- * * SAMPLE RUN: * * * * Please enter m, n, x and y: 5 5 0.5 0.2 * * m = 5, n = 5, x = 0.500000, y = 0.200000 * * * * m n Re[Pmn(z)] Im[Pmn(z)] Re[Pmn'(z)] Im[Pmn'(z)] * * ----------------------------------------------------------------- * * 0 5 0.252594 -0.530293 -3.476063 -1.942500 * * 1 5 3.330709 1.352057 11.764333 -14.432917 * * 2 5 -10.276875 12.575850 76.571250 59.850000 * * 3 5 -57.287858 -52.274368 -343.414490 147.389427 * * 4 5 335.711250 -38.915100 -226.327500 -737.100000 * * 5 5 -461.124511 329.122362 1871.802220 16.049431 * * * * ---------------------------------------------------------------------- * * 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 25 typedef double COMPLEX[2]; typedef COMPLEX CTab[NMAX][NMAX]; CTab CPM, CPD; int J,M,N; double X,Y; //ABSOLUTE VALUE OF A COMPLEX NUMBER Z=X+I*Y double CABS(COMPLEX Z) { // Labels: e10, e20 double U,V,Q,S; U = fabs(Z[0]); V = fabs(Z[1]); S = U+V; /*------------------------------------------------------------------- S*1.0 MAKES AN UNNORMALIZED UNDERFLOW ON CDC MACHINES INTO A TRUE FLOATING ZERO -------------------------------------------------------------------*/ S = S*1.0; if (S == 0.0) goto e20; if (U > V) goto e10; Q = U/V; return (V*sqrt(1.0+Q*Q)); e10: Q = V/U; return (U*sqrt(1.0+Q*Q)); e20: return 0.0; } // Z=Z1/Z2 void CDIV(COMPLEX Z1, COMPLEX Z2, COMPLEX Z) { double D; D=Z2[0]*Z2[0]+Z2[1]*Z2[1]; if (D<1e-12) return; Z[0]=(Z1[0]*Z2[0]+Z1[1]*Z2[1])/D; Z[1]=(Z1[1]*Z2[0]-Z1[0]*Z2[1])/D; } // Z=Z1*Z2 void CMUL(COMPLEX Z1, COMPLEX Z2, COMPLEX Z) { Z[0]=Z1[0]*Z2[0] - Z1[1]*Z2[1]; Z[1]=Z1[0]*Z2[1] + Z1[1]*Z2[0]; } void CSQRT(COMPLEX ZA, COMPLEX ZB) { /***BEGIN PROLOGUE CSQRT ! ! DOUBLE PRECISION COMPLEX SQUARE ROOT, CSQRT(ZA,ZB) ! !***FUNCTIONS CALLED: CABS() !***END PROLOGUE CSQRT */ // Labels: e10,e20,e30,e40,e50,e60,e70,e80 double AR,AI,BR,BI,ZM, DTHETA, DPI, DRT; DRT=7.071067811865475244008443621e-001; DPI=3.141592653589793238462643383e+000; AR=ZA[0]; AI=ZA[1]; ZM = CABS(ZA); ZM = sqrt(ZM); if (AR == 0.0) goto e10; if (AI == 0.0) goto e20; DTHETA = atan(AI/AR); if (DTHETA <= 0.0) goto e40; if (AR < 0.0) DTHETA -= DPI; goto e50; e10: if (AI > 0.0) goto e60; if (AI < 0.0) goto e70; BR = 0.0; BI = 0.0; goto e80; e20: if (AR > 0.0) goto e30; BR = 0.0; BI = sqrt(fabs(AR)); goto e80; e30: BR = sqrt(AR); BI = 0.0; goto e80; e40: if (AR < 0.0) DTHETA += DPI; e50: DTHETA *= 0.5; BR = ZM*cos(DTHETA); BI = ZM*sin(DTHETA); goto e80; e60: BR = ZM*DRT; BI = ZM*DRT; goto e80; e70: BR = ZM*DRT; BI = -ZM*DRT; e80 : ZB[0]=BR; ZB[1]=BI; } //ZSQRT() void CLPMN(int M, int N, double X, double Y, CTab CPM, CTab CPD) { /* ================================================================ Purpose: Compute the associated Legendre functions Pmn(z) and their derivatives Pmn'(z) for a complex argument Input : x --- Real part of z y --- Imaginary part of z m --- Order of Pmn(z), m = 0,1,2,...,n n --- Degree of Pmn(z), n = 0,1,2,...,N mm --- Physical dimension of CPM and CPD Output: Complex CPM(m,n) --- Pmn(z) Complex CPD(m,n) --- Pmn'(z) ============================================================= */ COMPLEX TMP,TMP1,TMP2,Z,ZQ,ZS; int I,J,LS; Z[0]=X; Z[1]=Y; for (I=0; I<=N; I++) for (J=0; J<=M; J++) { CPM[J][I][0]=0.0; CPM[J][I][1]=0.0; CPD[J][I][0]=0.0; CPD[J][I][1]=0.0; } CPM[0][0][0]=1.0; if (fabs(X) == 1.0 && Y == 0.0) { for (I=1; I<=N; I++) { CPM[0][I][0]=pow(X,I); CPM[0][I][1]=0.0; CPD[0][I][0]=0.5*I*(I+1)*pow(X,(I+1)); CPD[0][I][1]=0.0; } for (J=1; J<=N; J++) for (I=1; I<=M; I++) if (I == 1) { CPD[I][J][0]=1e+200; CPD[I][J][1]=0.0; } else if (I == 2) { CPD[I][J][0]=-0.25*(J+2)*(J+1)*J*(J-1)*pow(X,(J+1)); CPD[I][J][1]=0.0; } return; } LS=1; if (CABS(Z) > 1.0) LS=-1; // ZQ=CSQRT(LS*(1.0-Z*Z)) CMUL(Z,Z,TMP); TMP[0]=LS*(1.0-TMP[0]); TMP[1]=-LS*TMP[1]; CSQRT(TMP,ZQ); // ZS=LS*(1.0D0-Z*Z) ZS[0]=TMP[0]; ZS[1]=TMP[1]; for (I=1; I<=M; I++) { // CPM(I,I)=-LS*(2.0*I-1.0)*ZQ*CPM(I-1,I-1) CMUL(ZQ,CPM[I-1][I-1],TMP); CPM[I][I][0]=-LS*(2.0*I-1.0)*TMP[0]; CPM[I][I][1]=-LS*(2.0*I-1.0)*TMP[1]; } for (I=0; I<=M; I++) { // CPM(I,I+1)=(2.0*I+1.0)*Z*CPM(I,I) CMUL(Z,CPM[I][I],TMP); CPM[I][I+1][0]=(2.0*I+1.0)*TMP[0]; CPM[I][I+1][1]=(2.0*I+1.0)*TMP[1]; } for (I=0; I<=M; I++) for (J=I+2; J<=N; J++) { //CPM(I,J)=((2.0*J-1.0)*Z*CPM(I,J-1)-(I+J-1.0)*CPM(I,J-2))/(J-I) CMUL(Z,CPM[I][J-1],TMP); CPM[I][J][0]=((2.0*J-1.0)*TMP[0]-(I+J-1.0)*CPM[I][J-2][0])/(J-I); CPM[I][J][1]=((2.0*J-1.0)*TMP[1]-(I+J-1.0)*CPM[I][J-2][1])/(J-I); } CPD[0][0][0]=0.0; CPD[0][0][1]=0.0; for (J=1; J<=N; J++) { //CPD(0,J)=LS*J*(CPM(0,J-1)-Z*CPM(0,J))/ZS CMUL(Z,CPM[0][J],TMP); TMP[0]=LS*J*(CPM[0][J-1][0]-TMP[0]); TMP[1]=LS*J*(CPM[0][J-1][1]-TMP[1]); CDIV(TMP,ZS,CPD[0][J]); } for (I=1; I<=M; I++) for (J=I; J<=N; J++) { //CPD(I,J)=LS*I*Z*CPM(I,J)/ZS+(J+I)*(J-I+1.0)*CPM(I-1,J)/ZQ CMUL(Z,CPM[I][J],TMP); TMP[0]=TMP[0]*LS*I; TMP[1]=TMP[1]*LS*I; CDIV(TMP,ZS,TMP2); //Now TMP2=LS*I*Z*CPM(I,J)/ZS TMP[0]=TMP2[0]; TMP[1]=TMP2[1]; TMP1[0]=(J+I)*(J-I+1.0)*CPM[I-1][J][0]; TMP1[1]=(J+I)*(J-I+1.0)*CPM[I-1][J][1]; CDIV(TMP1,ZQ,TMP2); //Now TMP2=(J+I)*(J-I+1.0)*CPM(I-1,J)/ZQ CPD[I][J][0]=TMP[0]+TMP2[0]; CPD[I][J][1]=TMP[1]+TMP2[1]; } } void main() { printf("\n Please enter m, n, x and y: "); scanf("%d %d %lf %lf", &M,&N,&X,&Y); printf(" m = %d, n = %d, x = %f, y = %f\n\n", M, N, X, Y); CLPMN(M,N,X,Y,CPM,CPD); printf(" m n Re[Pmn(z)] Im[Pmn(z)] Re[Pmn'(z)] Im[Pmn'(z)]\n"); printf(" -------------------------------------------------------------------\n"); for (J=0; J<=M; J++) { printf(" %d %d %12.6f %12.6f %12.6f %12.6f\n", J, N, CPM[J][N][0], CPM[J][N][1], CPD[J][N][0], CPD[J][N][1]); } printf("\n"); } // end of file mclpmn.cpp