/**********************************************************************
*      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.0,  y =  2.0                                              *
*                                                                     *
*    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".               *
*                                                                     *
*                                   C++ Release By J-P Moreau, Paris. *
**********************************************************************/
#include <stdio.h>
#include <math.h>

#define NMAX  100

typedef double COMPLEX[2];

COMPLEX CPN[NMAX], CPD[NMAX];
int     K,N;
double  X,Y;


      // Z=Z1/Z2
      void CDIV(COMPLEX Z1, COMPLEX Z2, COMPLEX Z) {
      double D;
        D=Z2[0]*Z2[0]+Z2[1]*Z2[1];
        if (D<1e-12) return;
        Z[0]=(Z1[0]*Z2[0]+Z1[1]*Z2[1])/D;
        Z[1]=(Z1[1]*Z2[0]-Z1[0]*Z2[1])/D;
      }

      // Z=Z1*Z2
      void CMUL(COMPLEX Z1, COMPLEX Z2, COMPLEX Z) {
        Z[0]=Z1[0]*Z2[0] - Z1[1]*Z2[1];
        Z[1]=Z1[0]*Z2[1] + Z1[1]*Z2[0];
      }

      void CLPN(int N, double X, double Y, COMPLEX *CPN, COMPLEX *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:  Complex CPN(n) --- Pn(z)
                 Complex CPD(n) --- Pn'(z)
      ========================================================= */
        COMPLEX  CP0,CP1,CPF,TMP,TMP1,Z;

        Z[0]=X; Z[1]=Y;
        CPN[0][0]=1.0; CPN[0][1]=0.0;
        CPN[1][0]=Z[1]; CPN[1][1]=Z[2];
        CPD[0][0]=0.0; CPD[0][1]=0.0;
        CPD[1][0]=1.0; CPD[1][1]=0.0;
        CP0[0]=1.0; CP0[1]=0.0;
        CP1[0]=Z[0]; CP1[1]=Z[1];
        for (K=2; K<=N; K++) {
           //CPF=(2.0*K-1.0)/K*Z*CP1-(K-1.0)/K*CP0
           CMUL(Z,CP1,TMP);
           TMP[0]=(2.0*K-1.0)/K*TMP[0]; TMP[1]=(2.0*K-1.0)/K*TMP[1];
           CPF[0]=TMP[0]-(K-1.0)/K*CP0[0];
           CPF[1]=TMP[1]-(K-1.0)/K*CP0[1];
           CPN[K][0]=CPF[0]; CPN[K][1]=CPF[1];
           if (fabs(X) == 1.0 && Y == 0.0) {
             //CPD(K)=0.5*X^(K+1)*K*(K+1.0)
             CPD[K][0]=0.5*pow(X,K+1)*K*(K+1.0);
             CPD[K][1]=0.0;
           }
           else {
             //D(K)=K*(CP1-Z*CPF)/(1.0-Z*Z)
             CMUL(Z,CPF,TMP);
             TMP[0]=K*(CP1[0]-TMP[0]); TMP[1]=K*(CP1[1]-TMP[1]);
             CMUL(Z,Z,TMP1); TMP1[0]=1.0-TMP1[0]; TMP1[1]=-TMP1[1];
             CDIV(TMP,TMP1,CPD[K]);
           }
           CP0[0]=CP1[0]; CP0[1]=CP1[1];
           CP1[0]=CPF[0]; CP1[1]=CPF[1];
        }
      }


      void main()  {

        printf("\n  Please enter Nmax, x and y (z=x+iy): "); 
		scanf("%d %lf %lf", &N, &X, &Y);
        printf("   x =%5.1f,  y =%5.1f\n\n", X, Y);

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

        printf("   n    Re[Pn(z)]     Im[Pn(z)]      Re[Pn'(z)]    Im[Pn'(z)]\n");
        printf("  ------------------------------------------------------------\n");

        for (K=0; K<=N; K++)
          printf("%4d%14.6f%14.6f %14.6f%14.6f\n", K, CPN[K][0],CPN[K][1],CPD[K][0],CPD[K][1]);
 
        printf("\n");

      }

// end of file mclpn.cpp