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