|
Packit |
67cb25 |
/* specfunc/legendre_poly.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002 Gerard Jungman
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Author: G. Jungman */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_bessel.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_exp.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_gamma.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_log.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_pow_int.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_legendre.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "error.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Calculate P_m^m(x) from the analytic result:
|
|
Packit |
67cb25 |
* P_m^m(x) = (-1)^m (2m-1)!! (1-x^2)^(m/2) , m > 0
|
|
Packit |
67cb25 |
* = 1 , m = 0
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static double legendre_Pmm(int m, double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(m == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return 1.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double p_mm = 1.0;
|
|
Packit |
67cb25 |
double root_factor = sqrt(1.0-x)*sqrt(1.0+x);
|
|
Packit |
67cb25 |
double fact_coeff = 1.0;
|
|
Packit |
67cb25 |
int i;
|
|
Packit |
67cb25 |
for(i=1; i<=m; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
p_mm *= -fact_coeff * root_factor;
|
|
Packit |
67cb25 |
fact_coeff += 2.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
return p_mm;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_legendre_P1_e(double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result->val = x;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_legendre_P2_e(double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result->val = 0.5*(3.0*x*x - 1.0);
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * (fabs(3.0*x*x) + 1.0);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_legendre_P3_e(double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result->val = 0.5*x*(5.0*x*x - 3.0);
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * (fabs(result->val) + 0.5 * fabs(x) * (fabs(5.0*x*x) + 3.0));
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_legendre_Pl_e(const int l, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(l < 0 || x < -1.0 || x > 1.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(l == 0) {
|
|
Packit |
67cb25 |
result->val = 1.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(l == 1) {
|
|
Packit |
67cb25 |
result->val = x;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(l == 2) {
|
|
Packit |
67cb25 |
result->val = 0.5 * (3.0*x*x - 1.0);
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * (fabs(3.0*x*x) + 1.0);
|
|
Packit |
67cb25 |
/*result->err = 3.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
removed this old bogus estimate [GJ]
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x == 1.0) {
|
|
Packit |
67cb25 |
result->val = 1.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x == -1.0) {
|
|
Packit |
67cb25 |
result->val = ( GSL_IS_ODD(l) ? -1.0 : 1.0 );
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(l < 100000) {
|
|
Packit |
67cb25 |
/* upward recurrence: l P_l = (2l-1) z P_{l-1} - (l-1) P_{l-2} */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double p_ellm2 = 1.0; /* P_0(x) */
|
|
Packit |
67cb25 |
double p_ellm1 = x; /* P_1(x) */
|
|
Packit |
67cb25 |
double p_ell = p_ellm1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double e_ellm2 = GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
double e_ellm1 = fabs(x)*GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
double e_ell = e_ellm1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int ell;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(ell=2; ell <= l; ell++){
|
|
Packit |
67cb25 |
p_ell = (x*(2*ell-1)*p_ellm1 - (ell-1)*p_ellm2) / ell;
|
|
Packit |
67cb25 |
p_ellm2 = p_ellm1;
|
|
Packit |
67cb25 |
p_ellm1 = p_ell;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
e_ell = 0.5*(fabs(x)*(2*ell-1.0) * e_ellm1 + (ell-1.0)*e_ellm2)/ell;
|
|
Packit |
67cb25 |
e_ellm2 = e_ellm1;
|
|
Packit |
67cb25 |
e_ellm1 = e_ell;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = p_ell;
|
|
Packit |
67cb25 |
result->err = e_ell + l*fabs(p_ell)*GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* Asymptotic expansion.
|
|
Packit |
67cb25 |
* [Olver, p. 473]
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
double u = l + 0.5;
|
|
Packit |
67cb25 |
double th = acos(x);
|
|
Packit |
67cb25 |
gsl_sf_result J0;
|
|
Packit |
67cb25 |
gsl_sf_result Jm1;
|
|
Packit |
67cb25 |
int stat_J0 = gsl_sf_bessel_J0_e(u*th, &J0;;
|
|
Packit |
67cb25 |
int stat_Jm1 = gsl_sf_bessel_Jn_e(-1, u*th, &Jm1);
|
|
Packit |
67cb25 |
double pre;
|
|
Packit |
67cb25 |
double B00;
|
|
Packit |
67cb25 |
double c1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* B00 = 1/8 (1 - th cot(th) / th^2
|
|
Packit |
67cb25 |
* pre = sqrt(th/sin(th))
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
if(th < GSL_ROOT4_DBL_EPSILON) {
|
|
Packit |
67cb25 |
B00 = (1.0 + th*th/15.0)/24.0;
|
|
Packit |
67cb25 |
pre = 1.0 + th*th/12.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
double sin_th = sqrt(1.0 - x*x);
|
|
Packit |
67cb25 |
double cot_th = x / sin_th;
|
|
Packit |
67cb25 |
B00 = 1.0/8.0 * (1.0 - th * cot_th) / (th*th);
|
|
Packit |
67cb25 |
pre = sqrt(th/sin_th);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
c1 = th/u * B00;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = pre * (J0.val + c1 * Jm1.val);
|
|
Packit |
67cb25 |
result->err = pre * (J0.err + fabs(c1) * Jm1.err);
|
|
Packit |
67cb25 |
result->err += GSL_SQRT_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_J0, stat_Jm1);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_legendre_Pl_array(const int lmax, const double x, double * result_array)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result_array) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(lmax < 0 || x < -1.0 || x > 1.0) {
|
|
Packit |
67cb25 |
GSL_ERROR ("domain error", GSL_EDOM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(lmax == 0) {
|
|
Packit |
67cb25 |
result_array[0] = 1.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(lmax == 1) {
|
|
Packit |
67cb25 |
result_array[0] = 1.0;
|
|
Packit |
67cb25 |
result_array[1] = x;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* upward recurrence: l P_l = (2l-1) z P_{l-1} - (l-1) P_{l-2} */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double p_ellm2 = 1.0; /* P_0(x) */
|
|
Packit |
67cb25 |
double p_ellm1 = x; /* P_1(x) */
|
|
Packit |
67cb25 |
double p_ell = p_ellm1;
|
|
Packit |
67cb25 |
int ell;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result_array[0] = 1.0;
|
|
Packit |
67cb25 |
result_array[1] = x;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(ell=2; ell <= lmax; ell++){
|
|
Packit |
67cb25 |
p_ell = (x*(2*ell-1)*p_ellm1 - (ell-1)*p_ellm2) / ell;
|
|
Packit |
67cb25 |
p_ellm2 = p_ellm1;
|
|
Packit |
67cb25 |
p_ellm1 = p_ell;
|
|
Packit |
67cb25 |
result_array[ell] = p_ell;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_legendre_Pl_deriv_array(const int lmax, const double x, double * result_array, double * result_deriv_array)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int stat_array = gsl_sf_legendre_Pl_array(lmax, x, result_array);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(lmax >= 0) result_deriv_array[0] = 0.0;
|
|
Packit |
67cb25 |
if(lmax >= 1) result_deriv_array[1] = 1.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(stat_array == GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int ell;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(x - 1.0)*(lmax+1.0)*(lmax+1.0) < GSL_SQRT_DBL_EPSILON)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* x is near 1 */
|
|
Packit |
67cb25 |
for(ell = 2; ell <= lmax; ell++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double pre = 0.5 * ell * (ell+1.0);
|
|
Packit |
67cb25 |
result_deriv_array[ell] = pre * (1.0 - 0.25 * (1.0-x) * (ell+2.0)*(ell-1.0));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(fabs(x + 1.0)*(lmax+1.0)*(lmax+1.0) < GSL_SQRT_DBL_EPSILON)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* x is near -1 */
|
|
Packit |
67cb25 |
for(ell = 2; ell <= lmax; ell++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double sgn = ( GSL_IS_ODD(ell) ? 1.0 : -1.0 ); /* derivative is odd in x for even ell */
|
|
Packit |
67cb25 |
const double pre = sgn * 0.5 * ell * (ell+1.0);
|
|
Packit |
67cb25 |
result_deriv_array[ell] = pre * (1.0 - 0.25 * (1.0+x) * (ell+2.0)*(ell-1.0));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double diff_a = 1.0 + x;
|
|
Packit |
67cb25 |
const double diff_b = 1.0 - x;
|
|
Packit |
67cb25 |
for(ell = 2; ell <= lmax; ell++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result_deriv_array[ell] = - ell * (x * result_array[ell] - result_array[ell-1]) / (diff_a * diff_b);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return stat_array;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_legendre_Plm_e(const int l, const int m, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* If l is large and m is large, then we have to worry
|
|
Packit |
67cb25 |
* about overflow. Calculate an approximate exponent which
|
|
Packit |
67cb25 |
* measures the normalization of this thing.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
const double dif = l-m;
|
|
Packit |
67cb25 |
const double sum = l+m;
|
|
Packit |
67cb25 |
const double t_d = ( dif == 0.0 ? 0.0 : 0.5 * dif * (log(dif)-1.0) );
|
|
Packit |
67cb25 |
const double t_s = ( dif == 0.0 ? 0.0 : 0.5 * sum * (log(sum)-1.0) );
|
|
Packit |
67cb25 |
const double exp_check = 0.5 * log(2.0*l+1.0) + t_d - t_s;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(m < 0 || l < m || x < -1.0 || x > 1.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(exp_check < GSL_LOG_DBL_MIN + 10.0){
|
|
Packit |
67cb25 |
/* Bail out. */
|
|
Packit |
67cb25 |
OVERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* Account for the error due to the
|
|
Packit |
67cb25 |
* representation of 1-x.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
const double err_amp = 1.0 / (GSL_DBL_EPSILON + fabs(1.0-fabs(x)));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* P_m^m(x) and P_{m+1}^m(x) */
|
|
Packit |
67cb25 |
double p_mm = legendre_Pmm(m, x);
|
|
Packit |
67cb25 |
double p_mmp1 = x * (2*m + 1) * p_mm;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(l == m){
|
|
Packit |
67cb25 |
result->val = p_mm;
|
|
Packit |
67cb25 |
result->err = err_amp * 2.0 * GSL_DBL_EPSILON * fabs(p_mm);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(l == m + 1) {
|
|
Packit |
67cb25 |
result->val = p_mmp1;
|
|
Packit |
67cb25 |
result->err = err_amp * 2.0 * GSL_DBL_EPSILON * fabs(p_mmp1);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else{
|
|
Packit |
67cb25 |
/* upward recurrence: (l-m) P(l,m) = (2l-1) z P(l-1,m) - (l+m-1) P(l-2,m)
|
|
Packit |
67cb25 |
* start at P(m,m), P(m+1,m)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double p_ellm2 = p_mm;
|
|
Packit |
67cb25 |
double p_ellm1 = p_mmp1;
|
|
Packit |
67cb25 |
double p_ell = 0.0;
|
|
Packit |
67cb25 |
int ell;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(ell=m+2; ell <= l; ell++){
|
|
Packit |
67cb25 |
p_ell = (x*(2*ell-1)*p_ellm1 - (ell+m-1)*p_ellm2) / (ell-m);
|
|
Packit |
67cb25 |
p_ellm2 = p_ellm1;
|
|
Packit |
67cb25 |
p_ellm1 = p_ell;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = p_ell;
|
|
Packit |
67cb25 |
result->err = err_amp * (0.5*(l-m) + 1.0) * GSL_DBL_EPSILON * fabs(p_ell);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_legendre_sphPlm_e(const int l, int m, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(m < 0 || l < m || x < -1.0 || x > 1.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(m == 0) {
|
|
Packit |
67cb25 |
gsl_sf_result P;
|
|
Packit |
67cb25 |
int stat_P = gsl_sf_legendre_Pl_e(l, x, &P);
|
|
Packit |
67cb25 |
double pre = sqrt((2.0*l + 1.0)/(4.0*M_PI));
|
|
Packit |
67cb25 |
result->val = pre * P.val;
|
|
Packit |
67cb25 |
result->err = pre * P.err;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_P;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x == 1.0 || x == -1.0) {
|
|
Packit |
67cb25 |
/* m > 0 here */
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* m > 0 and |x| < 1 here */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Starting value for recursion.
|
|
Packit |
67cb25 |
* Y_m^m(x) = sqrt( (2m+1)/(4pi m) gamma(m+1/2)/gamma(m) ) (-1)^m (1-x^2)^(m/2) / pi^(1/4)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result lncirc;
|
|
Packit |
67cb25 |
gsl_sf_result lnpoch;
|
|
Packit |
67cb25 |
double lnpre_val;
|
|
Packit |
67cb25 |
double lnpre_err;
|
|
Packit |
67cb25 |
gsl_sf_result ex_pre;
|
|
Packit |
67cb25 |
double sr;
|
|
Packit |
67cb25 |
const double sgn = ( GSL_IS_ODD(m) ? -1.0 : 1.0);
|
|
Packit |
67cb25 |
const double y_mmp1_factor = x * sqrt(2.0*m + 3.0);
|
|
Packit |
67cb25 |
double y_mm, y_mm_err;
|
|
Packit |
67cb25 |
double y_mmp1, y_mmp1_err;
|
|
Packit |
67cb25 |
gsl_sf_log_1plusx_e(-x*x, &lncirc);
|
|
Packit |
67cb25 |
gsl_sf_lnpoch_e(m, 0.5, &lnpoch); /* Gamma(m+1/2)/Gamma(m) */
|
|
Packit |
67cb25 |
lnpre_val = -0.25*M_LNPI + 0.5 * (lnpoch.val + m*lncirc.val);
|
|
Packit |
67cb25 |
lnpre_err = 0.25*M_LNPI*GSL_DBL_EPSILON + 0.5 * (lnpoch.err + fabs(m)*lncirc.err);
|
|
Packit |
67cb25 |
/* Compute exp(ln_pre) with error term, avoiding call to gsl_sf_exp_err BJG */
|
|
Packit |
67cb25 |
ex_pre.val = exp(lnpre_val);
|
|
Packit |
67cb25 |
ex_pre.err = 2.0*(sinh(lnpre_err) + GSL_DBL_EPSILON)*ex_pre.val;
|
|
Packit |
67cb25 |
sr = sqrt((2.0+1.0/m)/(4.0*M_PI));
|
|
Packit |
67cb25 |
y_mm = sgn * sr * ex_pre.val;
|
|
Packit |
67cb25 |
y_mm_err = 2.0 * GSL_DBL_EPSILON * fabs(y_mm) + sr * ex_pre.err;
|
|
Packit |
67cb25 |
y_mm_err *= 1.0 + 1.0/(GSL_DBL_EPSILON + fabs(1.0-x));
|
|
Packit |
67cb25 |
y_mmp1 = y_mmp1_factor * y_mm;
|
|
Packit |
67cb25 |
y_mmp1_err=fabs(y_mmp1_factor) * y_mm_err;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(l == m){
|
|
Packit |
67cb25 |
result->val = y_mm;
|
|
Packit |
67cb25 |
result->err = y_mm_err;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(y_mm);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(l == m + 1) {
|
|
Packit |
67cb25 |
result->val = y_mmp1;
|
|
Packit |
67cb25 |
result->err = y_mmp1_err;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(y_mmp1);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else{
|
|
Packit |
67cb25 |
double y_ell = 0.0;
|
|
Packit |
67cb25 |
double y_ell_err = 0.0;
|
|
Packit |
67cb25 |
int ell;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute Y_l^m, l > m+1, upward recursion on l. */
|
|
Packit |
67cb25 |
for(ell=m+2; ell <= l; ell++){
|
|
Packit |
67cb25 |
const double rat1 = (double)(ell-m)/(double)(ell+m);
|
|
Packit |
67cb25 |
const double rat2 = (ell-m-1.0)/(ell+m-1.0);
|
|
Packit |
67cb25 |
const double factor1 = sqrt(rat1*(2.0*ell+1.0)*(2.0*ell-1.0));
|
|
Packit |
67cb25 |
const double factor2 = sqrt(rat1*rat2*(2.0*ell+1.0)/(2.0*ell-3.0));
|
|
Packit |
67cb25 |
y_ell = (x*y_mmp1*factor1 - (ell+m-1.0)*y_mm*factor2) / (ell-m);
|
|
Packit |
67cb25 |
y_mm = y_mmp1;
|
|
Packit |
67cb25 |
y_mmp1 = y_ell;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
y_ell_err = 0.5*(fabs(x*factor1)*y_mmp1_err + fabs((ell+m-1.0)*factor2)*y_mm_err) / fabs(ell-m);
|
|
Packit |
67cb25 |
y_mm_err = y_mmp1_err;
|
|
Packit |
67cb25 |
y_mmp1_err = y_ell_err;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = y_ell;
|
|
Packit |
67cb25 |
result->err = y_ell_err + (0.5*(l-m) + 1.0) * GSL_DBL_EPSILON * fabs(y_ell);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#ifndef GSL_DISABLE_DEPRECATED
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_legendre_Plm_array(const int lmax, const int m, const double x, double * result_array)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* If l is large and m is large, then we have to worry
|
|
Packit |
67cb25 |
* about overflow. Calculate an approximate exponent which
|
|
Packit |
67cb25 |
* measures the normalization of this thing.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
const double dif = lmax-m;
|
|
Packit |
67cb25 |
const double sum = lmax+m;
|
|
Packit |
67cb25 |
const double t_d = ( dif == 0.0 ? 0.0 : 0.5 * dif * (log(dif)-1.0) );
|
|
Packit |
67cb25 |
const double t_s = ( dif == 0.0 ? 0.0 : 0.5 * sum * (log(sum)-1.0) );
|
|
Packit |
67cb25 |
const double exp_check = 0.5 * log(2.0*lmax+1.0) + t_d - t_s;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result_array) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(m < 0 || lmax < m || x < -1.0 || x > 1.0) {
|
|
Packit |
67cb25 |
GSL_ERROR ("domain error", GSL_EDOM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(m > 0 && (x == 1.0 || x == -1.0)) {
|
|
Packit |
67cb25 |
int ell;
|
|
Packit |
67cb25 |
for(ell=m; ell<=lmax; ell++) result_array[ell-m] = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(exp_check < GSL_LOG_DBL_MIN + 10.0){
|
|
Packit |
67cb25 |
/* Bail out. */
|
|
Packit |
67cb25 |
GSL_ERROR ("overflow", GSL_EOVRFLW);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
double p_mm = legendre_Pmm(m, x);
|
|
Packit |
67cb25 |
double p_mmp1 = x * (2.0*m + 1.0) * p_mm;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(lmax == m){
|
|
Packit |
67cb25 |
result_array[0] = p_mm;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(lmax == m + 1) {
|
|
Packit |
67cb25 |
result_array[0] = p_mm;
|
|
Packit |
67cb25 |
result_array[1] = p_mmp1;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
double p_ellm2 = p_mm;
|
|
Packit |
67cb25 |
double p_ellm1 = p_mmp1;
|
|
Packit |
67cb25 |
double p_ell = 0.0;
|
|
Packit |
67cb25 |
int ell;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result_array[0] = p_mm;
|
|
Packit |
67cb25 |
result_array[1] = p_mmp1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(ell=m+2; ell <= lmax; ell++){
|
|
Packit |
67cb25 |
p_ell = (x*(2.0*ell-1.0)*p_ellm1 - (ell+m-1)*p_ellm2) / (ell-m);
|
|
Packit |
67cb25 |
p_ellm2 = p_ellm1;
|
|
Packit |
67cb25 |
p_ellm1 = p_ell;
|
|
Packit |
67cb25 |
result_array[ell-m] = p_ell;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_legendre_Plm_deriv_array(
|
|
Packit |
67cb25 |
const int lmax, const int m, const double x,
|
|
Packit |
67cb25 |
double * result_array,
|
|
Packit |
67cb25 |
double * result_deriv_array)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(m < 0 || m > lmax)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("m < 0 or m > lmax", GSL_EDOM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(m == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* It is better to do m=0 this way, so we can more easily
|
|
Packit |
67cb25 |
* trap the divergent case which can occur when m == 1.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return gsl_sf_legendre_Pl_deriv_array(lmax, x, result_array, result_deriv_array);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int stat_array = gsl_sf_legendre_Plm_array(lmax, m, x, result_array);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(stat_array == GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int ell;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(m == 1 && (1.0 - fabs(x) < GSL_DBL_EPSILON))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* This divergence is real and comes from the cusp-like
|
|
Packit |
67cb25 |
* behaviour for m = 1. For example, P[1,1] = - Sqrt[1-x^2].
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
GSL_ERROR("divergence near |x| = 1.0 since m = 1", GSL_EOVRFLW);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(m == 2 && (1.0 - fabs(x) < GSL_DBL_EPSILON))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* m = 2 gives a finite nonzero result for |x| near 1 */
|
|
Packit |
67cb25 |
if(fabs(x - 1.0) < GSL_DBL_EPSILON)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for(ell = m; ell <= lmax; ell++) result_deriv_array[ell-m] = -0.25 * x * (ell - 1.0)*ell*(ell+1.0)*(ell+2.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(fabs(x + 1.0) < GSL_DBL_EPSILON)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for(ell = m; ell <= lmax; ell++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double sgn = ( GSL_IS_ODD(ell) ? 1.0 : -1.0 );
|
|
Packit |
67cb25 |
result_deriv_array[ell-m] = -0.25 * sgn * x * (ell - 1.0)*ell*(ell+1.0)*(ell+2.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* m > 2 is easier to deal with since the endpoints always vanish */
|
|
Packit |
67cb25 |
if(1.0 - fabs(x) < GSL_DBL_EPSILON)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for(ell = m; ell <= lmax; ell++) result_deriv_array[ell-m] = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double diff_a = 1.0 + x;
|
|
Packit |
67cb25 |
const double diff_b = 1.0 - x;
|
|
Packit |
67cb25 |
result_deriv_array[0] = - m * x / (diff_a * diff_b) * result_array[0];
|
|
Packit |
67cb25 |
if(lmax-m >= 1) result_deriv_array[1] = (2.0 * m + 1.0) * (x * result_deriv_array[0] + result_array[0]);
|
|
Packit |
67cb25 |
for(ell = m+2; ell <= lmax; ell++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result_deriv_array[ell-m] = - (ell * x * result_array[ell-m] - (ell+m) * result_array[ell-1-m]) / (diff_a * diff_b);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return stat_array;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_legendre_sphPlm_array(const int lmax, int m, const double x, double * result_array)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result_array) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(m < 0 || lmax < m || x < -1.0 || x > 1.0) {
|
|
Packit |
67cb25 |
GSL_ERROR ("error", GSL_EDOM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(m > 0 && (x == 1.0 || x == -1.0)) {
|
|
Packit |
67cb25 |
int ell;
|
|
Packit |
67cb25 |
for(ell=m; ell<=lmax; ell++) result_array[ell-m] = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
double y_mm;
|
|
Packit |
67cb25 |
double y_mmp1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(m == 0) {
|
|
Packit |
67cb25 |
y_mm = 0.5/M_SQRTPI; /* Y00 = 1/sqrt(4pi) */
|
|
Packit |
67cb25 |
y_mmp1 = x * M_SQRT3 * y_mm;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* |x| < 1 here */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_sf_result lncirc;
|
|
Packit |
67cb25 |
gsl_sf_result lnpoch;
|
|
Packit |
67cb25 |
double lnpre;
|
|
Packit |
67cb25 |
const double sgn = ( GSL_IS_ODD(m) ? -1.0 : 1.0);
|
|
Packit |
67cb25 |
gsl_sf_log_1plusx_e(-x*x, &lncirc);
|
|
Packit |
67cb25 |
gsl_sf_lnpoch_e(m, 0.5, &lnpoch); /* Gamma(m+1/2)/Gamma(m) */
|
|
Packit |
67cb25 |
lnpre = -0.25*M_LNPI + 0.5 * (lnpoch.val + m*lncirc.val);
|
|
Packit |
67cb25 |
y_mm = sqrt((2.0+1.0/m)/(4.0*M_PI)) * sgn * exp(lnpre);
|
|
Packit |
67cb25 |
y_mmp1 = x * sqrt(2.0*m + 3.0) * y_mm;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(lmax == m){
|
|
Packit |
67cb25 |
result_array[0] = y_mm;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(lmax == m + 1) {
|
|
Packit |
67cb25 |
result_array[0] = y_mm;
|
|
Packit |
67cb25 |
result_array[1] = y_mmp1;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else{
|
|
Packit |
67cb25 |
double y_ell;
|
|
Packit |
67cb25 |
int ell;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result_array[0] = y_mm;
|
|
Packit |
67cb25 |
result_array[1] = y_mmp1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute Y_l^m, l > m+1, upward recursion on l. */
|
|
Packit |
67cb25 |
for(ell=m+2; ell <= lmax; ell++){
|
|
Packit |
67cb25 |
const double rat1 = (double)(ell-m)/(double)(ell+m);
|
|
Packit |
67cb25 |
const double rat2 = (ell-m-1.0)/(ell+m-1.0);
|
|
Packit |
67cb25 |
const double factor1 = sqrt(rat1*(2*ell+1)*(2*ell-1));
|
|
Packit |
67cb25 |
const double factor2 = sqrt(rat1*rat2*(2*ell+1)/(2*ell-3));
|
|
Packit |
67cb25 |
y_ell = (x*y_mmp1*factor1 - (ell+m-1)*y_mm*factor2) / (ell-m);
|
|
Packit |
67cb25 |
y_mm = y_mmp1;
|
|
Packit |
67cb25 |
y_mmp1 = y_ell;
|
|
Packit |
67cb25 |
result_array[ell-m] = y_ell;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_legendre_sphPlm_deriv_array(
|
|
Packit |
67cb25 |
const int lmax, const int m, const double x,
|
|
Packit |
67cb25 |
double * result_array,
|
|
Packit |
67cb25 |
double * result_deriv_array)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(m < 0 || lmax < m || x < -1.0 || x > 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("domain", GSL_EDOM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(m == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* m = 0 is easy to trap */
|
|
Packit |
67cb25 |
const int stat_array = gsl_sf_legendre_Pl_deriv_array(lmax, x, result_array, result_deriv_array);
|
|
Packit |
67cb25 |
int ell;
|
|
Packit |
67cb25 |
for(ell = 0; ell <= lmax; ell++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double prefactor = sqrt((2.0 * ell + 1.0)/(4.0*M_PI));
|
|
Packit |
67cb25 |
result_array[ell] *= prefactor;
|
|
Packit |
67cb25 |
result_deriv_array[ell] *= prefactor;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
return stat_array;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(m == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Trapping m = 1 is necessary because of the possible divergence.
|
|
Packit |
67cb25 |
* Recall that this divergence is handled properly in ..._Plm_deriv_array(),
|
|
Packit |
67cb25 |
* and the scaling factor is not large for small m, so we just scale.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
const int stat_array = gsl_sf_legendre_Plm_deriv_array(lmax, m, x, result_array, result_deriv_array);
|
|
Packit |
67cb25 |
int ell;
|
|
Packit |
67cb25 |
for(ell = 1; ell <= lmax; ell++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double prefactor = sqrt((2.0 * ell + 1.0)/(ell + 1.0) / (4.0*M_PI*ell));
|
|
Packit |
67cb25 |
result_array[ell-1] *= prefactor;
|
|
Packit |
67cb25 |
result_deriv_array[ell-1] *= prefactor;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
return stat_array;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* as for the derivative of P_lm, everything is regular for m >= 2 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int stat_array = gsl_sf_legendre_sphPlm_array(lmax, m, x, result_array);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(stat_array == GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int ell;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(1.0 - fabs(x) < GSL_DBL_EPSILON)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for(ell = m; ell <= lmax; ell++) result_deriv_array[ell-m] = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double diff_a = 1.0 + x;
|
|
Packit |
67cb25 |
const double diff_b = 1.0 - x;
|
|
Packit |
67cb25 |
result_deriv_array[0] = - m * x / (diff_a * diff_b) * result_array[0];
|
|
Packit |
67cb25 |
if(lmax-m >= 1) result_deriv_array[1] = sqrt(2.0 * m + 3.0) * (x * result_deriv_array[0] + result_array[0]);
|
|
Packit |
67cb25 |
for(ell = m+2; ell <= lmax; ell++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double c1 = sqrt(((2.0*ell+1.0)/(2.0*ell-1.0)) * ((double)(ell-m)/(double)(ell+m)));
|
|
Packit |
67cb25 |
result_deriv_array[ell-m] = - (ell * x * result_array[ell-m] - c1 * (ell+m) * result_array[ell-1-m]) / (diff_a * diff_b);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return stat_array;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_legendre_array_size(const int lmax, const int m)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return lmax-m+1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#endif /* !GSL_DISABLE_DEPRECATED */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "eval.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_legendre_P1(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_legendre_P1_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_legendre_P2(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_legendre_P2_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_legendre_P3(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_legendre_P3_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_legendre_Pl(const int l, const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_legendre_Pl_e(l, x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_legendre_Plm(const int l, const int m, const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_legendre_Plm_e(l, m, x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_legendre_sphPlm(const int l, const int m, const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_legendre_sphPlm_e(l, m, x, &result));
|
|
Packit |
67cb25 |
}
|