/* specfunc/legendre_source.c

 * 

 * Copyright (C) 2009-2013 Patrick Alken

 * 

 * This program is free software; you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation; either version 3 of the License, or (at

 * your option) any later version.

 * 

 * This program is distributed in the hope that it will be useful, but

 * WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

 * General Public License for more details.

 * 

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

 */

/* define various macros for functions below */

#define CONCAT2x(a,b)   a ## _ ## b 

#define CONCAT3x(a,b,c) a ## _ ## b ## _ ## c

#if defined(LEGENDRE)

#define FUNCTION(dir,name) CONCAT2x(dir,name)

#define OUTPUT result_array

#define OUTPUT_ARG double result_array[]

#elif defined(LEGENDRE_DERIV)

#define FUNCTION(dir,name) CONCAT3x(dir,deriv,name)

#define OUTPUT result_array, result_deriv_array

#define OUTPUT_ARG double result_array[], double result_deriv_array[]

#elif defined(LEGENDRE_DERIV_ALT)

#define FUNCTION(dir,name) CONCAT3x(dir,deriv_alt,name)

#define OUTPUT result_array, result_deriv_array

#define OUTPUT_ARG double result_array[], double result_deriv_array[]

#define LEGENDRE_DERIV

#elif defined(LEGENDRE_DERIV2)

#define FUNCTION(dir,name) CONCAT3x(dir,deriv2,name)

#define OUTPUT result_array, result_deriv_array, result_deriv2_array

#define OUTPUT_ARG double result_array[], double result_deriv_array[], double result_deriv2_array[]

#define LEGENDRE_DERIV

#elif defined(LEGENDRE_DERIV2_ALT)

#define FUNCTION(dir,name) CONCAT3x(dir,deriv2_alt,name)

#define OUTPUT result_array, result_deriv_array, result_deriv2_array

#define OUTPUT_ARG double result_array[], double result_deriv_array[], double result_deriv2_array[]

#define LEGENDRE_DERIV

#define LEGENDRE_DERIV2

#define LEGENDRE_DERIV_ALT

#endif

static int FUNCTION (legendre, array_schmidt_e)

(const size_t lmax, const double x, const double csphase, OUTPUT_ARG);

static int FUNCTION(legendre, array_none_e)

(const size_t lmax, const double x, const double csphase, OUTPUT_ARG);

/*

gsl_sf_legendre_array()

Inputs: norm                - normlization type

        lmax                - maximum degree

        x                   - input argument

        result_array        - (output) normalized P_{lm}

        result_deriv_array  - (output) normalized P'_{lm}

        result_deriv2_array - (output) normalized P''_{lm}

*/

int

FUNCTION (gsl_sf_legendre, array)

(const gsl_sf_legendre_t norm, const size_t lmax, const double x,

 OUTPUT_ARG)

{

  int s = FUNCTION (gsl_sf_legendre, array_e)(norm, lmax, x, 1.0, OUTPUT);

  return s;

}

/*

gsl_sf_legendre_array_e()

Inputs: norm                - normlization type

        lmax                - maximum degree

        x                   - input argument

        csphase             - Condon-Shortley phase

        result_array        - (output) normalized P_{lm}

        result_deriv_array  - (output) normalized P'_{lm}

        result_deriv2_array - (output) normalized P''_{lm}

*/

int

FUNCTION (gsl_sf_legendre, array_e)

(const gsl_sf_legendre_t norm, const size_t lmax, const double x,

 const double csphase, OUTPUT_ARG)

{

  int s;

  const size_t nlm = gsl_sf_legendre_nlm(lmax);

#if !defined(LEGENDRE_DERIV_ALT)

  size_t i;

#if defined(LEGENDRE_DERIV) || defined(LEGENDRE_DERIV2)

  const double u = sqrt((1.0 - x) * (1.0 + x));

  const double uinv = 1.0 / u;

#endif

#if defined(LEGENDRE_DERIV2)

  const double uinv2 = uinv * uinv;

#endif

#endif

  double fac1 = 0.0, fac2 = 0.0; /* normalization factors */

  if (norm == GSL_SF_LEGENDRE_NONE)

    {

      /* compute P_{lm}(x) */

      s = FUNCTION(legendre,array_none_e)(lmax, x, csphase, OUTPUT);

    }

  else

    {

      /* compute S_{lm}(x) */

      s = FUNCTION(legendre,array_schmidt_e)(lmax, x, csphase, OUTPUT);

    }

#if !defined(LEGENDRE_DERIV_ALT)

  /* scale derivative arrays to recover P'(x) and P''(x) */

  for (i = 0; i < nlm; ++i)

    {

#if defined(LEGENDRE_DERIV2)

      double dp = result_deriv_array[i];

      double d2p = result_deriv2_array[i];

      result_deriv2_array[i] = (d2p - x * uinv * dp) * uinv2;

#endif

#if defined(LEGENDRE_DERIV)

      result_deriv_array[i] *= -uinv;

#endif

    }

#endif

  /* apply scaling for requested normalization */

  if (norm == GSL_SF_LEGENDRE_SCHMIDT || norm == GSL_SF_LEGENDRE_NONE)

    {

      return s;

    }

  else if (norm == GSL_SF_LEGENDRE_SPHARM)

    {

      fac1 = 1.0 / sqrt(4.0 * M_PI);

      fac2 = 1.0 / sqrt(8.0 * M_PI);

    }

  else if (norm == GSL_SF_LEGENDRE_FULL)

    {

      fac1 = 1.0 / sqrt(2.0);

      fac2 = 1.0 / sqrt(4.0);

    }

  /*

   * common code for different normalizations

   * P_{l0} = fac1 * sqrt(2l + 1) * S_{l0}

   * P_{lm} = fac2 * sqrt(2l + 1) * S_{lm}, m > 0

   */

  {

    size_t l, m;

    size_t twoellp1 = 1; /* 2l + 1 */

    double *sqrts = &(result_array[nlm]);

    for (l = 0; l <= lmax; ++l)

      {

        result_array[gsl_sf_legendre_array_index(l, 0)] *=

          sqrts[twoellp1] * fac1;

#if defined(LEGENDRE_DERIV)

        result_deriv_array[gsl_sf_legendre_array_index(l, 0)] *=

          sqrts[twoellp1] * fac1;

#endif

#if defined(LEGENDRE_DERIV2)

        result_deriv2_array[gsl_sf_legendre_array_index(l, 0)] *=

          sqrts[twoellp1] * fac1;

#endif

        for (m = 1; m <= l; ++m)

          {

            result_array[gsl_sf_legendre_array_index(l, m)] *=

              sqrts[twoellp1] * fac2;

#if defined(LEGENDRE_DERIV)

            result_deriv_array[gsl_sf_legendre_array_index(l, m)] *=

              sqrts[twoellp1] * fac2;

#endif

#if defined(LEGENDRE_DERIV2)

            result_deriv2_array[gsl_sf_legendre_array_index(l, m)] *=

              sqrts[twoellp1] * fac2;

#endif

          }

        twoellp1 += 2;

      }

  }

  return s;

}

/*

legendre,array_schmidt_e()

  This routine computes Schmidt semi-normalized associated

Legendre polynomials and their first and second derivatives.

Inputs: lmax                - maximum order

        x                   - legendre argument in [-1,1]

        csphase             - -1.0 to include CS phase (-1)^m,

                               1.0 to not include

        result_array        - (output) where to store P_{lm}(x) values

        result_deriv_array  - (output) where to store

                              d/dtheta P_{lm}(x) values

        result_deriv2_array - (output) where to store

                              d^2/dtheta^2 P_{lm}(x) values

*/

static int

FUNCTION(legendre, array_schmidt_e)

(const size_t lmax, const double x, const double csphase, OUTPUT_ARG)

{

  if (x > 1.0 || x < -1.0)

    {

      GSL_ERROR("x is outside [-1,1]", GSL_EDOM);

    }

#if defined(LEGENDRE_DERIV) || defined(LEGENDRE_DERIV2)

  else if (fabs(x) == 1.0)

    {

      GSL_ERROR("x cannot equal 1 or -1 for derivative computation", GSL_EDOM);

    }

#endif

  else if (csphase != 1.0 && csphase != -1.0)

    {

      GSL_ERROR("csphase has invalid value", GSL_EDOM);

    }

  else

    {

      const double eps = 1.0e-280;

      const double u = sqrt((1.0 - x) * (1.0 + x)); /* sin(theta) */

#if defined(LEGENDRE_DERIV)

      const double uinv = 1.0 / u;

#endif

#if defined(LEGENDRE_DERIV2)

      const double uinv2 = 1.0 / u / u;

#endif

#if defined(LEGENDRE_DERIV) || defined(LEGENDRE_DERIV2)

      const double xbyu = x * uinv; /* x / u */

#endif

      size_t l, m;

      size_t k, idxmm;

      double plm, /* eps * S(l,m) / u^m */

             pmm; /* eps * S(m,m) / u^m */

      double rescalem;

      double pm1, /* S(l-1,m) */

             pm2; /* S(l-2,m) */

      size_t nlm = gsl_sf_legendre_nlm(lmax);

      double *sqrts = &(result_array[nlm]);

      /* precompute square root factors for recurrence */

      legendre_sqrts(lmax, sqrts);

      /* initial values S(0,0) and S(1,0) */

      pm2 = 1.0; /* S(0,0) */

      pm1 = x;   /* S(1,0) */

      result_array[0] = pm2;

#if defined(LEGENDRE_DERIV)

      result_deriv_array[0] = 0.0;

#endif

#if defined(LEGENDRE_DERIV2)

      result_deriv2_array[0] = 0.0;

#endif

      if (lmax == 0)

        return GSL_SUCCESS;

      result_array[1] = pm1;

#if defined(LEGENDRE_DERIV)

      result_deriv_array[1] = -u;

#endif

#if defined(LEGENDRE_DERIV2)

      result_deriv2_array[1] = -x;

#endif

      /* Compute S(l,0) for l=2..lmax, no scaling required */

      k = 1; /* idx(1,0) */

      for (l = 2; l <= lmax; ++l)

        {

          double linv = 1.0 / (double)l;

          k += l;  /* idx(l,m) = idx(l-1,m) + l */

          plm = (2.0 - linv) * x * pm1 - (1.0 - linv) * pm2;

          result_array[k] = plm;

#if defined(LEGENDRE_DERIV)

          result_deriv_array[k] = uinv * l * (x * plm - pm1);

#endif

#if defined(LEGENDRE_DERIV2)

          result_deriv2_array[k] = -(double) l * (l + 1.0) * plm -

                                   xbyu * result_deriv_array[k];

#endif

          pm2 = pm1;

          pm1 = plm;

        }

      /* Compute S(m,m), S(m+1,m) and S(l,m) */

      /*

       * pi_m = Prod_{i=2}^m sqrt[ (2m - 1) / (2m) ]

       * but pi_1 = 1.0, so initialize to sqrt(2) so that

       * the first m = 1 iteration of the loop will reset it

       * to 1.0. Starting with m = 2 it will begin accumulating

       * the correct terms.

       *

       * pmm = S(m,m) * eps / u^m = pi_m

       */

      pmm = sqrt(2.0) * eps;

      rescalem = 1.0 / eps;

      idxmm = 0; /* tracks idx(m,m), initialize to idx(0,0) */

      for (m = 1; m < lmax; ++m)

        {

          /* rescalem = u^m / eps */

          rescalem *= u;

          /*

           * compute:

           * S(m,m) = u * sqrt((2m - 1) / (2m)) S(m-1,m-1) = u^m * pi_m

           * d_t S(m,m) = m * x / u * S(m,m)

           */

          idxmm += m + 1; /* idx(m,m) = idx(m-1,m-1) + m + 1 */

          pmm *= csphase * sqrts[2 * m - 1] / sqrts[2 * m]; /* S(m,m) * eps / u^m */

          result_array[idxmm] = pmm * rescalem;

#if defined(LEGENDRE_DERIV)

          result_deriv_array[idxmm] = m * xbyu * result_array[idxmm];

#endif

#if defined(LEGENDRE_DERIV2)

          result_deriv2_array[idxmm] =

            m * (uinv2 * m - (m + 1.0)) * result_array[idxmm] -

            xbyu * result_deriv_array[idxmm];

#endif

          pm2 = pmm;

          /*

           * compute:

           * S(m+1,m) = sqrt(2 * m + 1) * x * S(m,m)

           * d_t S(m+1,m) = 1/u * ((m+1)*x*S(m+1,m) - sqrt(2*m+1)*S(m,m))

           */

          k = idxmm + m + 1; /* idx(m+1,m) = idx(m,m) + m + 1 */

          pm1 = x * sqrts[2 * m + 1] * pm2;

          result_array[k] = pm1 * rescalem;

#if defined(LEGENDRE_DERIV)

          result_deriv_array[k] =

            uinv * ((m + 1.0) * x * result_array[k] -

                    sqrts[2 * m + 1] * result_array[idxmm]);

#endif

#if defined(LEGENDRE_DERIV2)

          result_deriv2_array[k] =

            (m * m * uinv2 - (m + 1.0) * (m + 2.0)) * result_array[k] -

            xbyu * result_deriv_array[k];

#endif

          /* compute S(l,m) for l=m+2...lmax */

          for (l = m + 2; l <= lmax; ++l)

            {

              k += l; /* idx(l,m) = idx(l-1,m) + l */

              plm =

                (2*l - 1) / sqrts[l + m] / sqrts[l - m] * x * pm1 -

                sqrts[l - m - 1] * sqrts[l + m - 1] /

                sqrts[l + m] / sqrts[l - m] * pm2;

              result_array[k] = plm * rescalem;

#if defined(LEGENDRE_DERIV)

              result_deriv_array[k] =

                uinv * (l * x * result_array[k] -

                        sqrts[l + m] * sqrts[l - m] * result_array[k - l]);

#endif

#if defined(LEGENDRE_DERIV2)

              result_deriv2_array[k] =

                (m * m * uinv2 - l * (l + 1.0)) * result_array[k] -

                xbyu * result_deriv_array[k];

#endif

              pm2 = pm1;

              pm1 = plm;

            }

        } /* for (m = 1; m < lmax; ++m) */

      /* compute S(lmax,lmax) */

      rescalem *= u;

      idxmm += m + 1; /* idx(lmax,lmax) */

      pmm *= csphase * sqrts[2 * lmax - 1] / sqrts[2 * lmax];

      result_array[idxmm] = pmm * rescalem;

#if defined(LEGENDRE_DERIV)

      result_deriv_array[idxmm] = lmax * xbyu * result_array[idxmm];

#endif

#if defined(LEGENDRE_DERIV2)

      result_deriv2_array[idxmm] =

        lmax * (uinv2 * lmax - (lmax + 1.0)) * result_array[idxmm] -

        xbyu * result_deriv_array[idxmm];

#endif

      return GSL_SUCCESS;

    }

}

/*

legendre_array_none_e()

  This routine computes unnormalized associated Legendre polynomials

and their derivatives.

Inputs: lmax                - maximum order

        x                   - legendre argument in [-1,1]

        csphase             - -1.0 to include CS phase (-1)^m,

                               1.0 to not include

        result_array        - (output) where to store P_{lm}(x) values

        result_deriv_array  - (output) where to store

                              d/dtheta P_{lm}(x) values

        result_deriv2_array - (output) where to store

                              d^2/dtheta^2 P_{lm}(x) values

*/

static int

FUNCTION(legendre, array_none_e)

(const size_t lmax, const double x, const double csphase, OUTPUT_ARG)

{

  if (x > 1.0 || x < -1.0)

    {

      GSL_ERROR("x is outside [-1,1]", GSL_EDOM);

    }

#if defined(LEGENDRE_DERIV) || defined(LEGENDRE_DERIV2)

  else if (fabs(x) == 1.0)

    {

      GSL_ERROR("x cannot equal 1 or -1 for derivative computation", GSL_EDOM);

    }

#endif

  else if (csphase != 1.0 && csphase != -1.0)

    {

      GSL_ERROR("csphase has invalid value", GSL_EDOM);

    }

  else

    {

      const double u = sqrt((1.0 - x) * (1.0 + x)); /* sin(theta) */

#if defined(LEGENDRE_DERIV)

      const double uinv = 1.0 / u;

#endif

#if defined(LEGENDRE_DERIV2)

      const double uinv2 = 1.0 / u / u;

#endif

#if defined(LEGENDRE_DERIV) || defined(LEGENDRE_DERIV2)

      const double xbyu = x * uinv; /* x / u */

#endif

      size_t l, m;

      size_t k, idxmm;

      double plm, pmm;

      double pm1,    /* P(l-1,m) */

             pm2;    /* P(l-2,m) */

      double twomm1; /* 2*m - 1 */

      /* initial values P(0,0) and P(1,0) */

      pm2 = 1.0; /* P(0,0) */

      pm1 = x;   /* P(1,0) */

      result_array[0] = pm2;

#if defined(LEGENDRE_DERIV)

      result_deriv_array[0] = 0.0;

#endif

#if defined(LEGENDRE_DERIV2)

      result_deriv2_array[0] = 0.0;

#endif

      if (lmax == 0)

        return 0;

      result_array[1] = pm1;

#if defined(LEGENDRE_DERIV)

      result_deriv_array[1] = -u;

#endif

#if defined(LEGENDRE_DERIV2)

      result_deriv2_array[1] = -x;

#endif

      /* Compute P(l,0) */

      k = 1;

      for (l = 2; l <= lmax; ++l)

        {

          k += l;

          plm = ((2*l - 1) * x * pm1 - (l - 1) * pm2) / (double) l;

          result_array[k] = plm;

#if defined(LEGENDRE_DERIV)

          result_deriv_array[k] = -(double)l * (pm1 - x * plm) * uinv;

#endif

#if defined(LEGENDRE_DERIV2)

          result_deriv2_array[k] = -(double) l * (l + 1.0) * plm -

                                   xbyu * result_deriv_array[k];

#endif

          pm2 = pm1;

          pm1 = plm;

        }

      /* Compute P(m,m), P(m+1,m) and P(l,m) */

      pmm = 1.0;

      twomm1 = -1.0; /* 2 * m - 1 */

      idxmm = 0; /* tracks idx(m,m), initialize to idx(0,0) */

      for (m = 1; m <= lmax - 1; ++m)

        {

          /*

           * compute

           *

           * P(m,m) = u * (2m - 1) P(m-1,m-1)

           * and

           * dP(m,m)/dtheta = m * x * P(m,m) / u

           */

          idxmm += m + 1;

          twomm1 += 2.0;

          pmm *= csphase * u * twomm1;

          result_array[idxmm] = pmm;

#if defined(LEGENDRE_DERIV)

          result_deriv_array[idxmm] = m * xbyu * pmm;

#endif

#if defined(LEGENDRE_DERIV2)

          result_deriv2_array[idxmm] =

            m * (uinv2 * m - (m + 1.0)) * result_array[idxmm] -

            xbyu * result_deriv_array[idxmm];

#endif

          pm2 = pmm;

          /*

           * compute

           *

           * P(m+1,m) = (2 * m + 1) * x * P(m,m)

           * and

           * dP(m+1,m)/dt = -[(2*m + 1) * P(m,m) - (m+1) * x * P(m+1,m)]/u

           */

          k = idxmm + m + 1;

          pm1 = x * pmm * (2*m + 1);

          result_array[k] = pm1;

#if defined(LEGENDRE_DERIV)

          result_deriv_array[k] = -uinv * ((2*m + 1) * pmm - (m + 1) * x * pm1);

#endif

#if defined(LEGENDRE_DERIV2)

          result_deriv2_array[k] =

            (m * m * uinv2 - (m + 1.0) * (m + 2.0)) * result_array[k] -

            xbyu * result_deriv_array[k];

#endif

          /* compute P(l,m) */

          for (l = m + 2; l <= lmax; ++l)

            {

              k += l;

              plm = ((2*l - 1) * x * pm1 - (l + m - 1) * pm2) /

                    (double) (l - m);

              result_array[k] = plm;

#if defined(LEGENDRE_DERIV)

              result_deriv_array[k] = -uinv * ((l + m) * pm1 - l * x * plm);

#endif

#if defined(LEGENDRE_DERIV2)

              result_deriv2_array[k] =

                (m * m * uinv2 - l * (l + 1.0)) * result_array[k] -

                xbyu * result_deriv_array[k];

#endif

              pm2 = pm1;

              pm1 = plm;

            }

        } /* for (m = 1; m <= lmax - 1; ++m) */

      /* compute P(lmax,lmax) */

      idxmm += m + 1;

      twomm1 += 2.0;

      pmm *= csphase * u * twomm1;

      result_array[idxmm] = pmm;

#if defined(LEGENDRE_DERIV)

      result_deriv_array[idxmm] = lmax * x * pmm * uinv;

#endif

#if defined(LEGENDRE_DERIV2)

      result_deriv2_array[idxmm] =

        lmax * (uinv2 * lmax - (lmax + 1.0)) * result_array[idxmm] -

        xbyu * result_deriv_array[idxmm];

#endif

      return GSL_SUCCESS;

    }

} /* legendre_array_none_e() */

#undef FUNCTION

#undef CONCAT2x

#undef CONCAT3x

#undef OUTPUT

#undef OUTPUT_ARG

#undef LEGENDRE_DERIV

#undef LEGENDRE_DERIV2

#undef LEGENDRE_DERIV_ALT