DECLARE SUB CLPMN (M%, N%, X#, Y#) DECLARE FUNCTION CABS# (Z#()) DECLARE SUB CDIV (Z1#(), Z2#(), Z#()) DECLARE SUB CMUL (Z1#(), Z2#(), Z#()) DECLARE SUB CSQRT (ZA#(), ZB#()) '************************************************************************** '* 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 = .5, y = .2 * '* * '* 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". * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '************************************************************************** DEFDBL A-H, O-Z DEFINT I-N OPTION BASE 0 ' int 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 DIM SHARED CPM(M + 1, N + 1, 2), CPD(M + 1, N + 1, 2) CALL CLPMN(M, N, X, Y) F\$ = " ######.######" PRINT " m n Re[Pmn(z)] Im[Pmn(z)] Re[Pmn'(z)] Im[Pmn'(z)] " PRINT " ---------------------------------------------------------------------" FOR J = 0 TO M PRINT " "; J; " "; N; PRINT USING F\$; CPM(J, N, 0); PRINT USING F\$; CPM(J, N, 1); PRINT USING F\$; CPD(J, N, 0); PRINT USING F\$; CPD(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 CLPMN (M, N, X, Y) ' ================================================================ ' 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) ' ================================================================ DIM TMP(2), TMP1(2), TMP2(2), Z(2), ZQ(2), ZS(2) 'Complex numbers ' int I,J,LS Z(0) = X: Z(1) = Y FOR I = 0 TO N FOR J = 0 TO M CPM(J, I, 0) = 0#: CPM(J, I, 1) = 0# CPD(J, I, 0) = 0#: CPD(J, I, 1) = 0# NEXT J NEXT I CPM(0, 0, 0) = 1# IF ABS(X) = 1# AND Y = 0# THEN FOR I = 1 TO N CPM(0, I, 0) = X ^ I: CPM(0, I, 1) = 0# CPD(0, I, 0) = .5 * I * (I + 1) * X ^ (I + 1) CPD(0, I, 1) = 0# NEXT I FOR J = 1 TO N FOR I = 1 TO M IF I = 1 THEN CPD(I, J, 0) = 1D+200 CPD(I, J, 1) = 0# ELSEIF I = 2 THEN CPD(I, J, 0) = -.25 * (J + 2) * (J + 1) * J * (J - 1) * X ^ (J + 1) CPD(I, J, 1) = 0# END IF NEXT I NEXT J RETURN END IF LS = 1 IF CABS(Z()) > 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.0D0-Z*Z) ZS(0) = TMP(0): ZS(1) = TMP(1) FOR I = 1 TO M ' CPM(I,I)=-LS*(2.0*I-1.0)*ZQ*CPM(I-1,I-1) TMP1(0) = CPM(I - 1, I - 1, 0): TMP1(1) = CPM(I - 1, I - 1, 1) CALL CMUL(ZQ(), TMP1(), TMP()) CPM(I, I, 0) = -LS * (2# * I - 1#) * TMP(0) CPM(I, I, 1) = -LS * (2# * I - 1#) * TMP(1) NEXT I FOR I = 0 TO M ' CPM(I,I+1)=(2.0*I+1.0)*Z*CPM(I,I) TMP1(0) = CPM(I, I, 0): TMP1(1) = CPM(I, I, 1) CALL CMUL(Z(), TMP1(), TMP()) CPM(I, I + 1, 0) = (2# * I + 1#) * TMP(0) CPM(I, I + 1, 1) = (2# * I + 1#) * TMP(1) NEXT I FOR I = 0 TO M FOR J = I + 2 TO N 'CPM(I,J)=((2.0*J-1.0)*Z*CPM(I,J-1)-(I+J-1.0)*CPM(I,J-2))/(J-I) TMP1(0) = CPM(I, J - 1, 0): TMP1(1) = CPM(I, J - 1, 1) CALL CMUL(Z(), TMP1(), TMP()) CPM(I, J, 0) = ((2# * J - 1#) * TMP(0) - (I + J - 1#) * CPM(I, J - 2, 0)) / (J - I) CPM(I, J, 1) = ((2# * J - 1#) * TMP(1) - (I + J - 1#) * CPM(I, J - 2, 1)) / (J - I) NEXT J NEXT I CPD(0, 0, 0) = 0#: CPD(0, 0, 1) = 0# FOR J = 1 TO N 'CPD(0,J)=LS*J*(CPM(0,J-1)-Z*CPM(0,J))/ZS TMP1(0) = CPM(0, J, 0): TMP1(1) = CPM(0, J, 1) CALL CMUL(Z(), TMP1(), TMP()) TMP(0) = LS * J * (CPM(0, J - 1, 0) - TMP(0)) TMP(1) = LS * J * (CPM(0, J - 1, 1) - TMP(1)) CALL CDIV(TMP(), ZS(), TMP1()) CPD(0, J, 0) = TMP1(0): CPD(0, J, 1) = TMP1(1) NEXT J FOR I = 1 TO M FOR J = I TO N 'CPD(I,J)=LS*I*Z*CPM(I,J)/ZS+(J+I)*(J-I+1.0)*CPM(I-1,J)/ZQ TMP1(0) = CPM(I, J, 0): TMP1(1) = CPM(I, J, 1) CALL CMUL(Z(), TMP1(), TMP()) TMP(0) = TMP(0) * LS * I: TMP(1) = TMP(1) * LS * I CALL CDIV(TMP(), ZS(), TMP2()) TMP(0) = TMP2(0): TMP(1) = TMP2(1) 'Now TMP=LS*I*Z*CPM(I,J)/ZS TMP1(0) = (J + I) * (J - I + 1#) * CPM(I - 1, J, 0) TMP1(1) = (J + I) * (J - I + 1#) * CPM(I - 1, J, 1) CALL 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) NEXT J NEXT I END SUB ' 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: 10, 20, 30, 40, 50, 60, 70, 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 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 IF AI < 0# THEN GOTO 70 BR = 0# BI = 0# GOTO 80 21 IF AR > 0# THEN GOTO 30 BR = 0# BI = SQR(ABS(AR)) GOTO 80 30 BR = SQR(AR) BI = 0# GOTO 80 40 IF AR < 0# THEN DTHETA = DTHETA + DPI 50 DTHETA = DTHETA * .5 BR = ZM * COS(DTHETA) BI = ZM * SIN(DTHETA) GOTO 80 60 BR = ZM * DRT BI = ZM * DRT GOTO 80 70 BR = ZM * DRT BI = -ZM * DRT 80 : ZB(0) = BR: ZB(1) = BI END SUB 'CSQRT