htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/legendre_poly.c

   1 /* specfunc/legendre_poly.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman */
  21
  22 #include <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_sf_bessel.h>
  26 #include <gsl/gsl_sf_exp.h>
  27 #include <gsl/gsl_sf_gamma.h>
  28 #include <gsl/gsl_sf_log.h>
  29 #include <gsl/gsl_sf_pow_int.h>
  30 #include <gsl/gsl_sf_legendre.h>
  31
  32 #include "error.h"
  33
  34
  35
  36 /* Calculate P_m^m(x) from the analytic result:
  37  *   P_m^m(x) = (-1)^m (2m-1)!! (1-x^2)^(m/2) , m > 0
  38  *            = 1 , m = 0
  39  */
  40 static double legendre_Pmm(int m, double x)
  41 {
  42   if(m == 0)
  43   {
  44     return 1.0;
  45   }
  46   else
  47   {
  48     double p_mm = 1.0;
  49     double root_factor = sqrt(1.0-x)*sqrt(1.0+x);
  50     double fact_coeff = 1.0;
  51     int i;
  52     for(i=1; i<=m; i++)
  53     {
  54       p_mm *= -fact_coeff * root_factor;
  55       fact_coeff += 2.0;
  56     }
  57     return p_mm;
  58   }
  59 }
  60
  61
  62
  63 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  64
  65 int
  66 gsl_sf_legendre_P1_e(double x, gsl_sf_result * result)
  67 {
  68   /* CHECK_POINTER(result) */
  69
  70   {
  71     result->val = x;
  72     result->err = 0.0;
  73     return GSL_SUCCESS;
  74   }
  75 }
  76
  77 int
  78 gsl_sf_legendre_P2_e(double x, gsl_sf_result * result)
  79 {
  80   /* CHECK_POINTER(result) */
  81
  82   {
  83     result->val = 0.5*(3.0*x*x - 1.0);
  84     result->err = GSL_DBL_EPSILON * (fabs(3.0*x*x) + 1.0);
  85     return GSL_SUCCESS;
  86   }
  87 }
  88
  89 int
  90 gsl_sf_legendre_P3_e(double x, gsl_sf_result * result)
  91 {
  92   /* CHECK_POINTER(result) */
  93
  94   {
  95     result->val = 0.5*x*(5.0*x*x - 3.0);
  96     result->err = GSL_DBL_EPSILON * (fabs(result->val) + 0.5 * fabs(x) * (fabs(5.0*x*x) + 3.0));
  97     return GSL_SUCCESS;
  98   }
  99 }
 100
 101
 102 int
 103 gsl_sf_legendre_Pl_e(const int l, const double x, gsl_sf_result * result)
 104 {
 105   /* CHECK_POINTER(result) */
 106
 107   if(l < 0 || x < -1.0 || x > 1.0) {
 108     DOMAIN_ERROR(result);
 109   }
 110   else if(l == 0) {
 111     result->val = 1.0;
 112     result->err = 0.0;
 113     return GSL_SUCCESS;
 114   }
 115   else if(l == 1) {
 116     result->val = x;
 117     result->err = 0.0;
 118     return GSL_SUCCESS;
 119   }
 120   else if(l == 2) {
 121     result->val = 0.5 * (3.0*x*x - 1.0);
 122     result->err = GSL_DBL_EPSILON * (fabs(3.0*x*x) + 1.0);
 123     /*result->err = 3.0 * GSL_DBL_EPSILON * fabs(result->val);
 124       removed this old bogus estimate [GJ]
 125       */
 126     return GSL_SUCCESS;
 127   }
 128   else if(x == 1.0) {
 129     result->val = 1.0;
 130     result->err = 0.0;
 131     return GSL_SUCCESS;
 132   }
 133   else if(x == -1.0) {
 134     result->val = ( GSL_IS_ODD(l) ? -1.0 : 1.0 );
 135     result->err = 0.0;
 136     return GSL_SUCCESS;
 137   }
 138   else if(l < 100000) {
 139     /* upward recurrence: l P_l = (2l-1) z P_{l-1} - (l-1) P_{l-2} */
 140
 141     double p_ellm2 = 1.0;    /* P_0(x) */
 142     double p_ellm1 = x;      /* P_1(x) */
 143     double p_ell = p_ellm1;
 144
 145     double e_ellm2 = GSL_DBL_EPSILON;
 146     double e_ellm1 = fabs(x)*GSL_DBL_EPSILON;
 147     double e_ell = e_ellm1;
 148
 149     int ell;
 150
 151     for(ell=2; ell <= l; ell++){
 152       p_ell = (x*(2*ell-1)*p_ellm1 - (ell-1)*p_ellm2) / ell;
 153       p_ellm2 = p_ellm1;
 154       p_ellm1 = p_ell;
 155
 156       e_ell = 0.5*(fabs(x)*(2*ell-1.0) * e_ellm1 + (ell-1.0)*e_ellm2)/ell;
 157       e_ellm2 = e_ellm1;
 158       e_ellm1 = e_ell;
 159     }
 160
 161     result->val = p_ell;
 162     result->err = e_ell + l*fabs(p_ell)*GSL_DBL_EPSILON;
 163     return GSL_SUCCESS;
 164   }
 165   else {
 166     /* Asymptotic expansion.
 167      * [Olver, p. 473]
 168      */
 169     double u  = l + 0.5;
 170     double th = acos(x);
 171     gsl_sf_result J0;
 172     gsl_sf_result Jm1;
 173     int stat_J0  = gsl_sf_bessel_J0_e(u*th, &J0);
 174     int stat_Jm1 = gsl_sf_bessel_Jn_e(-1, u*th, &Jm1);
 175     double pre;
 176     double B00;
 177     double c1;
 178
 179     /* B00 = 1/8 (1 - th cot(th) / th^2
 180      * pre = sqrt(th/sin(th))
 181      */
 182     if(th < GSL_ROOT4_DBL_EPSILON) {
 183       B00 = (1.0 + th*th/15.0)/24.0;
 184       pre = 1.0 + th*th/12.0;
 185     }
 186     else {
 187       double sin_th = sqrt(1.0 - x*x);
 188       double cot_th = x / sin_th;
 189       B00 = 1.0/8.0 * (1.0 - th * cot_th) / (th*th);
 190       pre = sqrt(th/sin_th);
 191     }
 192
 193     c1 = th/u * B00;
 194
 195     result->val  = pre * (J0.val + c1 * Jm1.val);
 196     result->err  = pre * (J0.err + fabs(c1) * Jm1.err);
 197     result->err += GSL_SQRT_DBL_EPSILON * fabs(result->val);
 198
 199     return GSL_ERROR_SELECT_2(stat_J0, stat_Jm1);
 200   }
 201 }
 202
 203
 204 int
 205 gsl_sf_legendre_Pl_array(const int lmax, const double x, double * result_array)
 206 {
 207   /* CHECK_POINTER(result_array) */
 208
 209   if(lmax < 0 || x < -1.0 || x > 1.0) {
 210     GSL_ERROR ("domain error", GSL_EDOM);
 211   }
 212   else if(lmax == 0) {
 213     result_array[0] = 1.0;
 214     return GSL_SUCCESS;
 215   }
 216   else if(lmax == 1) {
 217     result_array[0] = 1.0;
 218     result_array[1] = x;
 219     return GSL_SUCCESS;
 220   }
 221   else {
 222     /* upward recurrence: l P_l = (2l-1) z P_{l-1} - (l-1) P_{l-2} */
 223
 224     double p_ellm2 = 1.0;    /* P_0(x) */
 225     double p_ellm1 = x;    /* P_1(x) */
 226     double p_ell = p_ellm1;
 227     int ell;
 228
 229     result_array[0] = 1.0;
 230     result_array[1] = x;
 231
 232     for(ell=2; ell <= lmax; ell++){
 233       p_ell = (x*(2*ell-1)*p_ellm1 - (ell-1)*p_ellm2) / ell;
 234       p_ellm2 = p_ellm1;
 235       p_ellm1 = p_ell;
 236       result_array[ell] = p_ell;
 237     }
 238
 239     return GSL_SUCCESS;
 240   }
 241 }
 242
 243
 244 int
 245 gsl_sf_legendre_Pl_deriv_array(const int lmax, const double x, double * result_array, double * result_deriv_array)
 246 {
 247   int stat_array = gsl_sf_legendre_Pl_array(lmax, x, result_array);
 248
 249   if(lmax >= 0) result_deriv_array[0] = 0.0;
 250   if(lmax >= 1) result_deriv_array[1] = 1.0;
 251
 252   if(stat_array == GSL_SUCCESS)
 253   {
 254     int ell;
 255
 256     if(fabs(x - 1.0)*(lmax+1.0)*(lmax+1.0) <  GSL_SQRT_DBL_EPSILON)
 257     {
 258       /* x is near 1 */
 259       for(ell = 2; ell <= lmax; ell++)
 260       {
 261         const double pre = 0.5 * ell * (ell+1.0);
 262         result_deriv_array[ell] = pre * (1.0 - 0.25 * (1.0-x) * (ell+2.0)*(ell-1.0));
 263       }
 264     }
 265     else if(fabs(x + 1.0)*(lmax+1.0)*(lmax+1.0) <  GSL_SQRT_DBL_EPSILON)
 266     {
 267       /* x is near -1 */
 268       for(ell = 2; ell <= lmax; ell++)
 269       {
 270         const double sgn = ( GSL_IS_ODD(ell) ? 1.0 : -1.0 ); /* derivative is odd in x for even ell */
 271         const double pre = sgn * 0.5 * ell * (ell+1.0);
 272         result_deriv_array[ell] = pre * (1.0 - 0.25 * (1.0+x) * (ell+2.0)*(ell-1.0));
 273       }
 274     }
 275     else
 276     {
 277       const double diff_a = 1.0 + x;
 278       const double diff_b = 1.0 - x;
 279       for(ell = 2; ell <= lmax; ell++)
 280       {
 281         result_deriv_array[ell] = - ell * (x * result_array[ell] - result_array[ell-1]) / (diff_a * diff_b);
 282       }
 283     }
 284
 285     return GSL_SUCCESS;
 286   }
 287   else
 288   {
 289     return stat_array;
 290   }
 291 }
 292
 293
 294 int
 295 gsl_sf_legendre_Plm_e(const int l, const int m, const double x, gsl_sf_result * result)
 296 {
 297   /* If l is large and m is large, then we have to worry
 298    * about overflow. Calculate an approximate exponent which
 299    * measures the normalization of this thing.
 300    */
 301   const double dif = l-m;
 302   const double sum = l+m;
 303   const double t_d = ( dif == 0.0 ? 0.0 : 0.5 * dif * (log(dif)-1.0) );
 304   const double t_s = ( dif == 0.0 ? 0.0 : 0.5 * sum * (log(sum)-1.0) );
 305   const double exp_check = 0.5 * log(2.0*l+1.0) + t_d - t_s;
 306
 307   /* CHECK_POINTER(result) */
 308
 309   if(m < 0 || l < m || x < -1.0 || x > 1.0) {
 310     DOMAIN_ERROR(result);
 311   }
 312   else if(exp_check < GSL_LOG_DBL_MIN + 10.0){
 313     /* Bail out. */
 314     OVERFLOW_ERROR(result);
 315   }
 316   else {
 317     /* Account for the error due to the
 318      * representation of 1-x.
 319      */
 320     const double err_amp = 1.0 / (GSL_DBL_EPSILON + fabs(1.0-fabs(x)));
 321
 322     /* P_m^m(x) and P_{m+1}^m(x) */
 323     double p_mm   = legendre_Pmm(m, x);
 324     double p_mmp1 = x * (2*m + 1) * p_mm;
 325
 326     if(l == m){
 327       result->val = p_mm;
 328       result->err = err_amp * 2.0 * GSL_DBL_EPSILON * fabs(p_mm);
 329       return GSL_SUCCESS;
 330     }
 331     else if(l == m + 1) {
 332       result->val = p_mmp1;
 333       result->err = err_amp * 2.0 * GSL_DBL_EPSILON * fabs(p_mmp1);
 334       return GSL_SUCCESS;
 335     }
 336     else{
 337       /* upward recurrence: (l-m) P(l,m) = (2l-1) z P(l-1,m) - (l+m-1) P(l-2,m)
 338        * start at P(m,m), P(m+1,m)
 339        */
 340
 341       double p_ellm2 = p_mm;
 342       double p_ellm1 = p_mmp1;
 343       double p_ell = 0.0;
 344       int ell;
 345
 346       for(ell=m+2; ell <= l; ell++){
 347         p_ell = (x*(2*ell-1)*p_ellm1 - (ell+m-1)*p_ellm2) / (ell-m);
 348         p_ellm2 = p_ellm1;
 349         p_ellm1 = p_ell;
 350       }
 351
 352       result->val = p_ell;
 353       result->err = err_amp * (0.5*(l-m) + 1.0) * GSL_DBL_EPSILON * fabs(p_ell);
 354
 355       return GSL_SUCCESS;
 356     }
 357   }
 358 }
 359
 360
 361 int
 362 gsl_sf_legendre_Plm_array(const int lmax, const int m, const double x, double * result_array)
 363 {
 364   /* If l is large and m is large, then we have to worry
 365    * about overflow. Calculate an approximate exponent which
 366    * measures the normalization of this thing.
 367    */
 368   const double dif = lmax-m;
 369   const double sum = lmax+m;
 370   const double t_d = ( dif == 0.0 ? 0.0 : 0.5 * dif * (log(dif)-1.0) );
 371   const double t_s = ( dif == 0.0 ? 0.0 : 0.5 * sum * (log(sum)-1.0) );
 372   const double exp_check = 0.5 * log(2.0*lmax+1.0) + t_d - t_s;
 373
 374   /* CHECK_POINTER(result_array) */
 375
 376   if(m < 0 || lmax < m || x < -1.0 || x > 1.0) {
 377     GSL_ERROR ("domain error", GSL_EDOM);
 378   }
 379   else if(m > 0 && (x == 1.0 || x == -1.0)) {
 380     int ell;
 381     for(ell=m; ell<=lmax; ell++) result_array[ell-m] = 0.0;
 382     return GSL_SUCCESS;
 383   }
 384   else if(exp_check < GSL_LOG_DBL_MIN + 10.0){
 385     /* Bail out. */
 386     GSL_ERROR ("overflow", GSL_EOVRFLW);
 387   }
 388   else {
 389     double p_mm   = legendre_Pmm(m, x);
 390     double p_mmp1 = x * (2.0*m + 1.0) * p_mm;
 391
 392     if(lmax == m){
 393       result_array[0] = p_mm;
 394       return GSL_SUCCESS;
 395     }
 396     else if(lmax == m + 1) {
 397       result_array[0] = p_mm;
 398       result_array[1] = p_mmp1;
 399       return GSL_SUCCESS;
 400     }
 401     else {
 402       double p_ellm2 = p_mm;
 403       double p_ellm1 = p_mmp1;
 404       double p_ell = 0.0;
 405       int ell;
 406
 407       result_array[0] = p_mm;
 408       result_array[1] = p_mmp1;
 409
 410       for(ell=m+2; ell <= lmax; ell++){
 411         p_ell = (x*(2.0*ell-1.0)*p_ellm1 - (ell+m-1)*p_ellm2) / (ell-m);
 412         p_ellm2 = p_ellm1;
 413         p_ellm1 = p_ell;
 414         result_array[ell-m] = p_ell;
 415       }
 416
 417       return GSL_SUCCESS;
 418     }
 419   }
 420 }
 421
 422
 423 int
 424 gsl_sf_legendre_Plm_deriv_array(
 425   const int lmax, const int m, const double x,
 426   double * result_array,
 427   double * result_deriv_array)
 428 {
 429   if(m < 0 || m > lmax)
 430   {
 431     GSL_ERROR("m < 0 or m > lmax", GSL_EDOM);
 432   }
 433   else if(m == 0)
 434   {
 435     /* It is better to do m=0 this way, so we can more easily
 436      * trap the divergent case which can occur when m == 1.
 437      */
 438     return gsl_sf_legendre_Pl_deriv_array(lmax, x, result_array, result_deriv_array);
 439   }
 440   else
 441   {
 442     int stat_array = gsl_sf_legendre_Plm_array(lmax, m, x, result_array);
 443
 444     if(stat_array == GSL_SUCCESS)
 445     {
 446       int ell;
 447
 448       if(m == 1 && (1.0 - fabs(x) < GSL_DBL_EPSILON))
 449       {
 450         /* This divergence is real and comes from the cusp-like
 451          * behaviour for m = 1. For example, P[1,1] = - Sqrt[1-x^2].
 452          */
 453         GSL_ERROR("divergence near |x| = 1.0 since m = 1", GSL_EOVRFLW);
 454       }
 455       else if(m == 2 && (1.0 - fabs(x) < GSL_DBL_EPSILON))
 456       {
 457         /* m = 2 gives a finite nonzero result for |x| near 1 */
 458         if(fabs(x - 1.0) < GSL_DBL_EPSILON)
 459         {
 460           for(ell = m; ell <= lmax; ell++) result_deriv_array[ell-m] = -0.25 * x * (ell - 1.0)*ell*(ell+1.0)*(ell+2.0);
 461         }
 462         else if(fabs(x + 1.0) < GSL_DBL_EPSILON)
 463         {
 464           for(ell = m; ell <= lmax; ell++)
 465           {
 466             const double sgn = ( GSL_IS_ODD(ell) ? 1.0 : -1.0 );
 467             result_deriv_array[ell-m] = -0.25 * sgn * x * (ell - 1.0)*ell*(ell+1.0)*(ell+2.0);
 468           }
 469         }
 470         return GSL_SUCCESS;
 471       }
 472       else
 473       {
 474         /* m > 2 is easier to deal with since the endpoints always vanish */
 475         if(1.0 - fabs(x) < GSL_DBL_EPSILON)
 476         {
 477           for(ell = m; ell <= lmax; ell++) result_deriv_array[ell-m] = 0.0;
 478           return GSL_SUCCESS;
 479         }
 480         else
 481         {
 482           const double diff_a = 1.0 + x;
 483           const double diff_b = 1.0 - x;
 484           result_deriv_array[0] = - m * x / (diff_a * diff_b) * result_array[0];
 485           if(lmax-m >= 1) result_deriv_array[1] = (2.0 * m + 1.0) * (x * result_deriv_array[0] + result_array[0]);
 486           for(ell = m+2; ell <= lmax; ell++)
 487           {
 488             result_deriv_array[ell-m] = - (ell * x * result_array[ell-m] - (ell+m) * result_array[ell-1-m]) / (diff_a * diff_b);
 489           }
 490           return GSL_SUCCESS;
 491         }
 492       }
 493     }
 494     else
 495     {
 496       return stat_array;
 497     }
 498   }
 499 }
 500
 501
 502 int
 503 gsl_sf_legendre_sphPlm_e(const int l, int m, const double x, gsl_sf_result * result)
 504 {
 505   /* CHECK_POINTER(result) */
 506
 507   if(m < 0 || l < m || x < -1.0 || x > 1.0) {
 508     DOMAIN_ERROR(result);
 509   }
 510   else if(m == 0) {
 511     gsl_sf_result P;
 512     int stat_P = gsl_sf_legendre_Pl_e(l, x, &P);
 513     double pre = sqrt((2.0*l + 1.0)/(4.0*M_PI));
 514     result->val  = pre * P.val;
 515     result->err  = pre * P.err;
 516     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 517     return stat_P;
 518   }
 519   else if(x == 1.0 || x == -1.0) {
 520     /* m > 0 here */
 521     result->val = 0.0;
 522     result->err = 0.0;
 523     return GSL_SUCCESS;
 524   }
 525   else {
 526     /* m > 0 and |x| < 1 here */
 527
 528     /* Starting value for recursion.
 529      * 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)
 530      */
 531     gsl_sf_result lncirc;
 532     gsl_sf_result lnpoch;
 533     double lnpre_val;
 534     double lnpre_err;
 535     gsl_sf_result ex_pre;
 536     double sr;
 537     const double sgn = ( GSL_IS_ODD(m) ? -1.0 : 1.0);
 538     const double y_mmp1_factor = x * sqrt(2.0*m + 3.0);
 539     double y_mm, y_mm_err;
 540     double y_mmp1, y_mmp1_err;
 541     gsl_sf_log_1plusx_e(-x*x, &lncirc);
 542     gsl_sf_lnpoch_e(m, 0.5, &lnpoch);  /* Gamma(m+1/2)/Gamma(m) */
 543     lnpre_val = -0.25*M_LNPI + 0.5 * (lnpoch.val + m*lncirc.val);
 544     lnpre_err = 0.25*M_LNPI*GSL_DBL_EPSILON + 0.5 * (lnpoch.err + fabs(m)*lncirc.err);
 545     /* Compute exp(ln_pre) with error term, avoiding call to gsl_sf_exp_err BJG */
 546     ex_pre.val = exp(lnpre_val);
 547     ex_pre.err = 2.0*(sinh(lnpre_err) + GSL_DBL_EPSILON)*ex_pre.val;
 548     sr     = sqrt((2.0+1.0/m)/(4.0*M_PI));
 549     y_mm   = sgn * sr * ex_pre.val;
 550     y_mm_err  = 2.0 * GSL_DBL_EPSILON * fabs(y_mm) + sr * ex_pre.err;
 551     y_mm_err *= 1.0 + 1.0/(GSL_DBL_EPSILON + fabs(1.0-x));
 552     y_mmp1 = y_mmp1_factor * y_mm;
 553     y_mmp1_err=fabs(y_mmp1_factor) * y_mm_err;
 554
 555     if(l == m){
 556       result->val  = y_mm;
 557       result->err  = y_mm_err;
 558       result->err += 2.0 * GSL_DBL_EPSILON * fabs(y_mm);
 559       return GSL_SUCCESS;
 560     }
 561     else if(l == m + 1) {
 562       result->val  = y_mmp1;
 563       result->err  = y_mmp1_err;
 564       result->err += 2.0 * GSL_DBL_EPSILON * fabs(y_mmp1);
 565       return GSL_SUCCESS;
 566     }
 567     else{
 568       double y_ell = 0.0;
 569       double y_ell_err;
 570       int ell;
 571
 572       /* Compute Y_l^m, l > m+1, upward recursion on l. */
 573       for(ell=m+2; ell <= l; ell++){
 574         const double rat1 = (double)(ell-m)/(double)(ell+m);
 575         const double rat2 = (ell-m-1.0)/(ell+m-1.0);
 576         const double factor1 = sqrt(rat1*(2.0*ell+1.0)*(2.0*ell-1.0));
 577         const double factor2 = sqrt(rat1*rat2*(2.0*ell+1.0)/(2.0*ell-3.0));
 578         y_ell = (x*y_mmp1*factor1 - (ell+m-1.0)*y_mm*factor2) / (ell-m);
 579         y_mm   = y_mmp1;
 580         y_mmp1 = y_ell;
 581
 582         y_ell_err = 0.5*(fabs(x*factor1)*y_mmp1_err + fabs((ell+m-1.0)*factor2)*y_mm_err) / fabs(ell-m);
 583         y_mm_err = y_mmp1_err;
 584         y_mmp1_err = y_ell_err;
 585       }
 586
 587       result->val  = y_ell;
 588       result->err  = y_ell_err + (0.5*(l-m) + 1.0) * GSL_DBL_EPSILON * fabs(y_ell);
 589
 590       return GSL_SUCCESS;
 591     }
 592   }
 593 }
 594
 595
 596 int
 597 gsl_sf_legendre_sphPlm_array(const int lmax, int m, const double x, double * result_array)
 598 {
 599   /* CHECK_POINTER(result_array) */
 600
 601   if(m < 0 || lmax < m || x < -1.0 || x > 1.0) {
 602     GSL_ERROR ("error", GSL_EDOM);
 603   }
 604   else if(m > 0 && (x == 1.0 || x == -1.0)) {
 605     int ell;
 606     for(ell=m; ell<=lmax; ell++) result_array[ell-m] = 0.0;
 607     return GSL_SUCCESS;
 608   }
 609   else {
 610     double y_mm;
 611     double y_mmp1;
 612
 613     if(m == 0) {
 614       y_mm   = 0.5/M_SQRTPI;          /* Y00 = 1/sqrt(4pi) */
 615       y_mmp1 = x * M_SQRT3 * y_mm;
 616     }
 617     else {
 618       /* |x| < 1 here */
 619
 620       gsl_sf_result lncirc;
 621       gsl_sf_result lnpoch;
 622       double lnpre;
 623       const double sgn = ( GSL_IS_ODD(m) ? -1.0 : 1.0);
 624       gsl_sf_log_1plusx_e(-x*x, &lncirc);
 625       gsl_sf_lnpoch_e(m, 0.5, &lnpoch);  /* Gamma(m+1/2)/Gamma(m) */
 626       lnpre = -0.25*M_LNPI + 0.5 * (lnpoch.val + m*lncirc.val);
 627       y_mm   = sqrt((2.0+1.0/m)/(4.0*M_PI)) * sgn * exp(lnpre);
 628       y_mmp1 = x * sqrt(2.0*m + 3.0) * y_mm;
 629     }
 630
 631     if(lmax == m){
 632       result_array[0] = y_mm;
 633       return GSL_SUCCESS;
 634     }
 635     else if(lmax == m + 1) {
 636       result_array[0] = y_mm;
 637       result_array[1] = y_mmp1;
 638       return GSL_SUCCESS;
 639     }
 640     else{
 641       double y_ell;
 642       int ell;
 643
 644       result_array[0] = y_mm;
 645       result_array[1] = y_mmp1;
 646
 647       /* Compute Y_l^m, l > m+1, upward recursion on l. */
 648       for(ell=m+2; ell <= lmax; ell++){
 649         const double rat1 = (double)(ell-m)/(double)(ell+m);
 650         const double rat2 = (ell-m-1.0)/(ell+m-1.0);
 651         const double factor1 = sqrt(rat1*(2*ell+1)*(2*ell-1));
 652         const double factor2 = sqrt(rat1*rat2*(2*ell+1)/(2*ell-3));
 653         y_ell = (x*y_mmp1*factor1 - (ell+m-1)*y_mm*factor2) / (ell-m);
 654         y_mm   = y_mmp1;
 655         y_mmp1 = y_ell;
 656         result_array[ell-m] = y_ell;
 657       }
 658     }
 659
 660     return GSL_SUCCESS;
 661   }
 662 }
 663
 664
 665 int
 666 gsl_sf_legendre_sphPlm_deriv_array(
 667   const int lmax, const int m, const double x,
 668   double * result_array,
 669   double * result_deriv_array)
 670 {
 671   if(m < 0 || lmax < m || x < -1.0 || x > 1.0)
 672   {
 673     GSL_ERROR ("domain", GSL_EDOM);
 674   }
 675   else if(m == 0)
 676   {
 677     /* m = 0 is easy to trap */
 678     const int stat_array = gsl_sf_legendre_Pl_deriv_array(lmax, x, result_array, result_deriv_array);
 679     int ell;
 680     for(ell = 0; ell <= lmax; ell++)
 681     {
 682       const double prefactor = sqrt((2.0 * ell + 1.0)/(4.0*M_PI));
 683       result_array[ell] *= prefactor;
 684       result_deriv_array[ell] *= prefactor;
 685     }
 686     return stat_array;
 687   }
 688   else if(m == 1)
 689   {
 690     /* Trapping m = 1 is necessary because of the possible divergence.
 691      * Recall that this divergence is handled properly in ..._Plm_deriv_array(),
 692      * and the scaling factor is not large for small m, so we just scale.
 693      */
 694     const int stat_array = gsl_sf_legendre_Plm_deriv_array(lmax, m, x, result_array, result_deriv_array);
 695     int ell;
 696     for(ell = 1; ell <= lmax; ell++)
 697     {
 698       const double prefactor = sqrt((2.0 * ell + 1.0)/(ell + 1.0) / (4.0*M_PI*ell));
 699       result_array[ell-1] *= prefactor;
 700       result_deriv_array[ell-1] *= prefactor;
 701     }
 702     return stat_array;
 703   }
 704   else
 705   {
 706     /* as for the derivative of P_lm, everything is regular for m >= 2 */
 707
 708     int stat_array = gsl_sf_legendre_sphPlm_array(lmax, m, x, result_array);
 709
 710     if(stat_array == GSL_SUCCESS)
 711     {
 712       int ell;
 713
 714       if(1.0 - fabs(x) < GSL_DBL_EPSILON)
 715       {
 716         for(ell = m; ell <= lmax; ell++) result_deriv_array[ell-m] = 0.0;
 717         return GSL_SUCCESS;
 718       }
 719       else
 720       {
 721         const double diff_a = 1.0 + x;
 722         const double diff_b = 1.0 - x;
 723         result_deriv_array[0] = - m * x / (diff_a * diff_b) * result_array[0];
 724         if(lmax-m >= 1) result_deriv_array[1] = sqrt(2.0 * m + 3.0) * (x * result_deriv_array[0] + result_array[0]);
 725         for(ell = m+2; ell <= lmax; ell++)
 726         {
 727           const double c1 = sqrt(((2.0*ell+1.0)/(2.0*ell-1.0)) * ((double)(ell-m)/(double)(ell+m)));
 728           result_deriv_array[ell-m] = - (ell * x * result_array[ell-m] - c1 * (ell+m) * result_array[ell-1-m]) / (diff_a * diff_b);
 729         }
 730         return GSL_SUCCESS;
 731       }
 732     }
 733     else
 734     {
 735       return stat_array;
 736     }
 737   }
 738 }
 739
 740
 741 #ifndef HIDE_INLINE_STATIC
 742 int
 743 gsl_sf_legendre_array_size(const int lmax, const int m)
 744 {
 745   return lmax-m+1;
 746 }
 747 #endif
 748
 749
 750 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 751
 752 #include "eval.h"
 753
 754 double gsl_sf_legendre_P1(const double x)
 755 {
 756   EVAL_RESULT(gsl_sf_legendre_P1_e(x, &result));
 757 }
 758
 759 double gsl_sf_legendre_P2(const double x)
 760 {
 761   EVAL_RESULT(gsl_sf_legendre_P2_e(x, &result));
 762 }
 763
 764 double gsl_sf_legendre_P3(const double x)
 765 {
 766   EVAL_RESULT(gsl_sf_legendre_P3_e(x, &result));
 767 }
 768
 769 double gsl_sf_legendre_Pl(const int l, const double x)
 770 {
 771   EVAL_RESULT(gsl_sf_legendre_Pl_e(l, x, &result));
 772 }
 773
 774 double gsl_sf_legendre_Plm(const int l, const int m, const double x)
 775 {
 776   EVAL_RESULT(gsl_sf_legendre_Plm_e(l, m, x, &result));
 777 }
 778
 779 double gsl_sf_legendre_sphPlm(const int l, const int m, const double x)
 780 {
 781   EVAL_RESULT(gsl_sf_legendre_sphPlm_e(l, m, x, &result));
 782 }
 783