/*************************************************************************
*   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 = 0.500000, y = 0.200000                              *
*                                                                        *
*    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".                  *
*                                                                        *
*                                      C++ Release By J-P Moreau, Paris. *
*                                            (www.jpmoreau.fr)           *
*************************************************************************/ 
#include <stdio.h>
#include <math.h>

#define NMAX  25

  typedef double COMPLEX[2];
  typedef COMPLEX CTab[NMAX][NMAX];

  CTab    CQM, CQD;
  int     J,M,N;
  double  X,Y;


      //ABSOLUTE VALUE OF A COMPLEX NUMBER Z=X+I*Y
      double CABS(COMPLEX Z) {
//    Labels: e10, e20
      double U,V,Q,S;
        U = fabs(Z[0]);
        V = fabs(Z[1]);
        S = U+V;
/*-------------------------------------------------------------------
      S*1.0 MAKES AN UNNORMALIZED UNDERFLOW ON CDC MACHINES INTO A
      TRUE FLOATING ZERO
-------------------------------------------------------------------*/
        S = S*1.0;
        if (S == 0.0) goto e20;
        if (U > V) goto e10;
        Q = U/V;
        return (V*sqrt(1.0+Q*Q));
e10:    Q = V/U;
        return (U*sqrt(1.0+Q*Q));
e20:    return 0.0;
      }

      // 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;
      }

      void CLOG(COMPLEX ZA, COMPLEX ZB, int *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    */
      // Labels: e10,e20,e30,e40,e50,e60,e70
        double AR,AI,BR,BI,ZM, DTHETA, DPI, DHPI;

        DPI = 3.141592653589793238462643383E+0000;
        DHPI= 1.570796326794896619231321696E+0000;

        *IERR=0;
        AR=ZA[0]; AI=ZA[1];
        if (AR == 0.0) goto e10;
        if (AI == 0.0) goto e20;
        DTHETA = atan(AI/AR);
        if (DTHETA <= 0.0) goto e40;
        if (AR < 0.0) DTHETA -= DPI;
        goto e50;
e10:    if (AI == 0.0) goto e60;
        BI = DHPI;
        BR = log(fabs(AI));
        if (AI < 0.0) BI = -BI;
        goto e70;
e20:    if (AR > 0.0) goto e30;
        BR = log(fabs(AR));
        BI = DPI;
        goto e70;
e30:    BR = log(AR);
        BI = 0.0;
        goto e70;
e40:    if (AR < 0.0) DTHETA += DPI;
e50:    ZM = CABS(ZA);
        BR = log(ZM);
        BI = DTHETA;
        goto e70;
e60:    *IERR=1;
e70:    ZB[0]=BR; ZB[1]=BI;
      } //CLOG()

      // 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 CSQRT(COMPLEX ZA, COMPLEX ZB) {
      /***BEGIN PROLOGUE  CSQRT
      !
      !   DOUBLE PRECISION COMPLEX SQUARE ROOT, CSQRT(ZA,ZB)
      !
      !***FUNCTIONS CALLED:  CABS()
      !***END PROLOGUE  CSQRT    */
//    Labels: e10,e20,e30,e40,e50,e60,e70,e80
      double AR,AI,BR,BI,ZM, DTHETA, DPI, DRT;

	DRT=7.071067811865475244008443621e-001;
        DPI=3.141592653589793238462643383e+000;
        AR=ZA[0]; AI=ZA[1];
        ZM = CABS(ZA);
        ZM = sqrt(ZM);
        if (AR == 0.0) goto e10;
        if (AI == 0.0) goto e20;
        DTHETA = atan(AI/AR);
        if (DTHETA <= 0.0) goto e40;
        if (AR < 0.0)  DTHETA -= DPI;
        goto e50;
e10:    if (AI > 0.0) goto e60;
        if (AI < 0.0) goto e70;
        BR = 0.0;
        BI = 0.0;
        goto e80;
e20:    if (AR > 0.0) goto e30;
        BR = 0.0;
        BI = sqrt(fabs(AR));
        goto e80;
e30:    BR = sqrt(AR);
        BI = 0.0;
        goto e80;
e40:    if (AR < 0.0)  DTHETA += DPI;
e50:    DTHETA *= 0.5;
        BR = ZM*cos(DTHETA);
        BI = ZM*sin(DTHETA);
        goto e80;
e60:    BR = ZM*DRT;
        BI = ZM*DRT;
        goto e80;
e70:    BR = ZM*DRT;
        BI = -ZM*DRT;
e80   : ZB[0]=BR; ZB[1]=BI;
      } //CSQRT()


      void CLQMN(int M, int N, double X, double Y, CTab CQM, CTab CQD) {
/*    ========================================================
        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:   */
        double CQ0[2],CQ1[2],CQ10[2],CQF[2],CQF0[2],CQF1[2];
	double CQF2[2],TMP[2],TMP1[2],TMP2[2],Z[2],ZQ[2],ZS[2];

        int Err,I,K,KM,LS;
        double XC;

        Z[0]=X; Z[1]=Y;
        if (fabs(X) == 1.0 && Y == 0.0) {
           for (I=0; I<=M; I++)
	     for (J=0; J<=N; J++) {
               CQM[I][J][0]=1e+200; CQM[I][J][1]=0.0;
               CQD[I][J][0]=1e+200; CQD[I][J][1]=0.0;
             }
           return;
        }
        XC=CABS(Z);
        if (Z[1] == 0.0 || XC < 1.0)  LS=1;
        if (XC > 1.0)  LS=-1;
        // ZQ=CSQRT(LS*(1.0-Z*Z))
        CMUL(Z,Z,TMP); TMP[0]=LS*(1.0-TMP[0]); TMP[1]=-LS*TMP[1];
        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.0+Z[0]); TMP[1]=LS*Z[1];
        TMP1[0]=1.0-Z[0]; TMP1[1]=-Z[1];
        CDIV(TMP,TMP1,TMP2);
        CLOG(TMP2,CQ0,&Err); CQ0[0]=0.5*CQ0[0]; CQ0[1]=0.5*CQ0[1];
	if (XC < 1.0001) {
           CQM[0][0][0]=CQ0[0]; CQM[0][0][1]=CQ0[1];
           // CQM(0,1)=Z*CQ0-1.0
           CMUL(Z,CQ0,TMP);
           CQM[0][1][0]=TMP[0]-1.0; CQM[0][1][1]=TMP[1];
           // CQM(1,0)=-1.0/ZQ
           TMP[0]=-1.0; TMP[1]=0.0; CDIV(TMP,ZQ,CQM[1][0]);
           // CQM(1,1)=-ZQ*(CQ0+Z/(1.0-Z*Z))
           CMUL(Z,Z,TMP); TMP[0]=1.0-TMP[0]; TMP[1]=-TMP[1];
           CDIV(Z,TMP,TMP2);
           TMP[0]=CQ0[0]+TMP2[0]; TMP[1]=CQ0[1]+TMP2[1];
           CMUL(ZQ,TMP,CQM[1][1]);
           CQM[1][1][0]=-CQM[1][1][0]; CQM[1][1][1]=-CQM[1][1][1];
           for (I=0; I<2; I++)
	     for (J=2; J<=N; J++) {
               // CQM(I,J)=((2.0*J-1.0)*Z*CQM(I,J-1)-(J+I-1.0)*CQM(I,J-2))/(J-I)
               CMUL(Z,CQM[I][J-1],TMP); TMP[0]=(2.0*J-1.0)*TMP[0]; TMP[1]=(2.0*J-1.0)*TMP[1];
               TMP1[0]=(J+I-1.0)*CQM[I][J-2][0]; TMP1[1]=(J+I-1.0)*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);
             }
           for (J=0; J<=N; J++)
	     for (I=2; I<=M; I++) {
               // 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)
               CDIV(Z,ZQ,TMP); TMP[0]=-2.0*(I-1.0)*TMP[0]; TMP[1]=-2.0*(I-1.0)*TMP[1];
               CMUL(TMP,CQM[I-1][J],TMP2); //Now TMP2=-2.0*(I-1.0)*Z/ZQ*CQM(I-1,J)
               TMP1[0]=LS*(J+I-1.0)*(J-I+2.0)*CQM[I-2][J][0];
               TMP1[1]=LS*(J+I-1.0)*(J-I+2.0)*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];
             }       
        }
        else {
           if (XC > 1.1)
              KM=40+M+N;
           else
              KM=(int) ((40+M+N)*floor(-1.0-1.8*log(XC-1.0)));
           CQF2[0]=0.0; CQF2[1]=0.0;
           CQF1[0]=1.0; CQF1[1]=0.0;
           for (K=KM; K>=0; K--) {
             // CQF0=((2*K+3.0)*Z*CQF1-(K+2.0)*CQF2)/(K+1.0)
             CMUL(Z,CQF1,TMP); TMP[0]=(2*K+3.0)*TMP[0]; TMP[1]=(2*K+3.0)*TMP[1];
             TMP1[0]=(K+2.0)*CQF2[0]; TMP1[1]=(K+2.0)*CQF2[1];
             CQF0[0]=(TMP[0]-TMP1[0])/(K+1.0); CQF0[1]=(TMP[1]-TMP1[1])/(K+1.0);
             if (K <= N) {
               CQM[0][K][0]=CQF0[0]; CQM[0][K][1]=CQF0[1];
             }
             CQF2[0]=CQF1[0]; CQF2[1]=CQF1[1];
             CQF1[0]=CQF0[0]; CQF1[1]=CQF0[1];
           }
           for (K=0; K<=N; K++) {
             // CQM(0,K)=CQ0*CQM(0,K)/CQF0
             CMUL(CQ0,CQM[0][K],TMP);
             CDIV(TMP,CQF0,CQM[0][K]);
           }
           CQF2[0]=0.0; CQF2[1]=0.0;
           CQF1[0]=1.0; CQF1[1]=0.0;
           for (K=KM; K>=0; K--) {
             // CQF0=((2*K+3.0)*Z*CQF1-(K+1.0)*CQF2)/(K+2.0)
             CMUL(Z,CQF1,TMP); TMP[0]=(2*K+3.0)*TMP[0]; TMP[1]=(2*K+3.0)*TMP[1];
             TMP1[0]=(K+1.0)*CQF2[0]; TMP1[1]=(K+1.0)*CQF2[1];
             CQF0[0]=(TMP[0]-TMP1[0])/(K+2.0); CQF0[1]=(TMP[1]-TMP1[1])/(K+2.0);
             if (K <= N) {
               CQM[1][K][0]=CQF0[0]; CQM[1][K][1]=CQF0[1];
             }
             CQF2[0]=CQF1[0]; CQF2[1]=CQF1[1];
             CQF1[0]=CQF0[0]; CQF1[1]=CQF0[1];
           } 
           // CQ10=-1.0/ZQ
           TMP[0]=-1.0; TMP[1]=0.0; CDIV(TMP,ZQ,CQ10);
           for (K=0; K<=N; K++) {
             // CQM(1,K)=CQ10*CQM(1,K)/CQF0
             CMUL(CQ10,CQM[1][K],TMP);
             CDIV(TMP,CQF0,CQM[1][K]);
           }
           for (J=0; J<=N; J++) {
             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; I<M-1; I++) {
               //CQF=-2.0*(I+1)*Z/ZQ*CQ1+(J-I)*(J+I+1.0)*CQ0
               CDIV(Z,ZQ,TMP); TMP[0]=-2.0*(I+1)*TMP[0]; TMP[1]=-2.0*(I+1)*TMP[1];
               CMUL(TMP,CQ1,TMP2);             //Now TMP2=-2.0*(I+1)*Z/ZQ*CQ1
               TMP1[0]=(J-I)*(J+I+1.0)*CQ0[0];
               TMP1[1]=(J-I)*(J+I+1.0)*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];
             }
           }
        }
        // CQD(0,0)=LS/ZS
        TMP[0]=1.0*LS; TMP[1]=0.0; CDIV(TMP,ZS,CQD[0][0]);
        for (J=1; J<=N; J++) {
          // CQD(0,J)=LS*J*(CQM(0,J-1)-Z*CQM(0,J))/ZS
          CMUL(Z,CQM[0][J],TMP);
          TMP[0]=LS*J*(CQM[0][J-1][0]-TMP[0]);
          TMP[1]=LS*J*(CQM[0][J-1][1]-TMP[1]);
          CDIV(TMP,ZS,CQD[0][J]);
        }
        for (J=0; J<=N; J++)
	  for (I=1; I<=M; I++) {
            // CQD(I,J)=LS*I*Z/ZS*CQM(I,J)+(I+J)*(J-I+1.0)/ZQ*CQM(I-1,J)
            CDIV(Z,ZS,TMP); TMP[0]=LS*I*TMP[0]; TMP[1]=LS*I*TMP[1];
            CMUL(TMP,CQM[I][J],TMP2);    //Now TMP2=LS*I*Z/ZS*CQM(I,J)
            TMP1[0]=(I+J)*(J-I+1.0); TMP1[1]=0.0;
            CDIV(TMP1,ZQ,TMP);
            CMUL(TMP,CQM[I-1][J],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];
          };
        } //CLQMN()


      void main()  {

        printf("\n  Please enter m, n, x and y: "); 
		scanf("%d %d %lf %lf", &M, &N, &X, &Y);
        printf("  m = %2d, n = %2d, x = %f, y = %f\n\n", M, N, X, Y);

        CLQMN(M,N,X,Y,CQM,CQD);

        printf("    m   n      Re[Qmn(z)]    Im[Qmn(z)]    Re[Qmn'(z)]   Im[Qmn'(z)]\n");
        printf("   ------------------------------------------------------------------\n");
        
	for (J=0; J<=M; J++)
          printf(" %4d%4d  %14.6f%14.6f %14.6f%14.6f\n", J, N, CQM[J][N][0], CQM[J][N][1],
                 CQD[J][N][0], CQD[J][N][1]);
        printf("\n");

      }

// end of file mclqmn.cpp