/********************************************************************** * Calculate the Legendre polynomials for a complex argument * * ------------------------------------------------------------------- * * EXPLANATION: * * * * ========================================================== * * Purpose: This program computes the Legendre polynomials * * Qn(z) and Qn'(z) for a complex argument using * * subroutine CLQN * * Input : x --- Real part of z * * y --- Imaginary part of z * * n --- Degree of Qn(z), n = 0,1,... * * Output: CQN(n) --- Qn(z) * * CQD(n) --- Qn'(z) * * Examples: * * * * z = 0.5 + 0.5 i * * n Re[Qn(z)] Im[Qn(z)] Re[Qn'(z)] Im[Qn'(z)] * * ----------------------------------------------------------- * * 0 .402359D+00 .553574D+00 .800000D+00 .400000D+00 * * 1 -.107561D+01 .477967D+00 .602359D+00 .115357D+01 * * 2 -.136636D+01 -.725018D+00 -.242682D+01 .183390D+01 * * 3 .182619D+00 -.206146D+01 -.622944D+01 -.247151D+01 * * 4 .298834D+01 -.110022D+01 -.114849D+01 -.125963D+02 * * 5 .353361D+01 .334847D+01 .206656D+02 -.123735D+02 * * * * z = 3.0 + 2.0 i * * n Re[Qn(z)] Im[Qn(z)] Re[Qn'(z)] Im[Qn'(z)] * * ----------------------------------------------------------- * * 0 .229073D+00 -.160875D+00 -.250000D-01 .750000D-01 * * 1 .896860D-02 -.244805D-01 .407268D-02 .141247D-01 * * 2 -.736230D-03 -.281865D-02 .190581D-02 .155860D-02 * * 3 -.264727D-03 -.227023D-03 .391535D-03 .314880D-04 * * 4 -.430648D-04 -.443187D-05 .527190D-04 -.305592D-04 * * 5 -.481362D-05 .265297D-05 .395108D-05 -.839883D-05 * * ========================================================== * * * * ------------------------------------------------------------------- * * SAMPLE RUN: * * * * Please enter Nmax, x and y (z=x+iy): 5 3 2 * * x = 3.0, y = 2.0 * * * * n Re[Qn(z)] Im[Qn(z)] Re[Qn'(z)] Im[Qn'(z)] * * ----------------------------------------------------------- * * 0 0.22907268 -0.16087528 -0.02500000 0.07500000 * * 1 0.00896860 -0.02448047 0.00407268 0.01412472 * * 2 -0.00073623 -0.00281865 0.00190581 0.00155860 * * 3 -0.00026473 -0.00022702 0.00039153 0.00003149 * * 4 -0.00004306 -0.00000443 0.00005272 -0.00003056 * * 5 -0.00000481 0.00000265 0.00000395 -0.00000840 * * * * ------------------------------------------------------------------- * * 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 100 typedef double COMPLEX[2]; COMPLEX CQN[NMAX], CQD[NMAX]; int K,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; } void CLOG(COMPLEX ZA, COMPLEX ZB, int *IERR) { /***BEGIN PROLOGUE CLOG DOUBLE PRECISION COMPLEX LOGARITHM B=CLOG(A) IERR=0,NORMAL RETURN IERR=1, Z=CMPLX(0.0,0.0) ***ROUTINES CALLED ZABS ***END PROLOGUE CLOG */ // Labels: e10,e20,e30,e40,e50,e60,e70 double AR,AI,BR,BI,ZM, DTHETA, DPI, DHPI; DPI = 3.141592653589793238462643383E+0000; DHPI= 1.570796326794896619231321696E+0000; *IERR=0; AR=ZA[0]; AI=ZA[1]; 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; BI = DHPI; BR = log(fabs(AI)); if (AI < 0.0) BI = -BI; goto e70; e20: if (AR > 0.0) goto e30; BR = log(fabs(AR)); BI = DPI; goto e70; e30: BR = log(AR); BI = 0.0; goto e70; e40: if (AR < 0.0) DTHETA += DPI; e50: ZM = CABS(ZA); BR = log(ZM); BI = DTHETA; goto e70; e60: *IERR=1; e70: ZB[0]=BR; ZB[1]=BI; } //CLOG() // 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 CLQN(int N, double X, double Y, COMPLEX *CQN, COMPLEX *CQD) { /* ===================================================== Purpose: Compute the Legendre functions Qn(z) and their derivatives Qn'(z) for a complex argument Input : x --- Real part of z y --- Imaginary part of z n --- Degree of Qn(z), n = 0,1,2,... Output: CQN(n) --- Qn(z) CQD(n) --- Qn'(z) ===================================================== */ COMPLEX CONE,CQ0,CQ1,CQF0,CQF1,CQF2,Z,TMP,TMP1,TMP2; int Err,K,KM,LS; CONE[0]=1.0; CONE[1]=0.0; Z[0]=X; Z[1]=Y; if (Z[0]==1.0 && Z[1]==0.0) { for (K=0; K<=N; K++) { CQN[K][0]=1.0e+200; CQN[K][1]=0.0; CQD[K][0]=1.0e+200; CQD[K][1]=0.0; } return; } LS=1; if (CABS(Z) > 1.0) LS=-1; //CQ0=0.5*CLOG(LS*(1.0+Z)/(1.0-Z)) TMP[0]=LS*(CONE[0]+Z[0]); TMP[1]=LS*(CONE[1]+Z[1]); TMP1[0]=CONE[0]-Z[0]; TMP1[1]=CONE[1]-Z[1]; CDIV(TMP,TMP1,TMP2); CLOG(TMP2,CQ0,&Err); CQ0[0]=CQ0[0]*0.5; CQ0[1]=CQ0[1]*0.5; //CQ1=Z*CQ0-1.0 CMUL(Z,CQ0,CQ1); CQ1[0]=CQ1[0]-CONE[0]; CQN[0][0]=CQ0[0]; CQN[0][1]=CQ0[1]; CQN[1][0]=CQ1[0]; CQN[1][1]=CQ1[1]; if (CABS(Z) < 1.0001) { CQF0[0]=CQ0[0]; CQF0[1]=CQ0[1]; CQF1[0]=CQ1[0]; CQF1[1]=CQ1[1]; for (K=2; K<=N; K++) { //CQF2=((2.0*K-1.0)*Z*CQF1-(K-1.0)*CQF0)/K CMUL(Z,CQF1,TMP); TMP[0]=TMP[0]*(2.0*K-1.0); TMP[1]=TMP[1]*(2.0*K-1.0); TMP1[0]=(K-1.0)*CQF0[0]; TMP1[1]=(K-1.0)*CQF0[1]; CQF2[0]=(TMP[0]-TMP1[0])/K; CQF2[1]=(TMP[1]-TMP1[1])/K; CQN[K][0]=CQF2[0]; CQN[K][1]=CQF2[1]; CQF0[0]=CQF1[0]; CQF0[1]=CQF1[1]; CQF1[0]=CQF2[0]; CQF1[1]=CQF2[1]; } } else { if (CABS(Z) > 1.1) KM=40+N; else { TMP[0]=Z[0]-1.0; TMP[1]=Z[1]; KM=(int) ((40+N)*floor(-1.0-1.8*log(CABS(TMP)))); } CQF2[0]=0.0; CQF2[1]=0.0; CQF1[0]=1.0; CQF1[1]=0.0; for (K=KM; K>=0; K--) { //CQF0=((2*K+3.0)*Z*CQF1-(K+2.0)*CQF2)/(K+1.0) CMUL(Z,CQF1,TMP); TMP[0]=TMP[0]*(2.0*K+3.0); TMP[1]=TMP[1]*(2.0*K+3.0); TMP1[0]=(K+2.0)*CQF2[0]; TMP1[1]=(K+2.0)*CQF2[1]; CQF0[0]=(TMP[0]-TMP1[0])/(1.0+K); CQF0[1]=(TMP[1]-TMP1[1])/(1.0+K); if (K <= N) { CQN[K][0]=CQF0[0]; CQN[K][1]=CQF0[1]; } CQF2[0]=CQF1[0]; CQF2[1]=CQF1[1]; CQF1[0]=CQF0[0]; CQF1[1]=CQF0[1]; } for (K=0; K<=N; K++) { //CQN(K)=CQN(K)*CQ0/CQF0 CMUL(CQN[K],CQ0,TMP); CDIV(TMP,CQF0,CQN[K]); } } //CQD(0)=(CQN(1)-Z*CQN(0))/(Z*Z-1.0) CMUL(Z,CQN[0],TMP); TMP[0]=CQN[1][0]-TMP[0]; TMP[1]=CQN[1][1]-TMP[1]; CMUL(Z,Z,TMP1); TMP1[0]=TMP1[0]-1.0; CDIV(TMP,TMP1,CQD[0]); for (K=1; K<=N; K++) { //CQD(K)=(K*Z*CQN(K)-K*CQN(K-1))/(Z*Z-1.0) CMUL(Z,CQN[K],TMP); TMP[0]=K*TMP[0]; TMP[1]=K*TMP[1]; TMP2[0]=K*CQN[K-1][0]; TMP2[1]=K*CQN[K-1][1]; TMP[0]=TMP[0]-TMP2[0]; TMP[1]=TMP[1]-TMP2[1]; CDIV(TMP,TMP1,CQD[K]); } } 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); CLQN(N,X,Y,CQN,CQD); printf(" n Re[Qn(z)] Im[Qn(z)] Re[Qn'(z)] Im[Qn'(z)]\n"); printf(" -----------------------------------------------------------\n"); for (K=0; K<=N; K++) printf("%4d%14.8f%14.8f%14.8f%14.8f\n", K,CQN[K][0],CQN[K][1],CQD[K][0],CQD[K][1]); printf("\n"); } // end of file mclqn.cpp