DECLARE SUB CLQN (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#()) '*********************************************************************** '* Calculate the Legendre polynomials for a complex argument * '* ------------------------------------------------------------------- * '* EXPLANATION: * '* * '* ========================================================== * '* Purpose: This program computes the Legendre polynomials * '* of second kind 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, y = 2 * '* * '* 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". * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '*********************************************************************** 'PROGRAM MCLQN DEFDBL A-H, O-Z DEFINT I-N OPTION BASE 0 F$ = "###.######## " CLS PRINT INPUT " Please enter Nmax, x and y (z=x+iy): ", N, X, Y 'Tables of N+1 complex numbers DIM SHARED CQN(N + 1, 2), CQD(N + 1, 2) PRINT " x ="; X; ", y ="; Y PRINT CALL CLQN(N, X, Y) PRINT " n Re[Qn(z)] Im[Qn(z)] Re[Qn'(z)] Im[Qn'(z)] " PRINT " -----------------------------------------------------------" FOR K = 0 TO N PRINT " "; K; PRINT USING F$; CQN(K, 0); PRINT USING F$; CQN(K, 1); PRINT USING F$; CQD(K, 0); PRINT USING F$; CQD(K, 1) NEXT K PRINT END 'of main program 'ABSOLUTE VALUE OF A COMPLEX NUMBER Z=X+I*Y FUNCTION CABS (Z()) U = ABS(Z(0)) V = ABS(Z(1)) S = U + V '-------------------------------------------------------------------- ' S*1.0 MAKES AN UNNORMALIZED UNDERFLOW ON CDC MACHINES INTO A ' TRUE FLOATING ZERO '-------------------------------------------------------------------- S = S * 1# 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()) 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 ZABS '***END PROLOGUE CLOG DPI = 3.141592653589793# DHPI = 1.570796326794897# 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 CLQN (N, X, Y) ' ======================================================= ' 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,...,N ' Output: CQN(n) --- Qn(z) ' CQD(n) --- Qn'(z) ' Note: CQN, CQD are shared with main program. ' ======================================================= ' Complex numbers DIM CONE(2), CQ0(2), CQ1(2), CQF0(2), CQF1(2), CQF2(2), Z(2) DIM TCQ(2), TMP(2), TMP1(2), TMP2(2) DIM ERR1 AS INTEGER CONE(0) = 1#: CONE(1) = 0# Z(0) = X: Z(1) = Y IF Z(0) = 1# AND Z(1) = 0# THEN FOR K = 0 TO N CQN(K, 0) = 1D+200: CQN(K, 1) = 0# CQD(K, 0) = 1D+200: CQD(K, 1) = 0# NEXT K EXIT SUB END IF LS = 1 IF CABS(Z()) > 1# THEN LS = -1 'CQ0=0.5*CLOG(LS*(1#+Z)/(1#-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) CALL CDIV(TMP(), TMP1(), TMP2()) CALL CLOG(TMP2(), CQ0(), ERR1) CQ0(0) = CQ0(0) * .5: CQ0(1) = CQ0(1) * .5 'CQ1=Z*CQ0-1# CALL 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# THEN CQF0(0) = CQ0(0): CQF0(1) = CQ0(1) CQF1(0) = CQ1(0): CQF1(1) = CQ1(1) FOR K = 2 TO N 'CQF2=((2.0*K-1#)*Z*CQF1-(K-1#)*CQF0)/K CALL CMUL(Z(), CQF1(), TMP()) TMP(0) = TMP(0) * (2# * K - 1#) TMP(1) = TMP(1) * (2# * K - 1#) TMP1(0) = (K - 1#) * CQF0(0) TMP1(1) = (K - 1#) * 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) NEXT K ELSE IF CABS(Z()) > 1.1 THEN KM = 40 + N ELSE TMP(0) = Z(0) - 1#: TMP(1) = Z(1) KM = (40 + N) * INT(-1# - 1.8 * LOG(CABS(TMP()))) 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#)*CQF2)/(K+1#) CALL CMUL(Z(), CQF1(), TMP()) TMP(0) = TMP(0) * (2# * K + 3#) TMP(1) = TMP(1) * (2# * K + 3#) TMP1(0) = (K + 2#) * CQF2(0) TMP1(1) = (K + 2#) * CQF2(1) CQF0(0) = (TMP(0) - TMP1(0)) / (1# + K) CQF0(1) = (TMP(1) - TMP1(1)) / (1# + K) IF K <= N THEN CQN(K, 0) = CQF0(0): CQN(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 'CQN(K)=CQN(K)*CQ0/CQF0 TCQ(0) = CQN(K, 0): TCQ(1) = CQN(K, 1) CALL CMUL(TCQ(), CQ0(), TMP()) CALL CDIV(TMP(), CQF0(), TCQ()) CQN(K, 0) = TCQ(0): CQN(K, 1) = TCQ(1) NEXT K END IF 'CQD(0)=(CQN(0)-Z*CQN(0))/(Z*Z-1#) TCQ(0) = CQN(0, 0): TCQ(1) = CQN(0, 1) CALL CMUL(Z(), TCQ(), TMP()): TMP(0) = CQN(1, 0) - TMP(0): TMP(1) = CQN(1, 1) - TMP(1) CALL CMUL(Z(), Z(), TMP1()): TMP1(0) = TMP1(0) - 1# CALL CDIV(TMP(), TMP1(), TCQ()) CQD(0, 0) = TCQ(0): CQD(0, 1) = TCQ(1) FOR K = 1 TO N 'CQD(K)=(K*Z*CQN(K)-K*CQN(K-1))/(Z*Z-1#) TCQ(0) = CQN(K, 0): TCQ(1) = CQN(K, 1) CALL CMUL(Z(), TCQ(), 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) CALL CDIV(TMP(), TMP1(), TCQ()) CQD(K, 0) = TCQ(0): CQD(K, 1) = TCQ(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