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