Context Navigation

legendre.lisp

Revision 47, 15.1 kB (checked in by lhealy, 8 months ago)
Unification of errors and warnings using a single class 'gsl-condition. Each numbered GSL conditions is a subclass of this condition, under the name given by GSL, e.g. 'EDOM.
Property svn:keywords set to `Id`

Line
1	;; Legendre functions
2	;; Liam Healy, Sat Apr 29 2006 - 19:16
3	;; Time-stamp: <2008-03-27 21:30:02EDT legendre.lisp>
4	;; $Id$
5
6	(in-package :gsl)
7
8	;;; legendre-Plm-deriv-array same answer as legendre-Plm-array?
9
10	;;;;****************************************************************************
11	;;;; Legendre polynomials
12	;;;;****************************************************************************
13
14	(defmfun legendre-P1 (x)
15	"gsl_sf_legendre_P1_e" ((x :double) (ret sf-result))
16	:documentation ; FDL
17	"The Legendre polynomials P_1(x) using an explicit
18	representation.")
19
20	(defmfun legendre-P2 (x)
21	"gsl_sf_legendre_P2_e" ((x :double) (ret sf-result))
22	:documentation ; FDL
23	"The Legendre polynomials P_2(x) using an explicit
24	representation.")
25
26	(defmfun legendre-P3 (x)
27	"gsl_sf_legendre_P3_e" ((x :double) (ret sf-result))
28	:documentation ; FDL
29	"The Legendre polynomials P_3(x) using an explicit
30	representation.")
31
32	(defmfun legendre-Pl (l x)
33	"gsl_sf_legendre_Pl_e" ((l :int) (x :double) (ret sf-result))
34	:documentation ; FDL
35	"The Legendre polynomial P_l(x) for a specific value of l,
36	x subject to l >= 0, \|x\| <= 1.")
37
38	(defmfun legendre-Pl-array (x array)
39	"gsl_sf_legendre_Pl_array"
40	(((1- (dim0 array)) :int) (x :double) ((gsl-array array) :pointer))
41	:documentation ; FDL
42	"Compute an array of Legendre polynomials
43	P_l(x) for l = 0, ..., length(array), \|x\| <= 1."
44	:invalidate (array))
45
46	(defmfun legendre-Pl-deriv-array (x array)
47	"gsl_sf_legendre_Pl_deriv_array"
48	(((1- (dim0 array)) :int) (x :double) ((gsl-array array) :pointer))
49	:documentation ; FDL
50	"Compute an array of Legendre polynomials derivatives
51	dP_l(x)/dx, for l = 0, ..., length(array), \|x\| <= 1."
52	:invalidate (array))
53
54	(defmfun legendre-Q0 (x)
55	"gsl_sf_legendre_Q0_e" ((x :double) (ret sf-result))
56	:documentation ; FDL
57	"The Legendre function Q_0(x) for x > -1,
58	x /= 1.")
59
60	(defmfun legendre-Q1 (x)
61	"gsl_sf_legendre_Q1_e" ((x :double) (ret sf-result))
62	:documentation ; FDL
63	"The Legendre function Q_1(x) for x > -1,
64	x /= 1.")
65
66	(defmfun legendre-Ql (l x)
67	"gsl_sf_legendre_Ql_e" ((l :int) (x :double) (ret sf-result))
68	:documentation ; FDL
69	"The Legendre function Q_l(x) for x > -1, x /= 1, l >= 0.")
70
71	;;;;****************************************************************************
72	;;;; Associated Legendre Polynomials and Spherical Harmonics
73	;;;;****************************************************************************
74
75	;;; FDL
76	;;; The following functions compute the associated Legendre Polynomials
77	;;; P_l^m(x). Note that this function grows combinatorially with
78	;;; l and can overflow for l larger than about 150. There is
79	;;; no trouble for small m, but overflow occurs when m and
80	;;; l are both large. Rather than allow overflows, these functions
81	;;; refuse to calculate P_l^m(x) and return :EOVRFLW when
82	;;; they can sense that l and m are too big.
83
84	;;; If you want to calculate a spherical harmonic, then do not use
85	;;; these functions. Instead use legendre-sphPlm below,
86	;;; which uses a similar recursion, but with the normalized functions.
87
88	(defmfun legendre-Plm (l m x)
89	"gsl_sf_legendre_Plm_e" ((l :int) (m :int) (x :double) (ret sf-result))
90	:documentation ; FDL
91	"The associated Legendre polynomial
92	P_l^m(x) for m >= 0, l >= m, \|x\| <= 1.")
93
94	(defmfun legendre-Plm-array (m x array)
95	"gsl_sf_legendre_Plm_array"
96	(((+ (dim0 array) m -1) :int) (m :int) (x :double)
97	((gsl-array array) :pointer))
98	:documentation ; FDL
99	"An array of Legendre polynomials
100	P_l^m(x), for m >= 0,
101	l = \|m\|, ..., \|m\|+length(array)-1} and \|x\| <= 1."
102	:invalidate (array))
103
104	(defmfun legendre-Plm-deriv-array (m x values derivatives)
105	"gsl_sf_legendre_Plm_deriv_array"
106	(((+ (dim0 values) m -1) :int) (m :int) (x :double)
107	((gsl-array values) :pointer) ((gsl-array derivatives) :pointer))
108	:documentation ; FDL
109	"An array of Legendre polynomials
110	values and derivatives dP_l^m(x)/dx for m >= 0,
111	l = \|m\|, ..., length(values) and \|x\| <= 1."
112	:invalidate (values derivatives))
113
114	(defmfun legendre-sphPlm (l m x)
115	"gsl_sf_legendre_sphPlm_e" ((l :int) (m :int) (x :double) (ret sf-result))
116	:documentation ; FDL
117	"The normalized associated Legendre polynomial
118	\sqrt{(2l+1)/(4\pi) \sqrt{(l-m)!/(l+m)!} P_l^m(x) suitable
119	for use in spherical harmonics. The parameters must satisfy
120	m >= 0, l >= m, \|x\| <= 1. These routines avoid the overflows
121	that occur for the standard normalization of P_l^m(x).")
122
123	(defmfun legendre-sphPlm-array (m x array)
124	"gsl_sf_legendre_sphPlm_array"
125	(((+ (dim0 array) m -1) :int) (m :int) (x :double)
126	((gsl-array array) :pointer))
127	:documentation ; FDL
128	"An array of normalized associated Legendre functions
129	\sqrt(2l+1)/(4\pi) \sqrt(l-m)!/(l+m)! P_l^m(x),
130	for m >= 0, l = \|m\|, ..., length(array)}, \|x\| <= 1.0."
131	:invalidate (array))
132
133	(defmfun legendre-sphPlm-deriv-array (m x values derivatives)
134	"gsl_sf_legendre_sphPlm_deriv_array"
135	(((+ (dim0 values) m -1) :int) (m :int) (x :double)
136	((gsl-array values) :pointer) ((gsl-array derivatives) :pointer))
137	:documentation ; FDL
138	"An array of normalized associated Legendre functions
139	values and derivatives for m >= 0,
140	l = \|m\|, ..., length(array)}, \|x\| <= 1.0."
141	:invalidate (values derivatives))
142
143	(defmfun legendre-array-size (lmax m)
144	"gsl_sf_legendre_array_size" ((lmax :int) (m :int))
145	:documentation ; FDL
146	"The size of result array needed for the array
147	versions of P_l^m(x), lmax - m + 1."
148	:c-return :int)
149
150	;;;;****************************************************************************
151	;;;; Conical Functions
152	;;;;****************************************************************************
153
154	;;; FDL
155	;;; The Conical Functions P^\mu_{-(1/2)+i\lambda}(x)} and
156	;;; Q^\mu_{-(1/2)+i\lambda}
157	;;; are described in Abramowitz & Stegun, Section 8.12.
158
159	(defmfun legendre-conicalP-half (lambda x)
160	"gsl_sf_conicalP_half_e" ((lambda :double) (x :double) (ret sf-result))
161	:documentation ; FDL
162	"The irregular Spherical Conical Function
163	P^{1/2}_{-1/2 + i \lambda}(x) for x > -1.")
164
165	(defmfun legendre-conicalP-mhalf (lambda x)
166	"gsl_sf_conicalP_mhalf_e" ((lambda :double) (x :double) (ret sf-result))
167	:documentation ; FDL
168	"The regular Spherical Conical Function
169	P^{-1/2}_{-1/2 + i \lambda}(x) for x > -1.")
170
171	(defmfun legendre-conicalP-0 (lambda x)
172	"gsl_sf_conicalP_0_e" ((lambda :double) (x :double) (ret sf-result))
173	:documentation ; FDL
174	"The conical function P^0_{-1/2 + i \lambda(x)} for x > -1.")
175
176	(defmfun legendre-conicalP-1 (lambda x)
177	"gsl_sf_conicalP_1_e" ((lambda :double) (x :double) (ret sf-result))
178	:documentation ; FDL
179	"The conical function
180	P^1_{-1/2 + i \lambda}(x)} for x > -1.")
181
182	(defmfun legendre-regular-spherical-conical (l lambda x)
183	"gsl_sf_conicalP_sph_reg_e"
184	((l :int) (lambda :double) (x :double) (ret sf-result))
185	:documentation ; FDL
186	"The Regular Spherical Conical Function
187	P^{-1/2-l}_{-1/2 + i \lambda}(x) for x > -1, l >= -1.")
188
189	(defmfun legendre-regular-cylindrical-conical (l lambda x)
190	"gsl_sf_conicalP_cyl_reg_e"
191	((l :int) (lambda :double) (x :double) (ret sf-result))
192	:documentation ; FDL
193	"The Regular Cylindrical Conical Function
194	P^{-m}_{-1/2 + i \lambda}(x) for x > -1, m >= -1.")
195
196	;;;;****************************************************************************
197	;;;; Radial Functions for Hyperbolic Space
198	;;;;****************************************************************************
199
200	;;; FDL
201	;;; The following spherical functions are specializations of Legendre
202	;;; functions which give the regular eigenfunctions of the Laplacian
203	;;; on a 3-dimensional hyperbolic space H3d. Of particular interest
204	;;; is the flat limit, \lambda \to \infty, \eta \to 0, \lambda\eta
205	;;; fixed.
206
207	(defmfun legendre-H3d-0 (lambda eta)
208	"gsl_sf_legendre_H3d_0_e"
209	((lambda :double) (eta :double) (ret sf-result))
210	:documentation ; FDL
211	"The zeroth radial eigenfunction of the Laplacian on the
212	3-dimensional hyperbolic space,
213	L^{H3d}_0(\lambda,\eta) := \sin(\lambda\eta)/(\lambda\sinh(\eta))
214	for \eta >= 0. In the flat limit this takes the form
215	L^{H3d}_0(\lambda,\eta) = j_0(\lambda\eta).")
216
217	(defmfun legendre-H3d-1 (lambda eta)
218	"gsl_sf_legendre_H3d_1_e"
219	((lambda :double) (eta :double) (ret sf-result))
220	:documentation ; FDL
221	"The first radial eigenfunction of the Laplacian on
222	the 3-dimensional hyperbolic space,
223	L^{H3d}_1(\lambda,\eta) := 1/\sqrt{\lambda^2 + 1}
224	\sin(\lambda \eta)/(\lambda \sinh(\eta)) (\coth(\eta) - \lambda \cot(\lambda\eta))}
225	for \eta >= 0. In the flat limit this takes the form
226	L^{H3d}_1(\lambda,\eta) = j_1(\lambda\eta)}.")
227
228	(defmfun legendre-H3d (l lambda eta)
229	"gsl_sf_legendre_H3d_e"
230	((l :int) (lambda :double) (eta :double) (ret sf-result))
231	:documentation ; FDL
232	"The l-th radial eigenfunction of the
233	Laplacian on the 3-dimensional hyperbolic space
234	\eta >= 0, l >= 0. In the flat limit this takes the form
235	L^{H3d}_l(\lambda,\eta) = j_l(\lambda\eta).")
236
237	(defmfun legendre-H3d-array (lambda eta array)
238	"gsl_sf_legendre_H3d_array"
239	(((1- (dim0 array)) :int) (lambda :double) (eta :double)
240	((gsl-array array) :pointer))
241	:invalidate (array)
242	:documentation ; FDL
243	"An array of radial eigenfunctions
244	L^{H3d}_l(\lambda, \eta) for 0 <= l <= length(array).")
245
246	;;; (defparameter hleg (make-data 'vector nil 3))
247	;;; (legendre-H3d-array 1.0d0 0.5d0 hleg)
248	;;; #<GSL-VECTOR #(0.9200342692589383d0 0.21694026450392123d0 0.047950660488307775d0) {C07CB51}>
249
250	;;;;****************************************************************************
251	;;;; Examples and unit test
252	;;;;****************************************************************************
253
254	#\|
255	(make-tests legendre
256	(legendre-P1 0.3d0)
257	(legendre-P2 0.3d0)
258	(legendre-P3 0.3d0)
259	(legendre-Pl -4 0.3d0)
260	(legendre-Pl 4 3.0d0)
261	(legendre-Pl 4 0.3d0)
262	(letm ((arr (vector-double-float 4)))
263	(legendre-Pl-array 0.5d0 arr)
264	(data arr))
265	(legendre-Q0 3.3d0)
266	(legendre-Q1 3.3d0)
267	(legendre-Ql 2 3.3d0)
268	(legendre-Plm 4 3 0.5d0)
269	(letm ((arr (vector-double-float 4)))
270	(legendre-Plm-array 2 0.5d0 arr)
271	(data arr))
272	(letm ((val (vector-double-float 4))
273	(deriv (vector-double-float 4)))
274	(legendre-Plm-deriv-array 2 0.5d0 val deriv)
275	(data deriv))
276	(legendre-sphplm 1200 1100 0.3d0)
277	(letm ((arr (vector-double-float 4)))
278	(legendre-sphPlm-array 4 0.5d0 arr)
279	(data arr))
280	(letm ((val (vector-double-float 4))
281	(deriv (vector-double-float 4)))
282	(legendre-sphPlm-deriv-array 4 0.5d0 val deriv)
283	(data deriv))
284	(legendre-conicalp-half 3.5d0 10.0d0)
285	(legendre-conicalp-mhalf 3.5d0 10.0d0)
286	(legendre-conicalp-0 3.5d0 10.0d0)
287	(legendre-conicalp-1 3.5d0 10.0d0)
288	(legendre-regular-spherical-conical 3 3.5d0 10.0d0)
289	(legendre-regular-cylindrical-conical 3 3.5d0 10.0d0)
290	(legendre-h3d-0 1.0d0 0.5d0)
291	(legendre-h3d-1 1.0d0 0.5d0)
292	(legendre-h3d 4 1.0d0 0.5d0)
293	(letm ((arr (vector-double-float 4)))
294	(legendre-h3d-array 1.0d0 0.5d0 arr)
295	(data arr)))
296	\|#
297
298	(LISP-UNIT:DEFINE-TEST LEGENDRE
299	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 0.3d0 0.0d0)
300	(MULTIPLE-VALUE-LIST
301	(LEGENDRE-P1
302	0.3d0)))
303	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
304	(LIST -0.365d0 2.8199664825478977d-16)
305	(MULTIPLE-VALUE-LIST (LEGENDRE-P2 0.3d0)))
306	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
307	(LIST -0.38249999999999995d0 1.9984014443252816d-16)
308	(MULTIPLE-VALUE-LIST (LEGENDRE-P3 0.3d0)))
309	(LISP-UNIT:ASSERT-ERROR 'GSL-CONDITION
310	(LEGENDRE-PL -4 0.3d0))
311	(LISP-UNIT:ASSERT-ERROR 'GSL-CONDITION
312	(LEGENDRE-PL 4 3.0d0))
313	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
314	(LIST 0.07293749999999999d0 5.668382430101814d-17)
315	(MULTIPLE-VALUE-LIST (LEGENDRE-PL 4 0.3d0)))
316	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
317	(LIST #(1.0d0 0.5d0 -0.125d0 -0.4375d0))
318	(MULTIPLE-VALUE-LIST
319	(LETM ((ARR (VECTOR-DOUBLE-FLOAT 4)))
320	(LEGENDRE-PL-ARRAY 0.5d0 ARR) (DATA ARR))))
321	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
322	(LIST 0.3128529498822064d0 1.3893461931245028d-16)
323	(MULTIPLE-VALUE-LIST (LEGENDRE-Q0 3.3d0)))
324	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
325	(LIST 0.03241473461128108d0 1.4395033881023292d-17)
326	(MULTIPLE-VALUE-LIST (LEGENDRE-Q1 3.3d0)))
327	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
328	(LIST 0.004026461384737812d0 1.788108054840004d-18)
329	(MULTIPLE-VALUE-LIST (LEGENDRE-QL 2 3.3d0)))
330	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
331	(LIST -34.099750274012266d0 3.0286662310541114d-14)
332	(MULTIPLE-VALUE-LIST (LEGENDRE-PLM 4 3 0.5d0)))
333	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
334	(LIST #(2.25d0 5.625d0 4.21875d0 -4.921875d0))
335	(MULTIPLE-VALUE-LIST
336	(LETM ((ARR (VECTOR-DOUBLE-FLOAT 4)))
337	(LEGENDRE-PLM-ARRAY 2 0.5d0 ARR) (DATA ARR))))
338	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
339	(LIST #(-3.0d0 3.75d0 33.75d0 55.78125d0))
340	(MULTIPLE-VALUE-LIST
341	(LETM
342	((VAL (VECTOR-DOUBLE-FLOAT 4)) (DERIV (VECTOR-DOUBLE-FLOAT 4)))
343	(LEGENDRE-PLM-DERIV-ARRAY 2 0.5d0 VAL DERIV)
344	(DATA DERIV))))
345	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
346	(LIST 0.30366280894310793d0 3.438761110552081d-15)
347	(MULTIPLE-VALUE-LIST
348	(LEGENDRE-SPHPLM 1200 1100 0.3d0)))
349	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
350	(LIST
351	#(0.24892463950030283d0 0.4127948151484927d0
352	0.35120655562190445d0 0.051599351893561574d0))
353	(MULTIPLE-VALUE-LIST
354	(LETM ((ARR (VECTOR-DOUBLE-FLOAT 4)))
355	(LEGENDRE-SPHPLM-ARRAY 4 0.5d0 ARR)
356	(DATA ARR))))
357	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
358	(LIST
359	#(-0.6637990386674741d0 -0.2751965434323283d0
360	1.2710332489173686d0 2.648766730536161d0))
361	(MULTIPLE-VALUE-LIST
362	(LETM
363	((VAL (VECTOR-DOUBLE-FLOAT 4)) (DERIV (VECTOR-DOUBLE-FLOAT 4)))
364	(LEGENDRE-SPHPLM-DERIV-ARRAY 4 0.5d0 VAL DERIV)
365	(DATA DERIV))))
366	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
367	(LIST -0.1255299048878925d0 1.3395992025650077d-15)
368	(MULTIPLE-VALUE-LIST
369	(LEGENDRE-CONICALP-HALF 3.5d0 10.0d0)))
370	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
371	(LIST -0.06274336292793128d0 2.504525777730328d-16)
372	(MULTIPLE-VALUE-LIST
373	(LEGENDRE-CONICALP-MHALF 3.5d0 10.0d0)))
374	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
375	(LIST -0.1316366189368757d0 3.0865532615549887d-15)
376	(MULTIPLE-VALUE-LIST
377	(LEGENDRE-CONICALP-0 3.5d0 10.0d0)))
378	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
379	(LIST 0.17407195560076358d0 4.024312727070929d-15)
380	(MULTIPLE-VALUE-LIST
381	(LEGENDRE-CONICALP-1 3.5d0 10.0d0)))
382	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
383	(LIST 8.980797952969897d-4 7.157154480497032d-17)
384	(MULTIPLE-VALUE-LIST
385	(LEGENDRE-REGULAR-SPHERICAL-CONICAL 3 3.5d0 10.0d0)))
386	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
387	(LIST 0.0023060250619859907d0 2.7345630541297237d-16)
388	(MULTIPLE-VALUE-LIST
389	(LEGENDRE-REGULAR-CYLINDRICAL-CONICAL 3 3.5d0
390	10.0d0)))
391	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
392	(LIST 0.9200342692589382d0 1.7018240558144874d-15)
393	(MULTIPLE-VALUE-LIST (LEGENDRE-H3D-0 1.0d0 0.5d0)))
394	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
395	(LIST 0.2169402645039212d0 1.418962379417407d-15)
396	(MULTIPLE-VALUE-LIST (LEGENDRE-H3D-1 1.0d0 0.5d0)))
397	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
398	(LIST 0.002400616238997978d0 4.3381136468045d-17)
399	(MULTIPLE-VALUE-LIST (LEGENDRE-H3D 4 1.0d0 0.5d0)))
400	(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
401	(LIST
402	#(0.9200342692589379d0 0.21694026450392115d0
403	0.04795066048830776d0 0.010663769096144337d0))
404	(MULTIPLE-VALUE-LIST
405	(LETM ((ARR (VECTOR-DOUBLE-FLOAT 4)))
406	(LEGENDRE-H3D-ARRAY 1.0d0 0.5d0 ARR)
407	(DATA ARR)))))

Note: See TracBrowser for help on using the browser.

Download in other formats:

Original Format