DECLARE SUB CLPN (N%, X#, Y#) DECLARE SUB CDIV (Z1#(), Z2#(), Z#()) DECLARE SUB CMUL (Z1#(), Z2#(), Z#()) '*********************************************************************** '* 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 , y = 2 * '* * '* 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". * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '*********************************************************************** 'PROGRAM MCLPN DEFDBL A-H, O-Z DEFINT I-N 'Tables of Complex numbers 'double CPN(N+1,2), CPD(N+1,2) 'integer K,N 'double X,Y CLS PRINT INPUT " Please enter Nmax, x and y (z=x+iy): ", N, X, Y PRINT " x ="; X; ", y ="; Y DIM SHARED CPN(N + 1, 2), CPD(N + 1, 2) CALL CLPN(N, X, Y) F\$ = " #####.######" PRINT PRINT " n Re(Pn(z)) Im(Pn(z)) Re(Pn'(z)) Im(Pn'(z)) " PRINT " -------------------------------------------------------------" FOR K = 0 TO N PRINT " "; K; PRINT USING F\$; CPN(K, 0); PRINT USING F\$; CPN(K, 1); PRINT " "; PRINT USING F\$; CPD(K, 0); PRINT USING F\$; CPD(K, 1) NEXT K PRINT END 'of main program ' 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 CLPN (N, X, Y) ' ========================================================= ' 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) (1) ' Complex CPD(n) --- Pn'(z) (1) ' (1) shared with main program. ' ========================================================= DIM CP0(2), CP1(2), CPF(2), TMP(2), TMP1(2), TMP2(2), Z(2)'Complex numbers Z(0) = X: Z(1) = Y CPN(0, 0) = 1#: CPN(0, 1) = 0# CPN(1, 0) = Z(1): CPN(1, 1) = Z(2) CPD(0, 0) = 0#: CPD(0, 1) = 0# CPD(1, 0) = 1#: CPD(1, 1) = 0# CP0(0) = 1#: CP0(1) = 0# CP1(0) = Z(0): CP1(1) = Z(1) FOR K = 2 TO N 'CPF=(2.0*K-1.0)/K*Z*CP1-(K-1.0)/K*CP0 CALL CMUL(Z(), CP1(), TMP()) TMP(0) = (2# * K - 1#) / K * TMP(0): TMP(1) = (2# * K - 1#) / K * TMP(1) CPF(0) = TMP(0) - (K - 1#) / K * CP0(0) CPF(1) = TMP(1) - (K - 1#) / K * CP0(1) CPN(K, 0) = CPF(0): CPN(K, 1) = CPF(1) IF ABS(X) = 1# AND Y = 0# THEN 'CPD(K)=0.5*X^(K+1)*K*(K+1.0) CPD(K, 0) = .5 * X ^ K + 1 * K * (K + 1#) CPD(K, 1) = 0# ELSE 'CPD(K)=K*(CP1-Z*CPF)/(1.0-Z*Z) CALL CMUL(Z(), CPF(), TMP()) TMP(0) = K * (CP1(0) - TMP(0)): TMP(1) = K * (CP1(1) - TMP(1)) CALL CMUL(Z(), Z(), TMP1()): TMP1(0) = 1# - TMP1(0): TMP1(1) = -TMP1(1) CALL CDIV(TMP(), TMP1(), TMP2()) CPD(K, 0) = TMP2(0): CPD(K, 1) = TMP2(1) END IF CP0(0) = CP1(0): CP0(1) = CP1(1) CP1(0) = CPF(0): CP1(1) = CPF(1) NEXT K 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