/* 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