/********************************************************************** * Calculate the Legendre Polynomials for a Complex Argument * * ------------------------------------------------------------------- * * EXPLANATION: * * * * ========================================================== * * Purpose: This program computes the Legendre polynomials * * Pn(z) and Pn'(z) for a complex argument using * * subroutine CLPN * * Input : x --- Real part of z * * y --- Imaginary part of z * * n --- Degree of Pn(z), n = 0,1,...,N * * Output: Complex CPN(n) --- Pn(z) * * Complex CPD(n) --- Pn'(z) * * * * Example: z = 3.0 + 2.0 i * * * * n Re[Pn(z)] Im[Pn(z)] Re[Pn'(z)] Im[Pn'(z)] * * ----------------------------------------------------------- * * 0 .100000D+01 .000000D+00 .000000D+00 .000000D+00 * * 1 .300000D+01 .200000D+01 .100000D+01 .000000D+00 * * 2 .700000D+01 .180000D+02 .900000D+01 .600000D+01 * * 3 -.270000D+02 .112000D+03 .360000D+02 .900000D+02 * * 4 -.539000D+03 .480000D+03 -.180000D+03 .790000D+03 * * 5 -.461700D+04 .562000D+03 -.481500D+04 .441000D+04 * * ========================================================== * * * * ------------------------------------------------------------------- * * SAMPLE RUN: * * * * Please enter Nmax, x and y (z=x+iy): 5 3 2 * * x = 3.0, y = 2.0 * * * * n Re[Pn(z)] Im[Pn(z)] Re[Pn'(z)] Im[Pn'(z)] * * ----------------------------------------------------------- * * 0 1.000000 0.000000 0.000000 0.000000 * * 1 3.000000 2.000000 1.000000 0.000000 * * 2 7.000000 18.000000 9.000000 6.000000 * * 3 -27.000000 112.000000 36.000000 90.000000 * * 4 -539.000000 480.000000 -180.000000 790.000000 * * 5 -4617.00D000 562.000000 -4815.000000 4410.000000 * * * * ------------------------------------------------------------------- * * REFERENCE: "Fortran Routines for Computation of Special Functions, * * jin.ece.uiuc.edu/routines/routines.html". * * * * C++ Release By J-P Moreau, Paris. * **********************************************************************/ #include #include #define NMAX 100 typedef double COMPLEX[2]; COMPLEX CPN[NMAX], CPD[NMAX]; int K,N; double X,Y; // 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 CLPN(int N, double X, double Y, COMPLEX *CPN, COMPLEX *CPD) { /* ========================================================= Purpose: Compute Legendre polynomials Pn(z) and their derivatives Pn'(z) for a complex argument Input : x --- Real part of z y --- Imaginary part of z n --- Degree of Pn(z), n = 0,1,2,... Output: Complex CPN(n) --- Pn(z) Complex CPD(n) --- Pn'(z) ========================================================= */ COMPLEX CP0,CP1,CPF,TMP,TMP1,Z; Z[0]=X; Z[1]=Y; CPN[0][0]=1.0; CPN[0][1]=0.0; CPN[1][0]=Z[1]; CPN[1][1]=Z[2]; CPD[0][0]=0.0; CPD[0][1]=0.0; CPD[1][0]=1.0; CPD[1][1]=0.0; CP0[0]=1.0; CP0[1]=0.0; CP1[0]=Z[0]; CP1[1]=Z[1]; for (K=2; K<=N; K++) { //CPF=(2.0*K-1.0)/K*Z*CP1-(K-1.0)/K*CP0 CMUL(Z,CP1,TMP); TMP[0]=(2.0*K-1.0)/K*TMP[0]; TMP[1]=(2.0*K-1.0)/K*TMP[1]; CPF[0]=TMP[0]-(K-1.0)/K*CP0[0]; CPF[1]=TMP[1]-(K-1.0)/K*CP0[1]; CPN[K][0]=CPF[0]; CPN[K][1]=CPF[1]; if (fabs(X) == 1.0 && Y == 0.0) { //CPD(K)=0.5*X^(K+1)*K*(K+1.0) CPD[K][0]=0.5*pow(X,K+1)*K*(K+1.0); CPD[K][1]=0.0; } else { //D(K)=K*(CP1-Z*CPF)/(1.0-Z*Z) CMUL(Z,CPF,TMP); TMP[0]=K*(CP1[0]-TMP[0]); TMP[1]=K*(CP1[1]-TMP[1]); CMUL(Z,Z,TMP1); TMP1[0]=1.0-TMP1[0]; TMP1[1]=-TMP1[1]; CDIV(TMP,TMP1,CPD[K]); } CP0[0]=CP1[0]; CP0[1]=CP1[1]; CP1[0]=CPF[0]; CP1[1]=CPF[1]; } } void main() { printf("\n Please enter Nmax, x and y (z=x+iy): "); scanf("%d %lf %lf", &N, &X, &Y); printf(" x =%5.1f, y =%5.1f\n\n", X, Y); CLPN(N,X,Y,CPN,CPD); printf(" n Re[Pn(z)] Im[Pn(z)] Re[Pn'(z)] Im[Pn'(z)]\n"); printf(" ------------------------------------------------------------\n"); for (K=0; K<=N; K++) printf("%4d%14.6f%14.6f %14.6f%14.6f\n", K, CPN[K][0],CPN[K][1],CPD[K][0],CPD[K][1]); printf("\n"); } // end of file mclpn.cpp