/************************************************************************* * Calculate the associated Legendre Functions of the Second Kind and * * their First Derivatives for a Complex Argument * * ---------------------------------------------------------------------- * * EXPLANATION: * * * * ============================================================ * * Purpose: This program computes the associated Legendre * * functions Qmn(z) and their derivatives Qmn'(z) for * * a complex argument using subroutine CLQMN * * Definition: Qmn(z)=(-1)^m*(1-z*z)^(m/2)*dm/dzm[Qn(z)] * * Q0(z)=1/2*LOG[(1+z)/(1-z)] ( for |z|<1 ) * * Qmn(z)=(z*z-1)^(m/2)*dm/dzm[Qn(z)] * * Q0(z)=1/2*LOG[(z+1)/(z-1)] ( for |z|>1 ) * * Input : x --- Real part of z * * y --- Imaginary part of z * * m --- Order of Qmn(z) ( m = 0,1,2,... ) * * n --- Degree of Qmn(z) ( n = 0,1,2,... ) * * Output: COMPLEX CQM(m,n) --- Qmn(z) * * COMPLEX CQD(m,n) --- Qmn'(z) * * Examples: * * n = 5, x = 0.5, y = 0.2 * * * * m Re[Qmn(z)] Im[Qmn(z)] Re[Qmn'(z)] Im[Qmn'(z)] * * ------------------------------------------------------------- * * 0 .987156D+00 .354345D+00 .324023D+01 -.447297D+01 * * 1 -.240328D+01 .436861D+01 .281158D+02 .171437D+02 * * 2 -.245853D+02 -.138072D+02 -.106283D+03 .913792D+02 * * 3 .102723D+03 -.651233D+02 -.362578D+03 -.429802D+03 * * 4 .155510D+03 .357712D+03 .196975D+04 -.287414D+02 * * 5 -.167357D+04 -.680954D+03 -.193093D+04 -.925757D+03 * * * * n = 5, x = 2.5, y = 1.0 * * * * m Re[Qmn(z)] Im[Qmn(z)] Re[Qmn'(z)] Im[Qmn'(z)] * * ------------------------------------------------------------- * * 0 -.274023D-04 -.227141D-04 .809834D-04 .210884D-04 * * 1 .165620D-03 .136108D-03 -.489095D-03 -.124400D-03 * * 2 -.118481D-02 -.948832D-03 .349090D-02 .825057D-03 * * 3 .982179D-02 .753264D-02 -.288271D-01 -.596384D-02 * * 4 -.927915D-01 -.669521D-01 .270840D+00 .451376D-01 * * 5 .985601D+00 .656737D+00 -.285567D+01 -.332533D+00 * * ============================================================ * * * * ---------------------------------------------------------------------- * * SAMPLE RUN: * * * * Please enter m, n, x and y: 5 5 0.5 0.2 * * n = 5, m = 5, x = 0.500000, y = 0.200000 * * * * m n Re[Qmn(z)] Im[Qmn(z)] Re[Qmn'(z)] Im[Qmn'(z)] * * ----------------------------------------------------------------- * * 0 5 0.987156 0.354345 3.240234 -4.472973 * * 1 5 -2.403283 4.368612 28.115787 17.143746 * * 2 5 -24.585261 -13.807231 -106.282948 91.379180 * * 3 5 102.723445 -65.123308 -362.578436 -429.802309 * * 4 5 155.510087 357.711653 1969.747136 -28.741415 * * 5 5 -1673.568876 -680.953539 -1930.929669 -925.757384 * * * * ---------------------------------------------------------------------- * * 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 CQM, CQD; 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; } 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 CABS ***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 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; } //CSQRT() void CLQMN(int M, int N, double X, double Y, CTab CQM, CTab CQD) { /* ======================================================== Purpose: Compute the associated Legendre functions of the second kind, Qmn(z) and Qmn'(z), for a complex argument Input : x --- Real part of z y --- Imaginary part of z m --- Order of Qmn(z) ( m = 0,1,2,... ) n --- Degree of Qmn(z) ( n = 0,1,2,... ) Output: COMPLEX CQM(m,n) --- Qmn(z) COMPLEX CQD(m,n) --- Qmn'(z) ======================================================== complex numbers: */ double CQ0[2],CQ1[2],CQ10[2],CQF[2],CQF0[2],CQF1[2]; double CQF2[2],TMP[2],TMP1[2],TMP2[2],Z[2],ZQ[2],ZS[2]; int Err,I,K,KM,LS; double XC; Z[0]=X; Z[1]=Y; if (fabs(X) == 1.0 && Y == 0.0) { for (I=0; I<=M; I++) for (J=0; J<=N; J++) { CQM[I][J][0]=1e+200; CQM[I][J][1]=0.0; CQD[I][J][0]=1e+200; CQD[I][J][1]=0.0; } return; } XC=CABS(Z); if (Z[1] == 0.0 || XC < 1.0) LS=1; if (XC > 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.0-Z*Z) ZS[0]=TMP[0]; ZS[1]=TMP[1]; // CQ0=0.5*CLOG(LS*(1.0+Z)/(1.0-Z)) TMP[0]=LS*(1.0+Z[0]); TMP[1]=LS*Z[1]; TMP1[0]=1.0-Z[0]; TMP1[1]=-Z[1]; CDIV(TMP,TMP1,TMP2); CLOG(TMP2,CQ0,&Err); CQ0[0]=0.5*CQ0[0]; CQ0[1]=0.5*CQ0[1]; if (XC < 1.0001) { CQM[0][0][0]=CQ0[0]; CQM[0][0][1]=CQ0[1]; // CQM(0,1)=Z*CQ0-1.0 CMUL(Z,CQ0,TMP); CQM[0][1][0]=TMP[0]-1.0; CQM[0][1][1]=TMP[1]; // CQM(1,0)=-1.0/ZQ TMP[0]=-1.0; TMP[1]=0.0; CDIV(TMP,ZQ,CQM[1][0]); // CQM(1,1)=-ZQ*(CQ0+Z/(1.0-Z*Z)) CMUL(Z,Z,TMP); TMP[0]=1.0-TMP[0]; TMP[1]=-TMP[1]; CDIV(Z,TMP,TMP2); TMP[0]=CQ0[0]+TMP2[0]; TMP[1]=CQ0[1]+TMP2[1]; CMUL(ZQ,TMP,CQM[1][1]); CQM[1][1][0]=-CQM[1][1][0]; CQM[1][1][1]=-CQM[1][1][1]; for (I=0; I<2; I++) for (J=2; J<=N; J++) { // CQM(I,J)=((2.0*J-1.0)*Z*CQM(I,J-1)-(J+I-1.0)*CQM(I,J-2))/(J-I) CMUL(Z,CQM[I][J-1],TMP); TMP[0]=(2.0*J-1.0)*TMP[0]; TMP[1]=(2.0*J-1.0)*TMP[1]; TMP1[0]=(J+I-1.0)*CQM[I][J-2][0]; TMP1[1]=(J+I-1.0)*CQM[I][J-2][1]; CQM[I][J][0]=(TMP[0]-TMP1[0])/(J-I); CQM[I][J][1]=(TMP[1]-TMP1[1])/(J-I); } for (J=0; J<=N; J++) for (I=2; I<=M; I++) { // CQM(I,J)=-2.0*(I-1.0)*Z/ZQ*CQM(I-1,J)-LS*(J+I-1.0)*(J-I+2.0)*CQM(I-2,J) CDIV(Z,ZQ,TMP); TMP[0]=-2.0*(I-1.0)*TMP[0]; TMP[1]=-2.0*(I-1.0)*TMP[1]; CMUL(TMP,CQM[I-1][J],TMP2); //Now TMP2=-2.0*(I-1.0)*Z/ZQ*CQM(I-1,J) TMP1[0]=LS*(J+I-1.0)*(J-I+2.0)*CQM[I-2][J][0]; TMP1[1]=LS*(J+I-1.0)*(J-I+2.0)*CQM[I-2][J][1]; //Now TMP1=LS*(J+I-1.0)*(J-I+2.0)*CQM(I-2,J) CQM[I][J][0]=TMP2[0]-TMP1[0]; CQM[I][J][1]=TMP2[1]-TMP1[1]; } } else { if (XC > 1.1) KM=40+M+N; else KM=(int) ((40+M+N)*floor(-1.0-1.8*log(XC-1.0))); 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]=(2*K+3.0)*TMP[0]; TMP[1]=(2*K+3.0)*TMP[1]; TMP1[0]=(K+2.0)*CQF2[0]; TMP1[1]=(K+2.0)*CQF2[1]; CQF0[0]=(TMP[0]-TMP1[0])/(K+1.0); CQF0[1]=(TMP[1]-TMP1[1])/(K+1.0); if (K <= N) { CQM[0][K][0]=CQF0[0]; CQM[0][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++) { // CQM(0,K)=CQ0*CQM(0,K)/CQF0 CMUL(CQ0,CQM[0][K],TMP); CDIV(TMP,CQF0,CQM[0][K]); } 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+1.0)*CQF2)/(K+2.0) CMUL(Z,CQF1,TMP); TMP[0]=(2*K+3.0)*TMP[0]; TMP[1]=(2*K+3.0)*TMP[1]; TMP1[0]=(K+1.0)*CQF2[0]; TMP1[1]=(K+1.0)*CQF2[1]; CQF0[0]=(TMP[0]-TMP1[0])/(K+2.0); CQF0[1]=(TMP[1]-TMP1[1])/(K+2.0); if (K <= N) { CQM[1][K][0]=CQF0[0]; CQM[1][K][1]=CQF0[1]; } CQF2[0]=CQF1[0]; CQF2[1]=CQF1[1]; CQF1[0]=CQF0[0]; CQF1[1]=CQF0[1]; } // CQ10=-1.0/ZQ TMP[0]=-1.0; TMP[1]=0.0; CDIV(TMP,ZQ,CQ10); for (K=0; K<=N; K++) { // CQM(1,K)=CQ10*CQM(1,K)/CQF0 CMUL(CQ10,CQM[1][K],TMP); CDIV(TMP,CQF0,CQM[1][K]); } for (J=0; J<=N; J++) { CQ0[0]=CQM[0][J][0]; CQ0[1]=CQM[0][J][1]; CQ1[0]=CQM[1][J][0]; CQ1[1]=CQM[1][J][1]; for (I=0; I