!***********************************************************************
!*      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:  CPN(n) --- Pn(z)                                      *
!*               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.0,  y =  2.0                                              *
!*                                                                     *
!*    n     Re[Pn(z)]     Im[Pn(z)]      Re[Pn'(z)]   Im[Pn'(z)]       *
!*   -------------------------------------------------------------     *
!*    0   0.100000D+01  0.000000D+00   0.000000D+00  0.000000D+00      *
!*    1   0.300000D+01  0.200000D+01   0.100000D+01  0.000000D+00      *
!*    2   0.700000D+01  0.180000D+02   0.900000D+01  0.600000D+01      *
!*    3  -0.270000D+02  0.112000D+03   0.360000D+02  0.900000D+02      *
!*    4  -0.539000D+03  0.480000D+03  -0.180000D+03  0.790000D+03      *
!*    5  -0.461700D+04  0.562000D+03  -0.481500D+04  0.441000D+04      *
!*                                                                     *
!* ------------------------------------------------------------------- *
!* REFERENCE: "Fortran Routines for Computation of Special Functions,  *
!*             jin.ece.uiuc.edu/routines/routines.html".               *
!*                                                                     *
!*                                   F90 Release By J-P Moreau, Paris. *
!*                                          (www.jpmoreau.fr)          *
!***********************************************************************
PROGRAM MCLPN

IMPLICIT DOUBLE PRECISION (X,Y)
IMPLICIT COMPLEX *16 (C,Z)
DIMENSION CPN(0:100),CPD(0:100)

   PRINT *,' '
   WRITE(*,10,advance='no')
   READ *, N, X, Y
   WRITE(*,30)X,Y
   PRINT *,' '

   CALL CLPN(N,X,Y,CPN,CPD)

   WRITE(*,*)'  n    Re[Pn(z)]     Im[Pn(z)]     Re[Pn''(z)]',  &
   '   Im[Pn''(z)]'
   WRITE(*,*)' ---------------------------------------------',  &
   '--------------'
   DO K=0,N
     WRITE(*,20)K,CPN(K),CPD(K)
   END DO
   PRINT *,' '
   STOP
10 FORMAT('  Please enter Nmax, x and y (z=x+iy): ')
20 FORMAT(1X,I3,4D14.6)
30 FORMAT(3X,'x =',F5.1,',  ','y =',F5.1)
   END

   SUBROUTINE CLPN(N,X,Y,CPN,CPD)
!  ==================================================
!  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:  CPN(n) --- Pn(z)
!           CPD(n) --- Pn'(z)
!  ==================================================
   IMPLICIT DOUBLE PRECISION (X,Y)
   IMPLICIT COMPLEX *16 (C,Z)
   DIMENSION CPN(0:N),CPD(0:N)

   Z=CMPLX(X,Y)
   CPN(0)=(1.0D0,0.0D0)
   CPN(1)=Z
   CPD(0)=(0.0D0,0.0D0)
   CPD(1)=(1.0D0,0.0D0)
   CP0=(1.0D0,0.0D0)
   CP1=Z

   DO 10 K=2,N
     CPF=(2.0D0*K-1.0D0)/K*Z*CP1-(K-1.0D0)/K*CP0
     CPN(K)=CPF
     IF (DABS(X).EQ.1.0D0.AND.Y.EQ.0.0D0) THEN
       CPD(K)=0.5D0*X**(K+1)*K*(K+1.0D0)
     ELSE
       CPD(K)=K*(CP1-Z*CPF)/(1.0D0-Z*Z)
     ENDIF
     CP0=CP1
10   CP1=CPF
   RETURN
   END

!end of file mclpn.f90