1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
|
/* legendre_P.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 2 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.
*/
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_legendre.h>
/*
* The routines in this module compute associated Legendre functions
* (ALFs) up to order and degree 2700, using the method described
* in
*
* [1] S. A. Holmes and W. E. Featherstone, A unified approach
* to the Clenshaw summation and the recursive computation of very
* high degree and order normalised associated Legendre functions,
* Journal of Geodesy, 76, pg. 279-299, 2002.
*
* Further information on ALFs can be found in
*
* [2] Abramowitz and Stegun, Handbook of Mathematical Functions,
* Chapter 8, 1972.
*/
static void legendre_sqrts(const size_t lmax, double *array);
#define LEGENDRE
#include "legendre_source.c"
#undef LEGENDRE
#define LEGENDRE_DERIV
#include "legendre_source.c"
#undef LEGENDRE_DERIV
#define LEGENDRE_DERIV_ALT
#include "legendre_source.c"
#undef LEGENDRE_DERIV_ALT
#define LEGENDRE_DERIV2
#include "legendre_source.c"
#undef LEGENDRE_DERIV2
#define LEGENDRE_DERIV2_ALT
#include "legendre_source.c"
#undef LEGENDRE_DERIV2_ALT
/* number of P_{lm} functions for a given lmax */
size_t
gsl_sf_legendre_nlm(const size_t lmax)
{
return ((lmax + 1) * (lmax + 2) / 2);
}
/*
gsl_sf_legendre_array_n()
This routine returns the minimum result_array[] size needed
for a given lmax
*/
size_t
gsl_sf_legendre_array_n(const size_t lmax)
{
size_t nlm = gsl_sf_legendre_nlm(lmax);
size_t nsqrt = 2 * lmax + 2; /* extra room to precompute sqrt factors */
return (nlm + nsqrt);
} /* gsl_sf_legendre_array_n() */
/*
gsl_sf_legendre_array_index()
This routine computes the index into a result_array[] corresponding
to a given (l,m)
*/
size_t
gsl_sf_legendre_array_index(const size_t l, const size_t m)
{
return (l * (l + 1) / 2 + m);
} /* gsl_sf_legendre_array_index() */
/*********************************************************
* INTERNAL ROUTINES *
*********************************************************/
/*
legendre_sqrts()
Precompute square root factors needed for Legendre recurrence.
On output, array[i] = sqrt(i)
*/
static void
legendre_sqrts(const size_t lmax, double *array)
{
size_t l;
for (l = 0; l <= 2 * lmax + 1; ++l)
array[l] = sqrt((double) l);
}
|
