include/gsl/gsl_sf_legendre.h

0001 /* specfunc/gsl_sf_legendre.h
0002  *
0003  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2004 Gerard Jungman
0004  * Copyright (C) 2019 Patrick Alken
0005  *
0006  * This program is free software; you can redistribute it and/or modify
0007  * it under the terms of the GNU General Public License as published by
0008  * the Free Software Foundation; either version 3 of the License, or (at
0009  * your option) any later version.
0010  *
0011  * This program is distributed in the hope that it will be useful, but
0012  * WITHOUT ANY WARRANTY; without even the implied warranty of
0013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0014  * General Public License for more details.
0015  *
0016  * You should have received a copy of the GNU General Public License
0017  * along with this program; if not, write to the Free Software
0018  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
0019  */
0020
0021 /* Author:  G. Jungman */
0022
0023 #ifndef __GSL_SF_LEGENDRE_H__
0024 #define __GSL_SF_LEGENDRE_H__
0025
0026 #include <stdlib.h>
0027 #include <gsl/gsl_inline.h>
0028 #include <gsl/gsl_sf_result.h>
0029
0030 #undef __BEGIN_DECLS
0031 #undef __END_DECLS
0032 #ifdef __cplusplus
0033 # define __BEGIN_DECLS extern "C" {
0034 # define __END_DECLS }
0035 #else
0036 # define __BEGIN_DECLS /* empty */
0037 # define __END_DECLS /* empty */
0038 #endif
0039
0040 __BEGIN_DECLS
0041
0042
0043 /* P_l(x)   l >= 0; |x| <= 1
0044  *
0045  * exceptions: GSL_EDOM
0046  */
0047 int     gsl_sf_legendre_Pl_e(const int l, const double x, gsl_sf_result * result);
0048 double  gsl_sf_legendre_Pl(const int l, const double x);
0049
0050
0051 /* P_l(x) for l=0,...,lmax; |x| <= 1
0052  *
0053  * exceptions: GSL_EDOM
0054  */
0055 int gsl_sf_legendre_Pl_array(
0056   const int lmax, const double x,
0057   double * result_array
0058   );
0059
0060
0061 /* P_l(x) and P_l'(x) for l=0,...,lmax; |x| <= 1
0062  *
0063  * exceptions: GSL_EDOM
0064  */
0065 int gsl_sf_legendre_Pl_deriv_array(
0066   const int lmax, const double x,
0067   double * result_array,
0068   double * result_deriv_array
0069   );
0070
0071
0072 /* P_l(x), l=1,2,3
0073  *
0074  * exceptions: none
0075  */
0076 int gsl_sf_legendre_P1_e(double x, gsl_sf_result * result);
0077 int gsl_sf_legendre_P2_e(double x, gsl_sf_result * result);
0078 int gsl_sf_legendre_P3_e(double x, gsl_sf_result * result);
0079 double gsl_sf_legendre_P1(const double x);
0080 double gsl_sf_legendre_P2(const double x);
0081 double gsl_sf_legendre_P3(const double x);
0082
0083
0084 /* Q_0(x), x > -1, x != 1
0085  *
0086  * exceptions: GSL_EDOM
0087  */
0088 int gsl_sf_legendre_Q0_e(const double x, gsl_sf_result * result);
0089 double gsl_sf_legendre_Q0(const double x);
0090
0091
0092 /* Q_1(x), x > -1, x != 1
0093  *
0094  * exceptions: GSL_EDOM
0095  */
0096 int gsl_sf_legendre_Q1_e(const double x, gsl_sf_result * result);
0097 double gsl_sf_legendre_Q1(const double x);
0098
0099
0100 /* Q_l(x), x > -1, x != 1, l >= 0
0101  *
0102  * exceptions: GSL_EDOM
0103  */
0104 int gsl_sf_legendre_Ql_e(const int l, const double x, gsl_sf_result * result);
0105 double gsl_sf_legendre_Ql(const int l, const double x);
0106
0107
0108 /* P_l^m(x)  m >= 0; l >= m; |x| <= 1.0
0109  *
0110  * Note that this function grows combinatorially with l.
0111  * Therefore we can easily generate an overflow for l larger
0112  * than about 150.
0113  *
0114  * There is no trouble for small m, but when m and l are both large,
0115  * then there will be trouble. Rather than allow overflows, these
0116  * functions refuse to calculate when they can sense that l and m are
0117  * too big.
0118  *
0119  * If you really want to calculate a spherical harmonic, then DO NOT
0120  * use this. Instead use legendre_sphPlm() below, which  uses a similar
0121  * recursion, but with the normalized functions.
0122  *
0123  * exceptions: GSL_EDOM, GSL_EOVRFLW
0124  */
0125 int     gsl_sf_legendre_Plm_e(const int l, const int m, const double x, gsl_sf_result * result);
0126 double  gsl_sf_legendre_Plm(const int l, const int m, const double x);
0127
0128
0129 /* P_l^m(x)  m >= 0; l >= m; |x| <= 1.0
0130  * l=|m|,...,lmax
0131  *
0132  * exceptions: GSL_EDOM, GSL_EOVRFLW
0133  */
0134 int gsl_sf_legendre_Plm_array(
0135   const int lmax, const int m, const double x,
0136   double * result_array
0137   );
0138
0139
0140 /* P_l^m(x)  and d(P_l^m(x))/dx;  m >= 0; lmax >= m; |x| <= 1.0
0141  * l=|m|,...,lmax
0142  *
0143  * exceptions: GSL_EDOM, GSL_EOVRFLW
0144  */
0145 int gsl_sf_legendre_Plm_deriv_array(
0146   const int lmax, const int m, const double x,
0147   double * result_array,
0148   double * result_deriv_array
0149   );
0150
0151
0152 /* P_l^m(x), normalized properly for use in spherical harmonics
0153  * m >= 0; l >= m; |x| <= 1.0
0154  *
0155  * There is no overflow problem, as there is for the
0156  * standard normalization of P_l^m(x).
0157  *
0158  * Specifically, it returns:
0159  *
0160  *        sqrt((2l+1)/(4pi)) sqrt((l-m)!/(l+m)!) P_l^m(x)
0161  *
0162  * exceptions: GSL_EDOM
0163  */
0164 int     gsl_sf_legendre_sphPlm_e(const int l, int m, const double x, gsl_sf_result * result);
0165 double  gsl_sf_legendre_sphPlm(const int l, const int m, const double x);
0166
0167
0168 /* sphPlm(l,m,x) values
0169  * m >= 0; l >= m; |x| <= 1.0
0170  * l=|m|,...,lmax
0171  *
0172  * exceptions: GSL_EDOM
0173  */
0174 int gsl_sf_legendre_sphPlm_array(
0175   const int lmax, int m, const double x,
0176   double * result_array
0177   );
0178
0179
0180 /* sphPlm(l,m,x) and d(sphPlm(l,m,x))/dx values
0181  * m >= 0; l >= m; |x| <= 1.0
0182  * l=|m|,...,lmax
0183  *
0184  * exceptions: GSL_EDOM
0185  */
0186 int gsl_sf_legendre_sphPlm_deriv_array(
0187   const int lmax, const int m, const double x,
0188   double * result_array,
0189   double * result_deriv_array
0190   );
0191
0192
0193
0194 /* size of result_array[] needed for the array versions of Plm
0195  * (lmax - m + 1)
0196  */
0197 int gsl_sf_legendre_array_size(const int lmax, const int m);
0198
0199 /* Irregular Spherical Conical Function
0200  * P^{1/2}_{-1/2 + I lambda}(x)
0201  *
0202  * x > -1.0
0203  * exceptions: GSL_EDOM
0204  */
0205 int gsl_sf_conicalP_half_e(const double lambda, const double x, gsl_sf_result * result);
0206 double gsl_sf_conicalP_half(const double lambda, const double x);
0207
0208
0209 /* Regular Spherical Conical Function
0210  * P^{-1/2}_{-1/2 + I lambda}(x)
0211  *
0212  * x > -1.0
0213  * exceptions: GSL_EDOM
0214  */
0215 int gsl_sf_conicalP_mhalf_e(const double lambda, const double x, gsl_sf_result * result);
0216 double gsl_sf_conicalP_mhalf(const double lambda, const double x);
0217
0218
0219 /* Conical Function
0220  * P^{0}_{-1/2 + I lambda}(x)
0221  *
0222  * x > -1.0
0223  * exceptions: GSL_EDOM
0224  */
0225 int gsl_sf_conicalP_0_e(const double lambda, const double x, gsl_sf_result * result);
0226 double gsl_sf_conicalP_0(const double lambda, const double x);
0227
0228
0229 /* Conical Function
0230  * P^{1}_{-1/2 + I lambda}(x)
0231  *
0232  * x > -1.0
0233  * exceptions: GSL_EDOM
0234  */
0235 int gsl_sf_conicalP_1_e(const double lambda, const double x, gsl_sf_result * result);
0236 double gsl_sf_conicalP_1(const double lambda, const double x);
0237
0238
0239 /* Regular Spherical Conical Function
0240  * P^{-1/2-l}_{-1/2 + I lambda}(x)
0241  *
0242  * x > -1.0, l >= -1
0243  * exceptions: GSL_EDOM
0244  */
0245 int gsl_sf_conicalP_sph_reg_e(const int l, const double lambda, const double x, gsl_sf_result * result);
0246 double gsl_sf_conicalP_sph_reg(const int l, const double lambda, const double x);
0247
0248
0249 /* Regular Cylindrical Conical Function
0250  * P^{-m}_{-1/2 + I lambda}(x)
0251  *
0252  * x > -1.0, m >= -1
0253  * exceptions: GSL_EDOM
0254  */
0255 int gsl_sf_conicalP_cyl_reg_e(const int m, const double lambda, const double x, gsl_sf_result * result);
0256 double gsl_sf_conicalP_cyl_reg(const int m, const double lambda, const double x);
0257
0258
0259 /* The following spherical functions are specializations
0260  * of Legendre functions which give the regular eigenfunctions
0261  * of the Laplacian on a 3-dimensional hyperbolic space.
0262  * Of particular interest is the flat limit, which is
0263  * Flat-Lim := {lambda->Inf, eta->0, lambda*eta fixed}.
0264  */
0265
0266 /* Zeroth radial eigenfunction of the Laplacian on the
0267  * 3-dimensional hyperbolic space.
0268  *
0269  * legendre_H3d_0(lambda,eta) := sin(lambda*eta)/(lambda*sinh(eta))
0270  *
0271  * Normalization:
0272  * Flat-Lim legendre_H3d_0(lambda,eta) = j_0(lambda*eta)
0273  *
0274  * eta >= 0.0
0275  * exceptions: GSL_EDOM
0276  */
0277 int gsl_sf_legendre_H3d_0_e(const double lambda, const double eta, gsl_sf_result * result);
0278 double gsl_sf_legendre_H3d_0(const double lambda, const double eta);
0279
0280
0281 /* First radial eigenfunction of the Laplacian on the
0282  * 3-dimensional hyperbolic space.
0283  *
0284  * legendre_H3d_1(lambda,eta) :=
0285  *    1/sqrt(lambda^2 + 1) sin(lam eta)/(lam sinh(eta))
0286  *    (coth(eta) - lambda cot(lambda*eta))
0287  *
0288  * Normalization:
0289  * Flat-Lim legendre_H3d_1(lambda,eta) = j_1(lambda*eta)
0290  *
0291  * eta >= 0.0
0292  * exceptions: GSL_EDOM
0293  */
0294 int gsl_sf_legendre_H3d_1_e(const double lambda, const double eta, gsl_sf_result * result);
0295 double gsl_sf_legendre_H3d_1(const double lambda, const double eta);
0296
0297
0298 /* l'th radial eigenfunction of the Laplacian on the
0299  * 3-dimensional hyperbolic space.
0300  *
0301  * Normalization:
0302  * Flat-Lim legendre_H3d_l(l,lambda,eta) = j_l(lambda*eta)
0303  *
0304  * eta >= 0.0, l >= 0
0305  * exceptions: GSL_EDOM
0306  */
0307 int gsl_sf_legendre_H3d_e(const int l, const double lambda, const double eta, gsl_sf_result * result);
0308 double gsl_sf_legendre_H3d(const int l, const double lambda, const double eta);
0309
0310
0311 /* Array of H3d(ell),  0 <= ell <= lmax
0312  */
0313 int gsl_sf_legendre_H3d_array(const int lmax, const double lambda, const double eta, double * result_array);
0314
0315 /* associated legendre P_{lm} routines */
0316
0317 typedef enum
0318 {
0319   GSL_SF_LEGENDRE_SCHMIDT,
0320   GSL_SF_LEGENDRE_SPHARM,
0321   GSL_SF_LEGENDRE_FULL,
0322   GSL_SF_LEGENDRE_NONE
0323 } gsl_sf_legendre_t;
0324
0325 int gsl_sf_legendre_array(const gsl_sf_legendre_t norm,
0326                           const size_t lmax, const double x,
0327                           double result_array[]);
0328 int gsl_sf_legendre_array_e(const gsl_sf_legendre_t norm,
0329                             const size_t lmax, const double x,
0330                             const double csphase,
0331                             double result_array[]);
0332 int gsl_sf_legendre_deriv_array(const gsl_sf_legendre_t norm,
0333                                 const size_t lmax, const double x,
0334                                 double result_array[],
0335                                 double result_deriv_array[]);
0336 int gsl_sf_legendre_deriv_array_e(const gsl_sf_legendre_t norm,
0337                                   const size_t lmax, const double x,
0338                                   const double csphase,
0339                                   double result_array[],
0340                                   double result_deriv_array[]);
0341 int gsl_sf_legendre_deriv_alt_array(const gsl_sf_legendre_t norm,
0342                                     const size_t lmax, const double x,
0343                                     double result_array[],
0344                                     double result_deriv_array[]);
0345 int gsl_sf_legendre_deriv_alt_array_e(const gsl_sf_legendre_t norm,
0346                                       const size_t lmax, const double x,
0347                                       const double csphase,
0348                                       double result_array[],
0349                                       double result_deriv_array[]);
0350 int gsl_sf_legendre_deriv2_array(const gsl_sf_legendre_t norm,
0351                                  const size_t lmax, const double x,
0352                                  double result_array[],
0353                                  double result_deriv_array[],
0354                                  double result_deriv2_array[]);
0355 int gsl_sf_legendre_deriv2_array_e(const gsl_sf_legendre_t norm,
0356                                    const size_t lmax, const double x,
0357                                    const double csphase,
0358                                    double result_array[],
0359                                    double result_deriv_array[],
0360                                    double result_deriv2_array[]);
0361 int gsl_sf_legendre_deriv2_alt_array(const gsl_sf_legendre_t norm,
0362                                      const size_t lmax, const double x,
0363                                      double result_array[],
0364                                      double result_deriv_array[],
0365                                      double result_deriv2_array[]);
0366 int gsl_sf_legendre_deriv2_alt_array_e(const gsl_sf_legendre_t norm,
0367                                        const size_t lmax, const double x,
0368                                        const double csphase,
0369                                        double result_array[],
0370                                        double result_deriv_array[],
0371                                        double result_deriv2_array[]);
0372 size_t gsl_sf_legendre_array_n(const size_t lmax);
0373 size_t gsl_sf_legendre_nlm(const size_t lmax);
0374
0375 INLINE_DECL size_t gsl_sf_legendre_array_index(const size_t l, const size_t m);
0376
0377 #ifdef HAVE_INLINE
0378
0379 /*
0380 gsl_sf_legendre_array_index()
0381 This routine computes the index into a result_array[] corresponding
0382 to a given (l,m)
0383 */
0384 INLINE_FUN
0385 size_t
0386 gsl_sf_legendre_array_index(const size_t l, const size_t m)
0387 {
0388   return (((l * (l + 1)) >> 1) + m);
0389 }
0390
0391 #endif /* HAVE_INLINE */
0392
0393 __END_DECLS
0394
0395 #endif /* __GSL_SF_LEGENDRE_H__ */