DECLARE SUB CLQMN (M%, N%, X#, Y#) DECLARE FUNCTION CABS# (Z#()) DECLARE SUB CDIV (Z1#(), Z2#(), Z#()) DECLARE SUB CLOG (ZA#(), ZB#(), IERR%) DECLARE SUB CMUL (Z1#(), Z2#(), Z#()) DECLARE SUB CSQRT (ZA#(), ZB#()) '************************************************************************** '* 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 = .5, y = .2 * '* * '* 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". * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '************************************************************************** 'PROGRAM MCLQMN DEFDBL A-H, O-Z DEFINT I-N OPTION BASE 0 ' Note: here complex numbers are represented by Tables(0..1) of doubles ' complex CQM(M+1,N+1,2), CQD(M+1,N+1,2) ' integer J, M, N ' double X, Y CLS PRINT INPUT " Please enter m, n, x and y: ", M, N, X, Y PRINT " m = "; M; ", n = "; N; ", x = "; X; ", y = "; Y PRINT DIM SHARED CQM(M + 1, N + 1, 2), CQD(M + 1, N + 1, 2) CALL CLQMN(M, N, X, Y) F\$ = " #####.######" PRINT " m n Re[Qmn(z)] Im[Qmn(z)] Re[Qmn'(z)] Im[Qmn'(z)]" PRINT " -------------------------------------------------------------------" FOR J = 0 TO M PRINT " "; J; " "; N; PRINT USING F\$; CQM(J, N, 0); PRINT USING F\$; CQM(J, N, 1); PRINT " "; PRINT USING F\$; CQD(J, N, 0); PRINT USING F\$; CQD(J, N, 1) NEXT J PRINT END 'of main program 'ABSOLUTE VALUE OF A COMPLEX NUMBER Z=X+I*Y FUNCTION CABS (Z()) 'double U,V,Q,S U = ABS(Z(0)) V = ABS(Z(1)) S = U + V IF S = 0# THEN GOTO 20 IF U > V THEN GOTO 10 Q = U / V CABS = V * SQR(1# + Q * Q) EXIT FUNCTION 10 Q = V / U CABS = U * SQR(1# + Q * Q) EXIT FUNCTION 20 CABS = 0# END FUNCTION ' Z=Z1/Z2 SUB CDIV (Z1(), Z2(), Z()) ' double D D = Z2(0) * Z2(0) + Z2(1) * Z2(1) IF D < .000000000001# THEN EXIT SUB Z(0) = (Z1(0) * Z2(0) + Z1(1) * Z2(1)) / D Z(1) = (Z1(1) * Z2(0) - Z1(0) * Z2(1)) / D END SUB SUB CLOG (ZA(), ZB(), 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 DPI = 4# * ATN(1#) DHPI = .5 * DPI IERR = 0 AR = ZA(0): AI = ZA(1) IF AR = 0# THEN GOTO 11 IF AI = 0# THEN GOTO 21 DTHETA = ATN(AI / AR) IF DTHETA <= 0# THEN GOTO 40 IF AR < 0# THEN DTHETA = DTHETA - DPI GOTO 50 11 IF AI = 0# THEN GOTO 60 BI = DHPI BR = LOG(ABS(AI)) IF AI < 0# THEN BI = -BI GOTO 70 21 IF AR > 0# THEN GOTO 30 BR = LOG(ABS(AR)) BI = DPI GOTO 70 30 BR = LOG(AR) BI = 0# GOTO 70 40 IF AR < 0# THEN DTHETA = DTHETA + DPI 50 ZM = CABS(ZA()) BR = LOG(ZM) BI = DTHETA GOTO 70 60 IERR = 1 70 ZB(0) = BR: ZB(1) = BI END SUB 'CLOG SUB CLQMN (M, N, X, Y) ' ======================================================== ' 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: DIM CQ0(2), CQ1(2), CQ10(2), CQF(2), CQF0(2), CQF1(2), CQF2(2) DIM TMP(2), TMP1(2), TMP2(2), TMP3(2), Z(2), ZQ(2), ZS(2) ' integer I,K,KM,LS ' double XC DIM Err1 AS INTEGER Z(0) = X: Z(1) = Y IF ABS(X) = 1# AND Y = 0# THEN FOR I = 0 TO M FOR J = 0 TO N CQM(I, J, 0) = 1D+200: CQM(I, J, 1) = 0# CQD(I, J, 0) = 1D+200: CQD(I, J, 1) = 0# NEXT J NEXT I EXIT SUB END IF XC = CABS(Z()) IF Z(1) = 0# OR XC < 1# THEN LS = 1 IF XC > 1# THEN LS = -1 ' ZQ=CSQRT(LS*(1.0-Z*Z)) CALL CMUL(Z(), Z(), TMP()): TMP(0) = LS * (1# - TMP(0)): TMP(1) = -LS * TMP(1) CALL 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# + Z(0)): TMP(1) = LS * Z(1) TMP1(0) = 1# - Z(0): TMP1(1) = -Z(1) CALL CDIV(TMP(), TMP1(), TMP2()) CALL CLOG(TMP2(), CQ0(), Err1): CQ0(0) = .5 * CQ0(0): CQ0(1) = .5 * CQ0(1) IF XC < 1.0001 THEN CQM(0, 0, 0) = CQ0(0): CQM(0, 0, 1) = CQ0(1) ' CQM(0,1)=Z*CQ0-1.0 CALL CMUL(Z(), CQ0(), TMP()) CQM(0, 1, 0) = TMP(0) - 1#: CQM(0, 1, 1) = TMP(1) ' CQM(1,0)=-1.0/ZQ TMP(0) = -1#: TMP(1) = 0#: CALL CDIV(TMP(), ZQ(), TMP1()) CQM(1, 0, 0) = TMP1(0): CQM(1, 0, 1) = TMP1(1) ' CQM(1,1)=-ZQ*(CQ0+Z/(1.0-Z*Z)) CALL CMUL(Z(), Z(), TMP()): TMP(0) = 1# - TMP(0): TMP(1) = -TMP(1) CALL CDIV(Z(), TMP(), TMP2()) TMP(0) = CQ0(0) + TMP2(0): TMP(1) = CQ0(1) + TMP2(1) CALL CMUL(ZQ(), TMP(), TMP1()) CQM(1, 1, 0) = -TMP1(0): CQM(1, 1, 1) = -TMP1(1) FOR I = 0 TO 1 FOR J = 2 TO N ' CQM(I,J)=((2.0*J-1.0)*Z*CQM(I,J-1)-(J+I-1.0)*CQM(I,J-2))/(J-I) TMP1(0) = CQM(I, J - 1, 0): TMP1(1) = CQM(I, J - 1, 1) CALL CMUL(Z(), TMP1(), TMP()) TMP(0) = (2# * J - 1#) * TMP(0): TMP(1) = (2# * J - 1#) * TMP(1) TMP1(0) = (J + I - 1!) * CQM(I, J - 2, 0): TMP1(1) = (J + I - 1!) * 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) NEXT J NEXT I FOR J = 0 TO N FOR I = 2 TO M ' 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) CALL CDIV(Z(), ZQ(), TMP()): TMP(0) = -2# * (I - 1!) * TMP(0): TMP(1) = -2# * (I - 1!) * TMP(1) TMP1(0) = CQM(I - 1, J, 0): TMP1(1) = CQM(I - 1, J, 1) CALL CMUL(TMP(), TMP1(), TMP2()) 'Now TMP2=-2.0*(I-1.0)*Z/ZQ*CQM(I-1,J) TMP1(0) = LS * (J + I - 1#) * (J - I + 2!) * CQM(I - 2, J, 0) TMP1(1) = LS * (J + I - 1#) * (J - I + 2!) * 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) NEXT I NEXT J ELSE IF XC > 1.1 THEN KM = 40 + M + N ELSE KM = (40 + M + N) * INT(-1# - 1.8 * LOG(XC - 1#)) END IF CQF2(0) = 0#: CQF2(1) = 0# CQF1(0) = 1#: CQF1(1) = 0# FOR K = KM TO 0 STEP -1 ' CQF0=((2*K+3.0)*Z*CQF1-(K+2.0)*CQF2)/(K+1.0) CALL CMUL(Z(), CQF1(), TMP()): TMP(0) = (2 * K + 3#) * TMP(0): TMP(1) = (2 * K + 3#) * TMP(1) TMP1(0) = (K + 2#) * CQF2(0): TMP1(1) = (K + 2#) * CQF2(1) CQF0(0) = (TMP(0) - TMP1(0)) / (K + 1#): CQF0(1) = (TMP(1) - TMP1(1)) / (K + 1#) IF K <= N THEN CQM(0, K, 0) = CQF0(0): CQM(0, K, 1) = CQF0(1) END IF CQF2(0) = CQF1(0): CQF2(1) = CQF1(1) CQF1(0) = CQF0(0): CQF1(1) = CQF0(1) NEXT K FOR K = 0 TO N ' CQM(0,K)=CQ0*CQM(0,K)/CQF0 TMP1(0) = CQM(0, K, 0): TMP1(1) = CQM(0, K, 1) CALL CMUL(CQ0(), TMP1(), TMP()) CALL CDIV(TMP(), CQF0(), TMP1()) CQM(0, K, 0) = TMP1(0): CQM(0, K, 1) = TMP1(1) NEXT K CQF2(0) = 0#: CQF2(1) = 0# CQF1(0) = 1#: CQF1(1) = 0# FOR K = KM TO 0 STEP -1 ' CQF0=((2*K+3.0)*Z*CQF1-(K+1.0)*CQF2)/(K+2.0) CALL CMUL(Z(), CQF1(), TMP()): TMP(0) = (2 * K + 3#) * TMP(0): TMP(1) = (2 * K + 3#) * TMP(1) TMP1(0) = (K + 1#) * CQF2(0): TMP1(1) = (K + 1#) * CQF2(1) CQF0(0) = (TMP(0) - TMP1(0)) / (K + 2#): CQF0(1) = (TMP(1) - TMP1(1)) / (K + 2#) IF K <= N THEN CQM(1, K, 0) = CQF0(0): CQM(1, K, 1) = CQF0(1) END IF CQF2(0) = CQF1(0): CQF2(1) = CQF1(1) CQF1(0) = CQF0(0): CQF1(1) = CQF0(1) NEXT K ' CQ10=-1.0/ZQ TMP(0) = -1#: TMP(1) = 0#: CALL CDIV(TMP(), ZQ(), CQ10()) FOR K = 0 TO N ' CQM(1,K)=CQ10*CQM(1,K)/CQF0 TMP1(0) = CQM(1, K, 0): TMP1(1) = CQM(1, K, 1) CALL CMUL(CQ10(), TMP1(), TMP()) CALL CDIV(TMP(), CQF0(), TMP1()) CQM(1, K, 0) = TMP1(0): CQM(1, K, 1) = TMP1(1) NEXT K FOR J = 0 TO N 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 TO M - 2 'CQF=-2.0*(I+1)*Z/ZQ*CQ1+(J-I)*(J+I+1.0)*CQ0 CALL CDIV(Z(), ZQ(), TMP()): TMP(0) = -2# * (I + 1) * TMP(0): TMP(1) = -2# * (I + 1) * TMP(1) CALL CMUL(TMP(), CQ1(), TMP2()) 'Now TMP2=-2.0*(I+1)*Z/ZQ*CQ1 TMP1(0) = (J - I) * (J + I + 1#) * CQ0(0) TMP1(1) = (J - I) * (J + I + 1#) * CQ0(1)'Now TMP1=(J-I)*(J+I+1.0)*CQ0 CQF(0) = TMP2(0) + TMP1(0): CQF(1) = TMP2(1) + TMP1(1) CQM(I + 2, J, 0) = CQF(0): CQM(I + 2, J, 1) = CQF(1) CQ0(0) = CQ1(0): CQ0(1) = CQ1(1) CQ1(0) = CQF(0): CQ1(1) = CQF(1) NEXT I NEXT J END IF ' CQD(0,0)=LS/ZS TMP(0) = 1# * LS: TMP(1) = 0#: CALL CDIV(TMP(), ZS(), TMP1()) CQD(0, 0, 0) = TMP1(0): CQD(0, 0, 1) = TMP1(1) FOR J = 1 TO N ' CQD(0,J)=LS*J*(CQM(0,J-1)-Z*CQM(0,J))/ZS TMP1(0) = CQM(0, J, 0): TMP1(1) = CQM(0, J, 1) CALL CMUL(Z(), TMP1(), TMP()) TMP(0) = LS * J * (CQM(0, J - 1, 0) - TMP(0)) TMP(1) = LS * J * (CQM(0, J - 1, 1) - TMP(1)) CALL CDIV(TMP(), ZS(), TMP1()) CQD(0, J, 0) = TMP1(0): CQD(0, J, 1) = TMP1(1) NEXT J FOR J = 0 TO N FOR I = 1 TO M ' CQD(I,J)=LS*I*Z/ZS*CQM(I,J)+(I+J)*(J-I+1.0)/ZQ*CQM(I-1,J) CALL CDIV(Z(), ZS(), TMP()): TMP(0) = LS * I * TMP(0): TMP(1) = LS * I * TMP(1) TMP1(0) = CQM(I, J, 0): TMP1(1) = CQM(I, J, 1) CALL CMUL(TMP(), TMP1(), TMP2()) 'Now TMP2=LS*I*Z/ZS*CQM(I,J) TMP1(0) = (I + J) * (J - I + 1#): TMP1(1) = 0# CALL CDIV(TMP1(), ZQ(), TMP()) TMP3(0) = CQM(I - 1, J, 0): TMP3(1) = CQM(I - 1, J, 1) CALL CMUL(TMP(), TMP3(), TMP1()) 'Now TMP1=(I+J)*(J-I+1.0)/ZQ*CQM(I-1,J) CQD(I, J, 0) = TMP2(0) + TMP1(0) CQD(I, J, 1) = TMP2(1) + TMP1(1) NEXT I NEXT J END SUB 'CLQMN() ' Z=Z1*Z2 SUB CMUL (Z1(), Z2(), Z()) Z(0) = Z1(0) * Z2(0) - Z1(1) * Z2(1) Z(1) = Z1(0) * Z2(1) + Z1(1) * Z2(0) END SUB SUB CSQRT (ZA(), ZB()) '***BEGIN PROLOGUE CSQRT ' ' DOUBLE PRECISION COMPLEX SQUARE ROOT, CSQRT(ZA,ZB) ' '***FUNCTIONS CALLED: CABS() '***END PROLOGUE CSQRT ' Labels: 12, 22, 31, 41, 51, 61, 71, 80 ' double AR,AI,BR,BI,ZM, DTHETA, DPI, DRT DRT = .7071067811865476# DPI = 4# * ATN(1#) AR = ZA(0): AI = ZA(1) ZM = CABS(ZA()) ZM = SQR(ZM) IF AR = 0# THEN GOTO 12 IF AI = 0# THEN GOTO 22 DTHETA = ATN(AI / AR) IF DTHETA <= 0# THEN GOTO 41 IF AR < 0# THEN DTHETA = DTHETA - DPI GOTO 51 12 IF AI > 0# THEN GOTO 61 IF AI < 0# THEN GOTO 71 BR = 0# BI = 0# GOTO 80 22 IF AR > 0# THEN GOTO 31 BR = 0# BI = SQR(ABS(AR)) GOTO 80 31 BR = SQR(AR) BI = 0# GOTO 80 41 IF AR < 0# THEN DTHETA = DTHETA + DPI 51 DTHETA = DTHETA * .5 BR = ZM * COS(DTHETA) BI = ZM * SIN(DTHETA) GOTO 80 61 BR = ZM * DRT BI = ZM * DRT GOTO 80 71 BR = ZM * DRT BI = -ZM * DRT 80 ZB(0) = BR: ZB(1) = BI END SUB 'CSQRT