APPENDIX 2 - SOURCE CODE - /usr/local/src/Math/legendre.c

/*
 *  legendre.c						    R S Bogart 91.08.22
 */
#include <math.h>
#include <stdio.h>
#include "Math.h"

#define PI4NORM (0.282094791773878143474039725780)

int check_Plm_args (x, nu, mu)
double x;
int nu, mu;
{
int mup;
/*
 *  Trap illegal index and argument values:
 *	nu < 0		return	 2
 *	|mu| > nu		 1
 *	|x| > 1.0		-1
 *			else	 0
 */
if (nu < 0)
  return (2);
mup = abs (mu);
if (mup > nu)
  return (1);
if (x < -1.0 || x > 1.0)
  return (-1);
return (0);
}

double legendre (cos_t, nu, mu, phase)
double cos_t;
int nu, mu, phase;
/*
 *  legendre
 *
 *  Generate normalized Associated Legendre Function value P_nu^mu (cos_t)
 *    using the stable upward recurrence on nu from P_mu^mu (Arfken 12.88,
 *    Magnus & Oberhettinger 4.3.3.2).  The normalization is such that
 *
 *      Y_l^m (cos_t, f) = P_l^m (cos_t) * exp (i*m*f),
 *    i.e.,
 *      P_l^m (normalized) = sqrt ((2*l + 1) / 4 pi) *
 *                           sqrt ((l - m)! / (l + m)!) * P_l^m
 *    (the customary (-)**m is subsumed in the choice of phase for the P_l^m).
 *
 *  Returns an IEEE quiet NaN for illegal argument values.
 *
 *  The phase is selected by the argument phase.  If (phase != 0) the
 *    Condon-Shortley phase is used, alternating sign for positive m.
 *    Otherwise the standard phase, alternating sign for negative m, is used.
 *
 */
{
double s, x;
double s0, s1, t0, t1, t2, f0, f1, f2;
double mu2;
int i, l, lp1, m, mup;

x = cos_t;
if (i = check_Plm_args (cos_t, nu, mu))
  if (i == 1)
    return (0.0);
  else
    return (A_Quiet_dNaN ());
mup = abs (mu);
/*
 *  Starting value P_m^m (x) = 1 / sqrt (4 Pi) *
 *                             sqrt [(2m + 1)!! / (2m)!!] * (sin_t ** m)
 *  (Note that P_m^(-m) = (-)**m * P_m^m)
 */
s = sqrt (1.0 - x*x);
t1 = PI4NORM;
for (m = 0; m < mup; m++)
  t1 *= s * sqrt ((m + m + 3.0) / (m + m + 2.0));		/*  P_mu^mu  */
/*
 *  Case l = |m|
 */
if (mup == nu)
{
  if (mu >= 0)
  {
    if (mup % 2 && phase)
      t1 *= -1;
  }
  else
  {
    if (mup % 2 && !phase)
      t1 *= -1;
  }
  return (t1);
}
/*
 *  Recurrence:
 *    P_(l+1)^m (cos_t) =
 *         sqrt(2l+3) / S_(l+1)^m *
 *         [sqrt(2l+1) * x * P_l^m - S_l^m * P_(l-1)^m / sqrt(2l-1)] ,
 *  where
 *    S_x^y == sqrt (x**2 - y**2)
 */
t2 = s1 = 0.0;
mu2 = mu * mu;
f2 = 1.0;
f1 = sqrt (mup + mup + 1.0);
for (l = mup, lp1 = mup + 1; l < nu; l++, lp1++)
{
  s0 = sqrt (lp1 * lp1 - mu2);
  f0 = sqrt (lp1 + lp1 + 1.0);
  t0 = f0 * (f1 * x * t1 - s1 * t2 / f2) / s0;
  t2 = t1;
  t1 = t0;
  s1 = s0;
  f2 = f1;
  f1 = f0;
}
if (mu >= 0)
{
  if (mup % 2 && phase)
    t0 *= -1;
}
else
{
  if (mup % 2 && !phase)
    t0 *= -1;
}
return (t0);
}

int legendre_by_l (Plm, cos_t, nu_max, mu, phase)
double *Plm;
double cos_t;
int nu_max, mu, phase;
/*
 *  legendre_by_l
 *
 *  Generates normalized Associated Legendre Function values P_nu^mu (cos_t)
 *    for all nu from mu to nu_max using the same algorithm as legendre ().
 *
 *  Values of P_nu^mu are returned in the array Plm, which must be properly
 *    malloc'd to size (nu_max - mu + 1) in the order Plm[0] = P_mu^mu,
 *    Plm[1] = P_(mu+1)^mu, ..., Plm[nu_max - mu] = P_(nu_max)^mu.  The
 *    function returns the value 0 unless the arguments are out of range,
 *    in which case it returns 1.  If the arguments permit reasonable sizing
 *    of the return array but are nonetheless out of range (e.g. |x| > 1),
 *    the array is filled with quiet NaN's.
 *
 */
{
double s, x;
double s0, s1, t0, t1, t2, f0, f1, f2;
double mu2;
int i, l, lp1, m, mup, sgn;

x = cos_t;
mup = abs (mu);
if (i = check_Plm_args (cos_t, nu_max, mu))
{
  if (i < 0)
    for (l = mup; l <= nu_max; l++)
      *Plm++ = A_Quiet_dNaN ();
  return (1);
}
/*
 *  Starting value P_m^m (x) = 1 / sqrt (4 Pi) *
 *                             sqrt [(2m + 1)!! / (2m)!!] * (sin_t ** m)
 *  (Note that P_m^(-m) = (-)**m * P_m^m)
 */
s = sqrt (1.0 - x*x);
t1 = PI4NORM;
for (m = 0; m < mup; m++)
  t1 *= s * sqrt ((m + m + 3.0) / (m + m + 2.0));		 /* P_mu^mu  */
if (mu >= 0)
  sgn = (mup % 2 && phase) ? -1 : 1;
else
  sgn = (mup % 2 && !phase) ? -1 : 1;
/*
 *  Recurrence:
 *    P_(l+1)^m (cos_t) =
 *         sqrt(2l+3) / S_(l+1)^m *
 *         [sqrt(2l+1) * x * P_l^m - S_l^m * P_(l-1)^m / sqrt(2l-1)] ,
 *  where
 *    S_x^y == sqrt (x**2 - y**2)
 */
t0 = t1;
t2 = s1 = 0.0;
mu2 = mu * mu;
f2 = 1.0;
f1 = sqrt (mup + mup + 1.0);
for (l = mup, lp1 = mup + 1; l < nu_max; l++, lp1++)
{
  *Plm++ = sgn * t0;
  s0 = sqrt (lp1 * lp1 - mu2);
  f0 = sqrt (lp1 + lp1 + 1.0);
  t0 = f0 * (f1 * x * t1 - s1 * t2 / f2) / s0;
  t2 = t1;
  t1 = t0;
  s1 = s0;
  f2 = f1;
  f1 = f0;
}
*Plm = sgn * t0;
return (0);
}

#define MIN_SIN	(0.0245)

int legendre_by_m (Plm, cos_t, nu, phase)
double *Plm;
double cos_t;
int nu, phase;
/*
 *  legendre_by_m
 *
 *  Generates normalized Associated Legendre Function values P_nu^mu (cos_t)
 *    for all mu from 0 to nu using the stable downward recurrence on mu
 *    from P_nu^nu (Arfken 12.87, M & O 4.3.3.1).  The normalization is
 *    such that
 *
 *      Y_l^m (cos_t, f) = P_l^m (cos_t) * exp (i*m*f),
 *    i.e.,
 *      P_l^m (normalized) = sqrt ((2*l + 1) / 4 pi) *
 *                           sqrt ((l - m)! / (l + m)!) * P_l^m
 *    (the customary (-)**m is subsumed in the choice of phase for the P_l^m).
 *
 *  Values of P_nu^mu are returned in the array Plm, which must be properly
 *    malloc'd to size mu in the order Plm[mu] = P_nu^mu.  The function
 *    returns the value 0 unless the arguments are out of range, in which
 *    it returns 1.  If the arguments permit reasonable sizing of the
 *    return array but are nonetheless out of range (e.g. |x| > 1), the
 *    array is filled with quiet NaN's.
 *
 *  The phase is selected by the argument phase as in legendre ().
 *
 *  N.B.  MIN_SIN is a cutoff value used exclusively by legendre_by_m as
 *    an escape from its standard recursion.  Because one of the recursion
 *    coefficients depends on cot(t) the recursion fails when x = cos(t)
 *    is close to 1 (sin(t) close to 0).  The point of failure depends
 *    not only on the value of l, but on the numerical accuracy of the
 *    math library square root function, presumably.  There is no point
 *    in making MIN_SIN smaller than about 1.5e-8, since that corresponds
 *    to the minimum non-zero value of 1-x in IEEE double precision.
 *    At that value, the recursion will work successfully on a NeXT for
 *    all values of l through at least 500, but fails on a a DEC MIPS
 *    system at l of 50.  With a value of 0.0245, the routine will work
 *    successfully at all l up to about 2000 on a NeXT, opting for direct
 *    calculation only when 1-x < 0.0003 (x > 0.9997).  It will fail
 *    on a DECsystem however for l >= 200 when x > 0.9997.  Choosing a
 *    value to allow valid calculation of high-order values on a 
 *    DECsystem is impractical. as it would require such a large
 *    fraction of calls to this routine to drop back to direct calculation
 *    that the time advantage would be lost.  As an example, a value
 *    of 0.13 would support work up to l = 1000, but would do case-by-case
 *    direct calculation for 13% of the cases.
 *   
 */
{
double s, x, cot_t;
double s0, s1, t0, t1, t2;
double llp1;
int i, l, m, sgn;

if (i = check_Plm_args (cos_t, nu, 0))
{
  if (i < 0)
    for (m = 0; m <= nu; m++)
      *Plm++ = A_Quiet_dNaN ();
  return (1);
}
x = cos_t;
/*
 *  Starting value P_l^l (x) = 1 / sqrt (4 Pi) *
 *                             sqrt [(2l + 1)!! / (2l)!!] * (sin_t ** l)
 *  Note that P_l^(-l) = (-)**l * P_l^l, so there is no reason to calculate it
 */
s = sqrt (1.0 - x*x);
t1 = PI4NORM;
Plm += nu;
if (s < MIN_SIN)
{
  l = nu;
  while (l--)
    *Plm-- = legendre (cos_t, nu, l, phase);
  *Plm = legendre (cos_t, nu, 0, phase);
  return (0);
}
cot_t = x / s;
for (l = 0; l < nu; l++)
  t1 *= s * sqrt ((l + l + 3.0) / (l + l + 2.0));		/*  P_nu^nu  */
sgn = (nu % 2 && phase) ? -1 : 1;
*Plm-- = sgn * t1;
/*
 *  Recurrence:
 *    P_l^(m-1)(cos_t) = [2m*cot_t*P_l^m - S_l^m * P_l^(m+1)] / S_l^(m-1),
 *  where
 *    S_x^y == sqrt [x*(x+1) - y*(y+1)]
 */
t2 = s1 = 0.0;
llp1 = nu * (nu + 1);
for (m = nu; m > 0; m--)
  {
    s0 = sqrt (llp1 - m*(m-1));
    t0 = (2*m*cot_t*t1 - s1*t2) / s0;
    t2 = t1;
    t1 = t0;
    s1 = s0;
    if (phase)
      sgn *= -1;
    *Plm-- = sgn * t0;
  }
return (0);
}

int legendre_by_x (Plm, cos_t, nx, nu, mu, phase)
double Plm[], cos_t[];
unsigned int nx;
int nu, mu, phase;
/*
 *  legendre_by_x
 *
 *  Generates normalized Associated Legendre Function values P_nu^mu (cos_t)
 *    for given nu and mu and for all values of cos_t in the array x, using
 *    the same algorithm as legendre ().
 *
 *  Values of P_nu^mu are returned in the array Plm, which must be properly
 *    malloc'd to size nx, in the same order as the values in the array x.
 *    The function returns the value 0 unless the arguments nu and mu are out
 *    of range, in which case it returns 1.  If nu and mu are in range but
 *    values of x are out of range (|cos_t| > 1), quiet NaN's are placed in
 *    the output array, but without additional signaling.
 *
 */
{
double s, x;
double f0, f1, f2, s0, s1, t0, t1, t2;
double mu2;
int i, l, lp1, m, mup, n;
double *vf0, *vf0_0, *vs0, *vs0_0, *vt1, *vt1_0;
double *Plm_0;

if (i = check_Plm_args (0.0, nu, mu))
{
  for (n = 0; n <= nx; n++)
    Plm[n] = A_Quiet_dNaN ();
  return (1);
}

mup = abs (mu);
mu2 = mu * mu;

vt1_0 = vt1 = (double *)calloc (mup, sizeof (double));
for (m = 0; m < mup; m++)
  *vt1++ = sqrt ((m + m + 3.0) / (m + m + 2.0));
  
vs0_0 = vs0 = (double *)calloc (nu - mup, sizeof (double));
vf0_0 = vf0 = (double *)calloc (nu - mup, sizeof (double));
for (l = mup, lp1 = mup+1; l < nu; l++, lp1++)
{
  *vs0++ = sqrt (lp1*lp1 - mu2);
  *vf0++ = sqrt (lp1 + lp1 + 1.0);
}

for (n = 0; n < nx; n++)
{
  x = cos_t[n];
  if (x < -1.0 || x > 1.0)
  {
    Plm[n] = A_Quiet_dNaN ();
    continue;
  }
/*
 *  Starting value P_m^m (x) = 1 / sqrt (4 Pi) *
 *                             sqrt [(2m + 1)!! / (2m)!!] * (sin_t ** m)
 *  (Note that P_m^(-m) = (-)**m * P_m^m)
 */
  s = sqrt (1.0 - x*x);
  t1 = PI4NORM;
  vt1 = vt1_0;
  for (m = 0; m < mup; m++)
    t1 *= s * (*vt1++);						/*  P_mu^mu  */
/*
 *  Case l = |m|
 */
  if (mup == nu)
  {
    if (mu >= 0)
    {
      if (mup%2 && phase)
        t1 *= -1;
    }
    else
    {
      if (mup%2 && !phase)
        t1 *= -1;
    }
    t0 = t1;
  }
  else
/*
 *  Recurrence:
 *    P_(l+1)^m (cos_t) = sqrt(2l+3) / S_(l+1)^m *
 *                        [sqrt(2l+1) * x * P_l^m - S_l^m * P_(l-1)^m /
 *                         sqrt(2l-1)] ,
 *  where
 *    S_x^y == sqrt (x**2 - y**2)
 */
  {
    t2 = 0.0;
    s1 = 0.0;
    f2 = 1.0;
    f1 = sqrt (mup + mup + 1.0);
    vs0 = vs0_0;
    vf0 = vf0_0;
    for (l = mup, lp1 = mup+1; l < nu; l++, lp1++)
    {
      s0 = *vs0++;
      f0 = *vf0++;
      t0 = f0 * (f1*x*t1 - s1*t2/f2) / s0;
      t2 = t1;
      t1 = t0;
      s1 = s0;
      f2 = f1;
      f1 = f0;
    }
  }
  if (mu >= 0)
  {
    if (mup%2 && phase)
      t0 *= -1;
  }
  else
  {
    if (mup%2 && !phase)
      t0 *= -1;
  }
  Plm[n] = t0;
}
free (vt1_0);
free (vs0_0);
free (vf0_0);
return (0);
}