Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculate spherical harmonics
10 : !> \note
11 : !> Spherical Harmonics
12 : !> Numerical Stability up to L=15
13 : !> Accuracy > 1.E-12 up to L=15 tested
14 : !> Definition is consistent with orbital_transformation_matrices
15 : !> \todo
16 : !> 1) Check if this definition is consistent with the
17 : !> Slater-Koster module
18 : !> \par History
19 : !> JGH 28-Feb-2002 : Change of sign convention (-1^m)
20 : !> JGH 1-Mar-2002 : Clebsch-Gordon Coefficients
21 : !> - Clebsch-Gordon coefficients and Wigner 3-j symbols added as generic routines using the
22 : !> standard normalization (19.09.2022, MK)
23 : !> - Clebsch-Gordon coefficients (which were really Gaunt coefficients of the
24 : !> real spherical harmonics) moved to gaunt_coefficients, and the unused
25 : !> generic Clebsch-Gordon and Wigner 3-j routines dropped, in favour of the
26 : !> exact evaluation in libwignernj (09.2026, Susi Lehtola)
27 : !> \author JGH 6-Oct-2000, MK
28 : ! **************************************************************************************************
29 : MODULE spherical_harmonics
30 :
31 : USE kinds, ONLY: dp
32 : USE mathconstants, ONLY: fac,&
33 : pi
34 : #include "../base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 :
38 : PRIVATE
39 :
40 : PUBLIC :: y_lm, dy_lm
41 : PUBLIC :: legendre, dlegendre
42 :
43 : INTERFACE y_lm
44 : MODULE PROCEDURE rvy_lm, rry_lm, cvy_lm, ccy_lm
45 : END INTERFACE
46 :
47 : INTERFACE dy_lm
48 : MODULE PROCEDURE dry_lm, dcy_lm
49 : END INTERFACE
50 :
51 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'spherical_harmonics'
52 :
53 : CONTAINS
54 :
55 : ! Calculation of Spherical Harmonics
56 :
57 : ! **************************************************************************************************
58 : !> \brief ...
59 : !> \param r ...
60 : !> \param y ...
61 : !> \param l ...
62 : !> \param m ...
63 : ! **************************************************************************************************
64 67537 : SUBROUTINE rvy_lm(r, y, l, m)
65 : !
66 : ! Real Spherical Harmonics
67 : ! _ _
68 : ! | [(2l+1)(l-|m|)!] |1/2 m cos(m p) m>=0
69 : ! Y_lm ( t, p ) = |---------------------| P_l(cos(t))
70 : ! |[2Pi(1+d_m0)(l+|m|)!]| sin(|m| p) m<0
71 : !
72 : ! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1
73 : !
74 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r
75 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: y
76 : INTEGER, INTENT(IN) :: l, m
77 :
78 : INTEGER :: i
79 : REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, sp, t, z
80 :
81 67537 : SELECT CASE (l)
82 : CASE (:-1)
83 0 : CPABORT("Negative l value")
84 : CASE (0)
85 3153 : pf = SQRT(1.0_dp/(4.0_dp*pi))
86 3153 : IF (m /= 0) CPABORT("l = 0 and m value out of bounds")
87 170927 : y(:) = pf
88 : CASE (1)
89 24056 : SELECT CASE (m)
90 : CASE DEFAULT
91 0 : CPABORT("l = 1 and m value out of bounds")
92 : CASE (1)
93 3007 : pf = SQRT(3.0_dp/(4.0_dp*pi))
94 163881 : y(:) = pf*r(1, :)
95 : CASE (0)
96 3007 : pf = SQRT(3.0_dp/(4.0_dp*pi))
97 163881 : y(:) = pf*r(3, :)
98 : CASE (-1)
99 3007 : pf = SQRT(3.0_dp/(4.0_dp*pi))
100 169895 : y(:) = pf*r(2, :)
101 : END SELECT
102 : CASE (2)
103 30694 : SELECT CASE (m)
104 : CASE DEFAULT
105 0 : CPABORT("l = 2 and m value out of bounds")
106 : CASE (2)
107 3007 : pf = SQRT(15.0_dp/(16.0_dp*pi))
108 163881 : y(:) = pf*(r(1, :)*r(1, :) - r(2, :)*r(2, :))
109 : CASE (1)
110 3007 : pf = SQRT(15.0_dp/(4.0_dp*pi))
111 163881 : y(:) = pf*r(3, :)*r(1, :)
112 : CASE (0)
113 3007 : pf = SQRT(5.0_dp/(16.0_dp*pi))
114 163881 : y(:) = pf*(3.0_dp*r(3, :)*r(3, :) - 1.0_dp)
115 : CASE (-1)
116 3007 : pf = SQRT(15.0_dp/(4.0_dp*pi))
117 163881 : y(:) = pf*r(3, :)*r(2, :)
118 : CASE (-2)
119 3007 : pf = SQRT(15.0_dp/(16.0_dp*pi))
120 175909 : y(:) = pf*2.0_dp*r(1, :)*r(2, :)
121 : END SELECT
122 : CASE (3)
123 40328 : SELECT CASE (m)
124 : CASE DEFAULT
125 0 : CPABORT("l = 3 and m value out of bounds")
126 : CASE (3)
127 2237 : pf = SQRT(35.0_dp/(32.0_dp*pi))
128 125043 : y(:) = pf*r(1, :)*(r(1, :)**2 - 3.0_dp*r(2, :)**2)
129 : CASE (2)
130 2237 : pf = SQRT(105.0_dp/(16.0_dp*pi))
131 125043 : y(:) = pf*r(3, :)*(r(1, :)**2 - r(2, :)**2)
132 : CASE (1)
133 2237 : pf = SQRT(21.0_dp/(32.0_dp*pi))
134 125043 : y(:) = pf*r(1, :)*(5.0_dp*r(3, :)*r(3, :) - 1.0_dp)
135 : CASE (0)
136 2237 : pf = SQRT(7.0_dp/(16.0_dp*pi))
137 125043 : y(:) = pf*r(3, :)*(5.0_dp*r(3, :)*r(3, :) - 3.0_dp)
138 : CASE (-1)
139 2237 : pf = SQRT(21.0_dp/(32.0_dp*pi))
140 125043 : y(:) = pf*r(2, :)*(5.0_dp*r(3, :)*r(3, :) - 1.0_dp)
141 : CASE (-2)
142 2237 : pf = SQRT(105.0_dp/(16.0_dp*pi))
143 125043 : y(:) = pf*2.0_dp*r(1, :)*r(2, :)*r(3, :)
144 : CASE (-3)
145 2237 : pf = SQRT(35.0_dp/(32.0_dp*pi))
146 138465 : y(:) = pf*r(2, :)*(3.0_dp*r(1, :)**2 - r(2, :)**2)
147 : END SELECT
148 : CASE DEFAULT
149 24669 : IF (m < -l .OR. m > l) CPABORT("m value out of bounds")
150 24669 : lpm = fac(l + ABS(m))
151 24669 : lmm = fac(l - ABS(m))
152 24669 : IF (m == 0) THEN
153 : t = 4.0_dp*pi
154 : ELSE
155 22106 : t = 2.0_dp*pi
156 : END IF
157 24669 : IF (ABS(lpm) < EPSILON(1.0_dp)) THEN
158 0 : pf = REAL(2*l + 1, KIND=dp)/t
159 : ELSE
160 24669 : pf = (REAL(2*l + 1, KIND=dp)*lmm)/(t*lpm)
161 : END IF
162 24669 : pf = SQRT(pf)
163 1846756 : DO i = 1, SIZE(r, 2)
164 1754550 : z = r(3, i)
165 : ! plm = legendre ( z, l, m )
166 1754550 : plm = legendre(z, l, ABS(m))
167 1779219 : IF (m == 0) THEN
168 165530 : y(i) = pf*plm
169 : ELSE
170 1589020 : rxy = SQRT(r(1, i)**2 + r(2, i)**2)
171 1589020 : IF (rxy < EPSILON(1.0_dp)) THEN
172 44212 : y(i) = 0.0_dp
173 : ELSE
174 1544808 : cp = r(1, i)/rxy
175 1544808 : sp = r(2, i)/rxy
176 1544808 : IF (m > 0) THEN
177 772404 : y(i) = pf*plm*cosn(m, cp, sp)
178 : ELSE
179 772404 : y(i) = pf*plm*sinn(-m, cp, sp)
180 : END IF
181 : END IF
182 : END IF
183 : END DO
184 : END SELECT
185 :
186 67537 : END SUBROUTINE rvy_lm
187 :
188 : ! **************************************************************************************************
189 : !> \brief ...
190 : !> \param r ...
191 : !> \param y ...
192 : !> \param l ...
193 : !> \param m ...
194 : ! **************************************************************************************************
195 1605711 : SUBROUTINE rry_lm(r, y, l, m)
196 : !
197 : ! Real Spherical Harmonics
198 : ! _ _
199 : ! | [(2l+1)(l-|m|)!] |1/2 m cos(m p) m>=0
200 : ! Y_lm ( t, p ) = |---------------------| P_l(cos(t))
201 : ! |[2Pi(1+d_m0)(l+|m|)!]| sin(|m| p) m<0
202 : !
203 : ! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1
204 : !
205 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r
206 : REAL(KIND=dp), INTENT(OUT) :: y
207 : INTEGER, INTENT(IN) :: l, m
208 :
209 : REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, sp, t, z
210 :
211 1605711 : SELECT CASE (l)
212 : CASE (:-1)
213 0 : CPABORT("Negative l value")
214 : CASE (0)
215 178735 : pf = SQRT(1.0_dp/(4.0_dp*pi))
216 178735 : IF (m /= 0) CPABORT("l = 0 and m value out of bounds")
217 178735 : y = pf
218 : CASE (1)
219 1219984 : SELECT CASE (m)
220 : CASE DEFAULT
221 0 : CPABORT("l = 1 and m value out of bounds")
222 : CASE (1)
223 218467 : pf = SQRT(3.0_dp/(4.0_dp*pi))
224 218467 : y = pf*r(1)
225 : CASE (0)
226 77351 : pf = SQRT(3.0_dp/(4.0_dp*pi))
227 77351 : y = pf*r(3)
228 : CASE (-1)
229 218467 : pf = SQRT(3.0_dp/(4.0_dp*pi))
230 514285 : y = pf*r(2)
231 : END SELECT
232 : CASE (2)
233 801210 : SELECT CASE (m)
234 : CASE DEFAULT
235 0 : CPABORT("l = 2 and m value out of bounds")
236 : CASE (2)
237 157087 : pf = SQRT(15.0_dp/(16.0_dp*pi))
238 157087 : y = pf*(r(1)*r(1) - r(2)*r(2))
239 : CASE (1)
240 157087 : pf = SQRT(15.0_dp/(4.0_dp*pi))
241 157087 : y = pf*r(3)*r(1)
242 : CASE (0)
243 77351 : pf = SQRT(5.0_dp/(16.0_dp*pi))
244 77351 : y = pf*(3.0_dp*r(3)*r(3) - 1.0_dp)
245 : CASE (-1)
246 157087 : pf = SQRT(15.0_dp/(4.0_dp*pi))
247 157087 : y = pf*r(3)*r(2)
248 : CASE (-2)
249 157087 : pf = SQRT(15.0_dp/(16.0_dp*pi))
250 705699 : y = pf*2.0_dp*r(1)*r(2)
251 : END SELECT
252 : CASE (3)
253 206992 : SELECT CASE (m)
254 : CASE DEFAULT
255 0 : CPABORT("l = 3 and m value out of bounds")
256 : CASE (3)
257 14981 : pf = SQRT(35.0_dp/(32.0_dp*pi))
258 14981 : y = pf*r(1)*(r(1)**2 - 3.0_dp*r(2)**2)
259 : CASE (2)
260 14981 : pf = SQRT(105.0_dp/(16.0_dp*pi))
261 14981 : y = pf*r(3)*(r(1)**2 - r(2)**2)
262 : CASE (1)
263 14981 : pf = SQRT(21.0_dp/(32.0_dp*pi))
264 14981 : y = pf*r(1)*(5.0_dp*r(3)*r(3) - 1.0_dp)
265 : CASE (0)
266 5625 : pf = SQRT(7.0_dp/(16.0_dp*pi))
267 5625 : y = pf*r(3)*(5.0_dp*r(3)*r(3) - 3.0_dp)
268 : CASE (-1)
269 14981 : pf = SQRT(21.0_dp/(32.0_dp*pi))
270 14981 : y = pf*r(2)*(5.0_dp*r(3)*r(3) - 1.0_dp)
271 : CASE (-2)
272 14981 : pf = SQRT(105.0_dp/(16.0_dp*pi))
273 14981 : y = pf*2.0_dp*r(1)*r(2)*r(3)
274 : CASE (-3)
275 14981 : pf = SQRT(35.0_dp/(32.0_dp*pi))
276 95511 : y = pf*r(2)*(3.0_dp*r(1)**2 - r(2)**2)
277 : END SELECT
278 : CASE DEFAULT
279 111481 : IF (m < -l .OR. m > l) CPABORT("m value out of bounds")
280 111481 : lpm = fac(l + ABS(m))
281 111481 : lmm = fac(l - ABS(m))
282 111481 : IF (m == 0) THEN
283 : t = 4.0_dp*pi
284 : ELSE
285 105856 : t = 2.0_dp*pi
286 : END IF
287 111481 : IF (ABS(lpm) < EPSILON(1.0_dp)) THEN
288 0 : pf = REAL(2*l + 1, KIND=dp)/t
289 : ELSE
290 111481 : pf = (REAL(2*l + 1, KIND=dp)*lmm)/(t*lpm)
291 : END IF
292 111481 : pf = SQRT(pf)
293 111481 : z = r(3)
294 111481 : plm = legendre(z, l, m)
295 1717192 : IF (m == 0) THEN
296 5625 : y = pf*plm
297 : ELSE
298 105856 : rxy = SQRT(r(1)**2 + r(2)**2)
299 105856 : IF (rxy < EPSILON(1.0_dp)) THEN
300 7368 : y = 0.0_dp
301 : ELSE
302 98488 : cp = r(1)/rxy
303 98488 : sp = r(2)/rxy
304 98488 : IF (m > 0) THEN
305 49244 : y = pf*plm*cosn(m, cp, sp)
306 : ELSE
307 49244 : y = pf*plm*sinn(-m, cp, sp)
308 : END IF
309 : END IF
310 : END IF
311 : END SELECT
312 :
313 1605711 : END SUBROUTINE rry_lm
314 :
315 : ! **************************************************************************************************
316 : !> \brief ...
317 : !> \param r ...
318 : !> \param y ...
319 : !> \param l ...
320 : !> \param m ...
321 : ! **************************************************************************************************
322 225 : SUBROUTINE cvy_lm(r, y, l, m)
323 : !
324 : ! Complex Spherical Harmonics
325 : ! _ _
326 : ! | [(2l+1)(l-|m|)!] |1/2 m
327 : ! Y_lm ( t, p ) = |------------------| P_l(cos(t)) Exp[ i m p ]
328 : ! | [4Pi(l+|m|)!]| |
329 : !
330 : ! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1
331 : !
332 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r
333 : COMPLEX(KIND=dp), DIMENSION(:), INTENT(OUT) :: y
334 : INTEGER, INTENT(IN) :: l, m
335 :
336 : INTEGER :: i, n
337 : REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, sp, t, ym, &
338 : yp, z
339 :
340 225 : n = SIZE(r, 2)
341 225 : SELECT CASE (l)
342 : CASE (:-1)
343 0 : CPABORT("Negative l value")
344 : CASE (0)
345 1 : pf = SQRT(1.0_dp/(4.0_dp*pi))
346 1 : IF (m /= 0) CPABORT("l = 0 and m value out of bounds")
347 591 : y(:) = pf
348 : CASE (1)
349 8 : SELECT CASE (m)
350 : CASE DEFAULT
351 0 : CPABORT("l = 1 and m value out of bounds")
352 : CASE (1)
353 591 : pf = SQRT(3.0_dp/(8.0_dp*pi))
354 591 : DO i = 1, n
355 590 : yp = r(1, i)
356 590 : ym = r(2, i)
357 591 : y(i) = pf*CMPLX(yp, ym, KIND=dp)
358 : END DO
359 : CASE (0)
360 1 : pf = SQRT(3.0_dp/(4.0_dp*pi))
361 591 : y(:) = pf*r(3, :)
362 : CASE (-1)
363 591 : pf = SQRT(3.0_dp/(8.0_dp*pi))
364 593 : DO i = 1, n
365 590 : yp = r(1, i)
366 590 : ym = r(2, i)
367 591 : y(i) = pf*CMPLX(yp, -ym, KIND=dp)
368 : END DO
369 : END SELECT
370 : CASE (2)
371 12 : SELECT CASE (m)
372 : CASE DEFAULT
373 0 : CPABORT("l = 2 and m value out of bounds")
374 : CASE (2)
375 591 : pf = SQRT(15.0_dp/(32.0_dp*pi))
376 591 : DO i = 1, n
377 590 : yp = (r(1, i)*r(1, i) - r(2, i)*r(2, i))
378 590 : ym = 2.0_dp*r(1, i)*r(2, i)
379 591 : y(i) = pf*CMPLX(yp, ym, KIND=dp)
380 : END DO
381 : CASE (1)
382 591 : pf = SQRT(15.0_dp/(8.0_dp*pi))
383 591 : DO i = 1, n
384 590 : yp = r(3, i)*r(1, i)
385 590 : ym = r(3, i)*r(2, i)
386 591 : y(i) = pf*CMPLX(yp, ym, KIND=dp)
387 : END DO
388 : CASE (0)
389 1 : pf = SQRT(5.0_dp/(16.0_dp*pi))
390 591 : y(:) = pf*(3.0_dp*r(3, :)*r(3, :) - 1.0_dp)
391 : CASE (-1)
392 591 : pf = SQRT(15.0_dp/(8.0_dp*pi))
393 591 : DO i = 1, n
394 590 : yp = r(3, i)*r(1, i)
395 590 : ym = r(3, i)*r(2, i)
396 591 : y(i) = pf*CMPLX(yp, -ym, KIND=dp)
397 : END DO
398 : CASE (-2)
399 591 : pf = SQRT(15.0_dp/(32.0_dp*pi))
400 595 : DO i = 1, n
401 590 : yp = (r(1, i)*r(1, i) - r(2, i)*r(2, i))
402 590 : ym = 2.0_dp*r(1, i)*r(2, i)
403 591 : y(i) = pf*CMPLX(yp, -ym, KIND=dp)
404 : END DO
405 : END SELECT
406 : CASE (3)
407 216 : SELECT CASE (m)
408 : CASE DEFAULT
409 0 : CPABORT("l = 3 and m value out of bounds")
410 : CASE (3)
411 591 : pf = SQRT(35.0_dp/(64.0_dp*pi))
412 591 : DO i = 1, n
413 590 : yp = r(1, i)*(r(1, i)**2 - 3.0_dp*r(2, i)**2)
414 590 : ym = r(2, i)*(3.0_dp*r(1, i)**2 - r(2, i)**2)
415 591 : y(i) = pf*CMPLX(yp, ym, KIND=dp)
416 : END DO
417 : CASE (2)
418 591 : pf = SQRT(105.0_dp/(32.0_dp*pi))
419 591 : DO i = 1, n
420 590 : yp = r(3, i)*(r(1, i)**2 - r(2, i)**2)
421 590 : ym = 2.0_dp*r(1, i)*r(2, i)*r(3, i)
422 591 : y(i) = pf*CMPLX(yp, ym, KIND=dp)
423 : END DO
424 : CASE (1)
425 591 : pf = SQRT(21.0_dp/(64.0_dp*pi))
426 591 : DO i = 1, n
427 590 : yp = r(1, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
428 590 : ym = r(2, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
429 591 : y(i) = pf*CMPLX(yp, ym, KIND=dp)
430 : END DO
431 : CASE (0)
432 1 : pf = SQRT(7.0_dp/(16.0_dp*pi))
433 591 : y(:) = pf*r(3, :)*(5.0_dp*r(3, :)*r(3, :) - 3.0_dp)
434 : CASE (-1)
435 591 : pf = SQRT(21.0_dp/(64.0_dp*pi))
436 591 : DO i = 1, n
437 590 : yp = r(1, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
438 590 : ym = r(2, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
439 591 : y(i) = pf*CMPLX(yp, -ym, KIND=dp)
440 : END DO
441 : CASE (-2)
442 591 : pf = SQRT(105.0_dp/(32.0_dp*pi))
443 591 : DO i = 1, n
444 590 : yp = r(3, i)*(r(1, i)**2 - r(2, i)**2)
445 590 : ym = 2.0_dp*r(1, i)*r(2, i)*r(3, i)
446 591 : y(i) = pf*CMPLX(yp, -ym, KIND=dp)
447 : END DO
448 : CASE (-3)
449 591 : pf = SQRT(35.0_dp/(64.0_dp*pi))
450 597 : DO i = 1, n
451 590 : yp = r(1, i)*(r(1, i)**2 - 3.0_dp*r(2, i)**2)
452 590 : ym = r(2, i)*(3.0_dp*r(1, i)**2 - r(2, i)**2)
453 591 : y(i) = pf*CMPLX(yp, -ym, KIND=dp)
454 : END DO
455 : END SELECT
456 : CASE DEFAULT
457 209 : IF (m < -l .OR. m > l) CPABORT("m value out of bounds")
458 209 : lpm = fac(l + ABS(m))
459 209 : lmm = fac(l - ABS(m))
460 209 : t = 4.0_dp*pi
461 209 : IF (ABS(lpm) < EPSILON(1.0_dp)) THEN
462 0 : pf = REAL(2*l + 1, KIND=dp)/t
463 : ELSE
464 209 : pf = (REAL(2*l + 1, KIND=dp)*lmm)/(t*lpm)
465 : END IF
466 209 : pf = SQRT(pf)
467 123744 : DO i = 1, n
468 123310 : z = r(3, i)
469 123310 : plm = legendre(z, l, m)
470 123519 : IF (m == 0) THEN
471 6490 : y(i) = pf*plm
472 : ELSE
473 116820 : rxy = SQRT(r(1, i)**2 + r(2, i)**2)
474 116820 : IF (rxy < EPSILON(1.0_dp)) THEN
475 396 : y(i) = 0.0_dp
476 : ELSE
477 116424 : cp = r(1, i)/rxy
478 116424 : sp = r(2, i)/rxy
479 116424 : IF (m > 0) THEN
480 58212 : yp = cosn(m, cp, sp)
481 58212 : ym = sinn(m, cp, sp)
482 58212 : y(i) = pf*plm*CMPLX(yp, ym, KIND=dp)
483 : ELSE
484 58212 : yp = cosn(-m, cp, sp)
485 58212 : ym = sinn(-m, cp, sp)
486 58212 : y(i) = pf*plm*CMPLX(yp, -ym, KIND=dp)
487 : END IF
488 : END IF
489 : END IF
490 : END DO
491 : END SELECT
492 :
493 225 : END SUBROUTINE cvy_lm
494 :
495 : ! **************************************************************************************************
496 : !> \brief ...
497 : !> \param r ...
498 : !> \param y ...
499 : !> \param l ...
500 : !> \param m ...
501 : ! **************************************************************************************************
502 0 : SUBROUTINE ccy_lm(r, y, l, m)
503 : !
504 : ! Complex Spherical Harmonics
505 : ! _ _
506 : ! | [(2l+1)(l-|m|)!] |1/2 m
507 : ! Y_lm ( t, p ) = |------------------| P_l(cos(t)) Exp[ i m p ]
508 : ! | [4Pi(l+|m|)!]| |
509 : !
510 : ! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1
511 : !
512 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r
513 : COMPLEX(KIND=dp), INTENT(OUT) :: y
514 : INTEGER, INTENT(IN) :: l, m
515 :
516 : REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, sp, t, ym, &
517 : yp, z
518 :
519 0 : SELECT CASE (l)
520 : CASE (:-1)
521 0 : CPABORT("Negative l value")
522 : CASE (0)
523 0 : pf = SQRT(1.0_dp/(4.0_dp*pi))
524 0 : IF (m /= 0) CPABORT("l = 0 and m value out of bounds")
525 0 : y = pf
526 : CASE (1)
527 0 : SELECT CASE (m)
528 : CASE DEFAULT
529 0 : CPABORT("l = 1 and m value out of bounds")
530 : CASE (1)
531 0 : pf = SQRT(3.0_dp/(8.0_dp*pi))
532 0 : yp = r(1)
533 0 : ym = r(2)
534 0 : y = pf*CMPLX(yp, ym, KIND=dp)
535 : CASE (0)
536 0 : pf = SQRT(3.0_dp/(4.0_dp*pi))
537 0 : y = pf*r(3)
538 : CASE (-1)
539 0 : pf = SQRT(3.0_dp/(8.0_dp*pi))
540 0 : yp = r(1)
541 0 : ym = r(2)
542 0 : y = pf*CMPLX(yp, -ym, KIND=dp)
543 : END SELECT
544 : CASE (2)
545 0 : SELECT CASE (m)
546 : CASE DEFAULT
547 0 : CPABORT("l = 2 and m value out of bounds")
548 : CASE (2)
549 0 : pf = SQRT(15.0_dp/(32.0_dp*pi))
550 0 : yp = (r(1)*r(1) - r(2)*r(2))
551 0 : ym = 2.0_dp*r(1)*r(2)
552 0 : y = pf*CMPLX(yp, ym, KIND=dp)
553 : CASE (1)
554 0 : pf = SQRT(15.0_dp/(8.0_dp*pi))
555 0 : yp = r(3)*r(1)
556 0 : ym = r(3)*r(2)
557 0 : y = pf*CMPLX(yp, ym, KIND=dp)
558 : CASE (0)
559 0 : pf = SQRT(5.0_dp/(16.0_dp*pi))
560 0 : y = pf*(3.0_dp*r(3)*r(3) - 1.0_dp)
561 : CASE (-1)
562 0 : pf = SQRT(15.0_dp/(8.0_dp*pi))
563 0 : yp = r(3)*r(1)
564 0 : ym = r(3)*r(2)
565 0 : y = pf*CMPLX(yp, -ym, KIND=dp)
566 : CASE (-2)
567 0 : pf = SQRT(15.0_dp/(32.0_dp*pi))
568 0 : yp = (r(1)*r(1) - r(2)*r(2))
569 0 : ym = 2.0_dp*r(1)*r(2)
570 0 : y = pf*CMPLX(yp, -ym, KIND=dp)
571 : END SELECT
572 : CASE (3)
573 0 : SELECT CASE (m)
574 : CASE DEFAULT
575 0 : CPABORT("l = 3 and m value out of bounds")
576 : CASE (3)
577 0 : pf = SQRT(35.0_dp/(64.0_dp*pi))
578 0 : yp = r(1)*(r(1)**2 - 3.0_dp*r(2)**2)
579 0 : ym = r(2)*(3.0_dp*r(1)**2 - r(2)**2)
580 0 : y = pf*CMPLX(yp, ym, KIND=dp)
581 : CASE (2)
582 0 : pf = SQRT(105.0_dp/(32.0_dp*pi))
583 0 : yp = r(3)*(r(1)**2 - r(2)**2)
584 0 : ym = 2.0_dp*r(1)*r(2)*r(3)
585 0 : y = pf*CMPLX(yp, ym, KIND=dp)
586 : CASE (1)
587 0 : pf = SQRT(21.0_dp/(64.0_dp*pi))
588 0 : yp = r(1)*(5.0_dp*r(3)*r(3) - 1.0_dp)
589 0 : ym = r(2)*(5.0_dp*r(3)*r(3) - 1.0_dp)
590 0 : y = pf*CMPLX(yp, ym, KIND=dp)
591 : CASE (0)
592 0 : pf = SQRT(7.0_dp/(16.0_dp*pi))
593 0 : y = pf*r(3)*(5.0_dp*r(3)*r(3) - 3.0_dp)
594 : CASE (-1)
595 0 : pf = SQRT(21.0_dp/(64.0_dp*pi))
596 0 : yp = r(1)*(5.0_dp*r(3)*r(3) - 1.0_dp)
597 0 : ym = r(2)*(5.0_dp*r(3)*r(3) - 1.0_dp)
598 0 : y = pf*CMPLX(yp, -ym, KIND=dp)
599 : CASE (-2)
600 0 : pf = SQRT(105.0_dp/(32.0_dp*pi))
601 0 : yp = r(3)*(r(1)**2 - r(2)**2)
602 0 : ym = 2.0_dp*r(1)*r(2)*r(3)
603 0 : y = pf*CMPLX(yp, -ym, KIND=dp)
604 : CASE (-3)
605 0 : pf = SQRT(35.0_dp/(64.0_dp*pi))
606 0 : yp = r(1)*(r(1)**2 - 3.0_dp*r(2)**2)
607 0 : ym = r(2)*(3.0_dp*r(1)**2 - r(2)**2)
608 0 : y = pf*CMPLX(yp, -ym, KIND=dp)
609 : END SELECT
610 : CASE DEFAULT
611 0 : IF (m < -l .OR. m > l) CPABORT("m value out of bounds")
612 0 : lpm = fac(l + ABS(m))
613 0 : lmm = fac(l - ABS(m))
614 0 : t = 4.0_dp*pi
615 0 : IF (ABS(lpm) < EPSILON(1.0_dp)) THEN
616 0 : pf = REAL(2*l + 1, KIND=dp)/t
617 : ELSE
618 0 : pf = (REAL(2*l + 1, KIND=dp)*lmm)/(t*lpm)
619 : END IF
620 0 : pf = SQRT(pf)
621 0 : z = r(3)
622 0 : plm = legendre(z, l, m)
623 0 : IF (m == 0) THEN
624 0 : y = pf*plm
625 : ELSE
626 0 : rxy = SQRT(r(1)**2 + r(2)**2)
627 0 : IF (rxy < EPSILON(1.0_dp)) THEN
628 0 : y = 0.0_dp
629 : ELSE
630 0 : cp = r(1)/rxy
631 0 : sp = r(2)/rxy
632 0 : IF (m > 0) THEN
633 0 : yp = cosn(m, cp, sp)
634 0 : ym = sinn(m, cp, sp)
635 0 : y = pf*plm*CMPLX(yp, ym, KIND=dp)
636 : ELSE
637 0 : yp = cosn(-m, cp, sp)
638 0 : ym = sinn(-m, cp, sp)
639 0 : y = pf*plm*CMPLX(yp, -ym, KIND=dp)
640 : END IF
641 : END IF
642 : END IF
643 : END SELECT
644 :
645 0 : END SUBROUTINE ccy_lm
646 :
647 : ! Calculation of derivatives of Spherical Harmonics
648 :
649 : ! **************************************************************************************************
650 : !> \brief ...
651 : !> \param c ...
652 : !> \param dy ...
653 : !> \param l ...
654 : !> \param m ...
655 : ! **************************************************************************************************
656 1122344 : SUBROUTINE dry_lm(c, dy, l, m)
657 : !
658 : ! Real Spherical Harmonics
659 : ! _ _
660 : ! | [(2l+1)(l-|m|)!] |1/2 m cos(m p) m>=0
661 : ! Y_lm ( t, p ) = |---------------------| P_l(cos(t))
662 : ! |[2Pi(1+d_m0)(l+|m|)!]| sin(|m| p) m<0
663 : !
664 : ! Input: c == (t,p)
665 : ! Output: dy == (dy/dt, dy/dp)
666 : !
667 : ! x == sin(t)*cos(p)
668 : ! y == sin(t)*sin(p)
669 : ! z == cos(t)
670 : !
671 :
672 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: c
673 : REAL(KIND=dp), DIMENSION(2), INTENT(OUT) :: dy
674 : INTEGER, INTENT(IN) :: l, m
675 :
676 : REAL(KIND=dp) :: cp, ct, dplm, lmm, lpm, p, pf, rxy, sp, &
677 : st, t, tt, y, z
678 : REAL(KIND=dp), DIMENSION(3) :: r
679 :
680 1122344 : t = c(1)
681 1122344 : ct = COS(t)
682 1122344 : st = SIN(t)
683 1122344 : p = c(2)
684 1122344 : cp = COS(p)
685 1122344 : sp = SIN(p)
686 1122344 : r(1) = st*cp
687 1122344 : r(2) = st*sp
688 1122344 : r(3) = ct
689 :
690 : ! dY/dp
691 1122344 : IF (m == 0) THEN
692 404176 : dy(2) = 0.0_dp
693 : ELSE
694 718168 : CALL y_lm(r, y, l, -m)
695 718168 : dy(2) = -REAL(m, KIND=dp)*y
696 : END IF
697 :
698 : ! dY/dt
699 1122344 : SELECT CASE (l)
700 : CASE (:-1)
701 0 : CPABORT("Negative l value")
702 : CASE (0)
703 167184 : IF (m /= 0) CPABORT("l = 0 and m value out of bounds")
704 167184 : dy(1) = 0.0_dp
705 : CASE (1)
706 398680 : SELECT CASE (m)
707 : CASE DEFAULT
708 0 : CPABORT("l = 1 and m value out of bounds")
709 : CASE (1)
710 141116 : pf = SQRT(3.0_dp/(4.0_dp*pi))
711 141116 : dy(1) = pf*ct*cp
712 : CASE (0)
713 141116 : pf = SQRT(3.0_dp/(4.0_dp*pi))
714 141116 : dy(1) = -pf*st
715 : CASE (-1)
716 141116 : pf = SQRT(3.0_dp/(4.0_dp*pi))
717 423348 : dy(1) = pf*ct*sp
718 : END SELECT
719 : CASE (2)
720 65492 : SELECT CASE (m)
721 : CASE DEFAULT
722 0 : CPABORT("l = 2 and m value out of bounds")
723 : CASE (2)
724 79736 : pf = SQRT(15.0_dp/(16.0_dp*pi))
725 79736 : dy(1) = pf*2.0_dp*st*ct*COS(2._dp*p)
726 : CASE (1)
727 79736 : pf = SQRT(15.0_dp/(4.0_dp*pi))
728 79736 : dy(1) = pf*cp*(ct*ct - st*st)
729 : CASE (0)
730 79736 : pf = SQRT(5.0_dp/(16.0_dp*pi))
731 79736 : dy(1) = -pf*6.0_dp*ct*st
732 : CASE (-1)
733 79736 : pf = SQRT(15.0_dp/(4.0_dp*pi))
734 79736 : dy(1) = pf*sp*(ct*ct - st*st)
735 : CASE (-2)
736 79736 : pf = SQRT(15.0_dp/(16.0_dp*pi))
737 398680 : dy(1) = pf*2.0_dp*st*ct*SIN(2._dp*p)
738 : END SELECT
739 : CASE (3)
740 67640 : SELECT CASE (m)
741 : CASE DEFAULT
742 0 : CPABORT("l = 3 and m value out of bounds")
743 : CASE (3)
744 9356 : pf = SQRT(35.0_dp/(32.0_dp*pi))
745 9356 : dy(1) = pf*3.0_dp*COS(3._dp*p)*ct*st*st
746 : CASE (2)
747 9356 : pf = SQRT(105.0_dp/(16.0_dp*pi))
748 9356 : dy(1) = pf*2.0_dp*COS(2._dp*p)*ct*st
749 : CASE (1)
750 9356 : pf = SQRT(21.0_dp/(32.0_dp*pi))
751 9356 : dy(1) = pf*cp*(ct*(5.0_dp*ct - 1.0_dp) - 5.0_dp*st*st)
752 : CASE (0)
753 9356 : pf = SQRT(7.0_dp/(16.0_dp*pi))
754 9356 : dy(1) = pf*r(3)*(3.0_dp - 15.0_dp*ct*ct)*st
755 : CASE (-1)
756 9356 : pf = SQRT(21.0_dp/(32.0_dp*pi))
757 9356 : dy(1) = pf*sp*(ct*(5.0_dp*ct - 1.0_dp) - 5.0_dp*st*st)
758 : CASE (-2)
759 9356 : pf = SQRT(105.0_dp/(16.0_dp*pi))
760 9356 : dy(1) = pf*2.0_dp*SIN(2._dp*p)*ct*st
761 : CASE (-3)
762 9356 : pf = SQRT(35.0_dp/(32.0_dp*pi))
763 65492 : dy(1) = pf*3.0_dp*SIN(3._dp*p)*ct*st*st
764 : END SELECT
765 : CASE DEFAULT
766 67640 : IF (m < -l .OR. m > l) CPABORT("m value out of bounds")
767 67640 : lpm = fac(l + ABS(m))
768 67640 : lmm = fac(l - ABS(m))
769 : IF (m == 0) THEN
770 : tt = 4.0_dp*pi
771 : ELSE
772 : tt = 2.0_dp*pi
773 : END IF
774 : IF (ABS(lpm) < EPSILON(1.0_dp)) THEN
775 : pf = REAL(2*l + 1, KIND=dp)/tt
776 : ELSE
777 67640 : pf = (REAL(2*l + 1, KIND=dp)*lmm)/(tt*lpm)
778 : END IF
779 67640 : pf = SQRT(pf)
780 67640 : z = ct
781 67640 : dplm = dlegendre(z, l, m)
782 1189984 : IF (m == 0) THEN
783 : y = pf*dplm
784 : ELSE
785 60856 : rxy = SQRT(r(1)**2 + r(2)**2)
786 : IF (rxy < EPSILON(1.0_dp)) THEN
787 : y = 0.0_dp
788 : ELSE
789 : IF (m > 0) THEN
790 : y = pf*dplm*cosn(m, cp, sp)
791 : ELSE
792 : y = pf*dplm*sinn(-m, cp, sp)
793 : END IF
794 : END IF
795 : END IF
796 : END SELECT
797 :
798 1122344 : END SUBROUTINE dry_lm
799 :
800 : ! **************************************************************************************************
801 : !> \brief ...
802 : !> \param c ...
803 : !> \param dy ...
804 : !> \param l ...
805 : !> \param m ...
806 : ! **************************************************************************************************
807 0 : SUBROUTINE dcy_lm(c, dy, l, m)
808 : !
809 : ! Complex Spherical Harmonics
810 : ! _ _
811 : ! | [(2l+1)(l-|m|)!] |1/2 m
812 : ! Y_lm ( t, p ) = |------------------| P_l(cos(t)) Exp[ i m p ]
813 : ! | [4Pi(l+|m|)!]| |
814 : !
815 : ! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1
816 : !
817 : ! Input: c == (t,p)
818 : ! Output: dy == (dy/dt, dy/dp)
819 : !
820 : ! x == sin(t)*cos(p)
821 : ! y == sin(t)*sin(p)
822 : ! z == cos(t)
823 : !
824 :
825 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: c
826 : COMPLEX(KIND=dp), DIMENSION(2), INTENT(OUT) :: dy
827 : INTEGER, INTENT(IN) :: l, m
828 :
829 : REAL(KIND=dp), DIMENSION(2) :: dd
830 :
831 0 : CPABORT("Not implemented")
832 0 : CALL dry_lm(c, dd, l, m)
833 0 : dy(1) = CMPLX(dd(1), 0.0_dp, KIND=dp)
834 0 : dy(2) = CMPLX(dd(2), 0.0_dp, KIND=dp)
835 :
836 0 : END SUBROUTINE dcy_lm
837 :
838 : ! **************************************************************************************************
839 : !> \brief ...
840 : !> \param x ...
841 : !> \param l ...
842 : !> \param m ...
843 : !> \return ...
844 : ! **************************************************************************************************
845 2481513 : FUNCTION legendre(x, l, m) RESULT(plm)
846 :
847 : REAL(KIND=dp), INTENT(IN) :: x
848 : INTEGER, INTENT(IN) :: l, m
849 : REAL(KIND=dp) :: plm
850 :
851 : INTEGER :: il, im, mm
852 : REAL(KIND=dp) :: fact, pll, pmm, pmmp1, somx2
853 :
854 2481513 : IF (ABS(x) > 1.0_dp) CPABORT("x value > 1")
855 2481513 : SELECT CASE (l)
856 : CASE (:-1)
857 0 : CPABORT("Negative l value")
858 : CASE (0)
859 198 : plm = 1.0_dp
860 : CASE (1)
861 396 : SELECT CASE (ABS(m))
862 : CASE DEFAULT
863 0 : CPABORT("l = 1 and m value out of bounds")
864 : CASE (1)
865 0 : plm = SQRT(1.0_dp - x*x)
866 : CASE (0)
867 198 : plm = x
868 : END SELECT
869 : CASE (2)
870 396 : SELECT CASE (ABS(m))
871 : CASE DEFAULT
872 0 : CPABORT("l = 2 and m value out of bounds")
873 : CASE (2)
874 0 : plm = 3.0_dp*(1.0_dp - x*x)
875 : CASE (1)
876 0 : plm = 3.0_dp*x*SQRT(1.0_dp - x*x)
877 : CASE (0)
878 198 : plm = 1.5_dp*x*x - 0.5_dp
879 : END SELECT
880 : CASE (3)
881 1189323 : SELECT CASE (ABS(m))
882 : CASE DEFAULT
883 0 : CPABORT("l = 3 and m value out of bounds")
884 : CASE (3)
885 0 : plm = 15.0_dp*(1.0_dp - x*x)**1.5_dp
886 : CASE (2)
887 0 : plm = 15.0_dp*x*(1.0_dp - x*x)
888 : CASE (1)
889 0 : plm = (7.5_dp*x*x - 1.5_dp)*SQRT(1.0_dp - x*x)
890 : CASE (0)
891 198 : plm = 2.5_dp*x**3 - 1.5_dp*x
892 : END SELECT
893 : CASE (4)
894 1355233 : SELECT CASE (ABS(m))
895 : CASE DEFAULT
896 0 : CPABORT("l = 4 and m value out of bounds")
897 : CASE (4)
898 265026 : plm = 105.0_dp*(1.0_dp - x*x)**2
899 : CASE (3)
900 265026 : plm = 105.0_dp*x*(1.0_dp - x*x)**1.5_dp
901 : CASE (2)
902 265026 : plm = (52.5_dp*x*x - 7.5_dp)*(1.0_dp - x*x)
903 : CASE (1)
904 265026 : plm = (17.5_dp*x**3 - 7.5_dp*x)*SQRT(1.0_dp - x*x)
905 : CASE (0)
906 1189125 : plm = 4.375_dp*x**4 - 3.75_dp*x**2 + 0.375_dp
907 : END SELECT
908 : CASE (5)
909 815016 : SELECT CASE (ABS(m))
910 : CASE DEFAULT
911 0 : CPABORT("l = 5 and m value out of bounds")
912 : CASE (5)
913 30800 : plm = 945.0_dp*(1.0_dp - x*x)**2.5_dp
914 : CASE (4)
915 30800 : plm = 945.0_dp*x*(1.0_dp - x*x)**2
916 : CASE (3)
917 30800 : plm = -(-472.5_dp*x*x + 52.5_dp)*(1.0_dp - x*x)**1.5_dp
918 : CASE (2)
919 30800 : plm = (157.5_dp*x**3 - 52.5_dp*x)*(1.0_dp - x*x)
920 : CASE (1)
921 : plm = -(-39.375_dp*x**4 + 26.25_dp*x*x - &
922 30800 : 1.875_dp)*SQRT(1.0_dp - x*x)
923 : CASE (0)
924 166108 : plm = 7.875_dp*x**5 - 8.75_dp*x**3 + 1.875_dp*x
925 : END SELECT
926 : CASE (6)
927 1125612 : SELECT CASE (ABS(m))
928 : CASE DEFAULT
929 0 : CPABORT("l = 6 and m value out of bounds")
930 : CASE (6)
931 99832 : plm = 10395.0_dp*(1.0_dp - x*x)**3
932 : CASE (5)
933 99832 : plm = 10395.0_dp*x*(1.0_dp - x*x)**2.5_dp
934 : CASE (4)
935 99832 : plm = (5197.5_dp*x*x - 472.5_dp)*(1.0_dp - x*x)**2
936 : CASE (3)
937 : plm = -(-1732.5_dp*x**3 + 472.5_dp*x)* &
938 99832 : (1.0_dp - x*x)**1.5_dp
939 : CASE (2)
940 : plm = (433.125_dp*x**4 - 236.25_dp*x**2 + &
941 99832 : 13.125_dp)*(1.0_dp - x*x)
942 : CASE (1)
943 : plm = -(-86.625_dp*x**5 + 78.75_dp*x**3 - &
944 99832 : 13.125_dp*x)*SQRT(1.0_dp - x*x)
945 : CASE (0)
946 : plm = 14.4375_dp*x**6 - 19.6875_dp*x**4 + &
947 648908 : 6.5625_dp*x**2 - 0.3125_dp
948 : END SELECT
949 : CASE DEFAULT
950 476704 : mm = ABS(m)
951 476704 : IF (mm > l) CPABORT("m out of bounds")
952 : ! use recurence from numerical recipies
953 476704 : pmm = 1.0_dp
954 476704 : IF (mm > 0) THEN
955 452296 : somx2 = SQRT((1.0_dp - x)*(1.0_dp + x))
956 452296 : fact = 1.0_dp
957 2865512 : DO im = 1, mm
958 2413216 : pmm = pmm*fact*somx2
959 2865512 : fact = fact + 2.0_dp
960 : END DO
961 : END IF
962 2958217 : IF (l == mm) THEN
963 : plm = pmm
964 : ELSE
965 427888 : pmmp1 = x*REAL(2*mm + 1, KIND=dp)*pmm
966 427888 : IF (l == mm + 1) THEN
967 : plm = pmmp1
968 : ELSE
969 2138252 : DO il = mm + 2, l
970 : pll = (x*REAL(2*il - 1, KIND=dp)*pmmp1 - &
971 1759180 : REAL(il + mm - 1, KIND=dp)*pmm)/REAL(il - mm, KIND=dp)
972 1759180 : pmm = pmmp1
973 2138252 : pmmp1 = pll
974 : END DO
975 : plm = pll
976 : END IF
977 : END IF
978 : END SELECT
979 :
980 2481513 : END FUNCTION legendre
981 :
982 : ! **************************************************************************************************
983 : !> \brief ...
984 : !> \param x ...
985 : !> \param l ...
986 : !> \param m ...
987 : !> \return ...
988 : ! **************************************************************************************************
989 559516 : FUNCTION dlegendre(x, l, m) RESULT(dplm)
990 : REAL(KIND=dp), INTENT(IN) :: x
991 : INTEGER, INTENT(IN) :: l, m
992 : REAL(KIND=dp) :: dplm
993 :
994 : INTEGER :: mm
995 :
996 559516 : IF (ABS(x) > 1.0_dp) CPABORT("x value > 1")
997 559516 : SELECT CASE (l)
998 : CASE (0)
999 124 : dplm = 0.0_dp
1000 : CASE (1)
1001 248 : SELECT CASE (ABS(m))
1002 : CASE DEFAULT
1003 0 : CPABORT("l = 1 and m value out of bounds")
1004 : CASE (1)
1005 0 : dplm = -x/SQRT(1.0_dp - x*x)
1006 : CASE (0)
1007 124 : dplm = 1.0_dp
1008 : END SELECT
1009 : CASE (2)
1010 248 : SELECT CASE (ABS(m))
1011 : CASE DEFAULT
1012 0 : CPABORT("l = 2 and m value out of bounds")
1013 : CASE (2)
1014 0 : dplm = -6.0_dp*x
1015 : CASE (1)
1016 0 : dplm = 3.0_dp*SQRT(1.0_dp - x*x) - 3.0_dp*x*x/SQRT(1.0_dp - x*x)
1017 : CASE (0)
1018 124 : dplm = 3.0_dp*x
1019 : END SELECT
1020 : CASE (3)
1021 31552 : SELECT CASE (ABS(m))
1022 : CASE DEFAULT
1023 0 : CPABORT("l = 3 and m value out of bounds")
1024 : CASE (3)
1025 0 : dplm = -45.0_dp*SQRT(1.0_dp - x*x)*x
1026 : CASE (2)
1027 0 : dplm = 15.0_dp*(1.0_dp - x*x) - 30.0_dp*x*x
1028 : CASE (1)
1029 0 : dplm = 15.0_dp*x*SQRT(1.0_dp - x*x) - (x*(7.5_dp*x*x - 1.5_dp))/SQRT(1.0_dp - x*x)
1030 : CASE (0)
1031 124 : dplm = 7.5_dp*x*x - 1.5_dp
1032 : END SELECT
1033 : CASE (4)
1034 67640 : SELECT CASE (ABS(m))
1035 : CASE DEFAULT
1036 0 : CPABORT("l = 4 and m value out of bounds")
1037 : CASE (4)
1038 6984 : dplm = -420*x*(1 - x*x)
1039 : CASE (3)
1040 6984 : dplm = 105.0_dp*((1.0_dp - x*x)**1.5_dp - 3.0_dp*x*x*(1.0_dp - x*x)**0.5_dp)
1041 : CASE (2)
1042 6984 : dplm = 105.0_dp*x*(1.0_dp - x*x) - 2.0_dp*x*(52.5_dp*x*x - 7.5_dp)
1043 : CASE (1)
1044 6984 : IF (ABS(x) - 1.0_dp < EPSILON(1.0_dp)) THEN
1045 : dplm = 0.0_dp
1046 : ELSE
1047 : dplm = (17.5_dp*3.0_dp*x**2 - 7.5_dp)*SQRT(1.0_dp - x*x) - &
1048 0 : x*(17.5_dp*x**3 - 7.5_dp*x)/SQRT(1.0_dp - x*x)
1049 : END IF
1050 : CASE (0)
1051 31428 : dplm = 4.375_dp*4.0_dp*x**3 - 3.75_dp*2.0_dp*x
1052 : END SELECT
1053 : CASE (5)
1054 527716 : SELECT CASE (ABS(m))
1055 : CASE DEFAULT
1056 0 : CPABORT("l = 5 and m value out of bounds")
1057 : CASE (5)
1058 6584 : dplm = -945.0_dp*5.0_dp*x*(1.0_dp - x*x)**1.5_dp
1059 : CASE (4)
1060 6584 : dplm = 945.0_dp*((1.0_dp - x*x)**2 - 2.0_dp*x*x*(1.0_dp - x*x))
1061 : CASE (3)
1062 : dplm = 945.0_dp*x*(1.0_dp - x*x)**1.5_dp - &
1063 6584 : 3.0_dp*x*(472.5_dp*x*x - 52.5_dp)*(1.0_dp - x*x)**0.5_dp
1064 : CASE (2)
1065 : dplm = (3.0_dp*157.5_dp*x**2 - 52.5_dp)*(1.0_dp - x*x) - &
1066 6584 : (157.5_dp*x**3 - 52.5_dp*x)*(-2.0_dp*x)
1067 : CASE (1)
1068 6584 : IF (ABS(x) - 1.0_dp < EPSILON(1.0_dp)) THEN
1069 : dplm = 0.0_dp
1070 : ELSE
1071 : dplm = -(-39.375_dp*4.0_dp*x*x*x + 2.0_dp*26.25_dp*x)*SQRT(1.0_dp - x*x) + &
1072 0 : x*(-39.375_dp*x**4 + 26.25_dp*x*x - 1.875_dp)/SQRT(1.0_dp - x*x)
1073 : END IF
1074 : CASE (0)
1075 36212 : dplm = 5.0_dp*7.875_dp*x**4 - 3.0_dp*8.75_dp*x**2 + 1.875_dp
1076 : END SELECT
1077 : CASE (6)
1078 491504 : SELECT CASE (ABS(m))
1079 : CASE DEFAULT
1080 0 : CPABORT("l = 6 and m value out of bounds")
1081 : CASE (6)
1082 75616 : dplm = -10395.0_dp*6.0_dp*x*(1.0_dp - x*x)**2
1083 : CASE (5)
1084 75616 : dplm = 10395.0_dp*((1.0_dp - x*x)**2.5_dp - 5.0_dp*x*x*(1.0_dp - x*x)**1.5_dp)
1085 : CASE (4)
1086 : dplm = 2.0_dp*5197.5_dp*x*(1.0_dp - x*x)**2 - &
1087 75616 : 4.0_dp*x*(5197.5_dp*x*x - 472.5_dp)*(1.0_dp - x*x)
1088 : CASE (3)
1089 : dplm = -(-3.0_dp*1732.5_dp*x*x + 472.5_dp)*(1.0_dp - x*x)**1.5_dp + &
1090 75616 : (-1732.5_dp*x**3 + 472.5_dp*x)*3.0_dp*x*(1.0_dp - x*x)**0.5_dp
1091 : CASE (2)
1092 : dplm = (433.125_dp*4.0_dp*x**3 - 2.0_dp*236.25_dp*x)*(1.0_dp - x*x) - &
1093 75616 : 2.0_dp*x*(433.125_dp*x**4 - 236.25_dp*x**2 + 13.125_dp)
1094 : CASE (1)
1095 75616 : IF (ABS(x) - 1.0_dp < EPSILON(1.0_dp)) THEN
1096 : dplm = 0.0_dp
1097 : ELSE
1098 : dplm = -(-5.0_dp*86.625_dp*x**4 + 3.0_dp*78.75_dp**2 - 13.125_dp)*SQRT(1.0_dp - x*x) + &
1099 0 : x*(-86.625_dp*x**5 + 78.75_dp*x**3 - 13.125_dp*x)/SQRT(1.0_dp - x*x)
1100 : END IF
1101 : CASE (0)
1102 : dplm = 14.4375_dp*6.0_dp*x**5 - 19.6875_dp*4.0_dp*x**3 + &
1103 491504 : 6.5625_dp*2.0_dp*x
1104 : END SELECT
1105 : CASE DEFAULT
1106 0 : mm = ABS(m)
1107 0 : IF (mm > l) CPABORT("m out of bounds")
1108 :
1109 : !From Wikipedia: dPlm(x) = -(l+1)x/(x^2-1)*Plm(x) + (l-m+1)/(x^2-1)Pl+1m(x)
1110 559516 : IF (ABS(x) - 1.0_dp < EPSILON(1.0_dp)) THEN
1111 : dplm = 0.0_dp
1112 : ELSE
1113 : dplm = -REAL(l + 1, dp)*x/(x**2 - 1.0_dp)*legendre(x, l, m) &
1114 0 : + REAL(l - m + 1, dp)/(x**2 - 1.0_dp)*legendre(x, l + 1, m)
1115 : END IF
1116 : END SELECT
1117 :
1118 559516 : END FUNCTION dlegendre
1119 :
1120 : ! **************************************************************************************************
1121 : !> \brief ...
1122 : !> \param x ...
1123 : !> \param l ...
1124 : !> \return ...
1125 : ! **************************************************************************************************
1126 0 : FUNCTION dPof1(x, l)
1127 :
1128 : REAL(KIND=dp), INTENT(IN) :: x
1129 : INTEGER, INTENT(IN) :: l
1130 : REAL(KIND=dp) :: dPof1
1131 :
1132 0 : IF (ABS(x) - 1.0_dp > EPSILON(1.0_dp)) THEN
1133 0 : CPABORT("|x| is not 1")
1134 : END IF
1135 0 : IF (x > 0.0_dp) THEN
1136 0 : SELECT CASE (l)
1137 : CASE (:-1)
1138 0 : CPABORT("Negative l value")
1139 : CASE (0)
1140 0 : dPof1 = 0.0_dp
1141 : CASE (1)
1142 0 : dPof1 = 1.0_dp
1143 : CASE (2)
1144 0 : dPof1 = 3.0_dp
1145 : CASE (3)
1146 0 : dPof1 = 6.0_dp
1147 : CASE (4)
1148 0 : dPof1 = 10.0_dp
1149 : CASE (5)
1150 0 : dPof1 = 15.0_dp
1151 : CASE (6)
1152 0 : dPof1 = 21.0_dp
1153 : CASE (7:)
1154 0 : CPABORT("Not implemented")
1155 : END SELECT
1156 : ELSE
1157 0 : SELECT CASE (l)
1158 : CASE (:-1)
1159 0 : CPABORT("Negative l value")
1160 : CASE (0)
1161 0 : dPof1 = 0.0_dp
1162 : CASE (1)
1163 0 : dPof1 = 1.0_dp
1164 : CASE (2)
1165 0 : dPof1 = -3.0_dp
1166 : CASE (3)
1167 0 : dPof1 = 6.0_dp
1168 : CASE (4)
1169 0 : dPof1 = -10.0_dp
1170 : CASE (5)
1171 0 : dPof1 = 15.0_dp
1172 : CASE (6)
1173 0 : dPof1 = -21.0_dp
1174 : CASE (7:)
1175 0 : CPABORT("Not implemented")
1176 : END SELECT
1177 : END IF
1178 :
1179 0 : END FUNCTION dPof1
1180 :
1181 : ! **************************************************************************************************
1182 : !> \brief ...
1183 : !> \param n ...
1184 : !> \param k ...
1185 : !> \return ...
1186 : ! **************************************************************************************************
1187 3683456 : FUNCTION choose(n, k)
1188 :
1189 : INTEGER, INTENT(IN) :: n, k
1190 : REAL(KIND=dp) :: choose
1191 :
1192 3683456 : IF (n >= k) THEN
1193 3683456 : choose = REAL(NINT(fac(n)/(fac(k)*fac(n - k))), KIND=dp)
1194 : ELSE
1195 : choose = 0.0_dp
1196 : END IF
1197 :
1198 3683456 : END FUNCTION choose
1199 :
1200 : ! **************************************************************************************************
1201 : !> \brief ...
1202 : !> \param n ...
1203 : !> \param c ...
1204 : !> \param s ...
1205 : !> \return ...
1206 : ! **************************************************************************************************
1207 938072 : FUNCTION cosn(n, c, s)
1208 :
1209 : INTEGER, INTENT(IN) :: n
1210 : REAL(KIND=dp), INTENT(IN) :: c, s
1211 : REAL(KIND=dp) :: cosn
1212 :
1213 : INTEGER :: i, j
1214 :
1215 938072 : cosn = 0.0_dp
1216 938072 : IF (ABS(c) < EPSILON(1.0_dp) .OR. n == 0) THEN
1217 88178 : IF (MOD(n, 2) == 0) THEN
1218 43034 : cosn = (-1.0_dp)**INT(n/2)
1219 : ELSE
1220 : cosn = 0.0_dp
1221 : END IF
1222 : ELSE
1223 2977770 : DO i = n, 0, -2
1224 2977770 : IF (i == 0) THEN
1225 410912 : j = n
1226 : IF (j == 0) THEN
1227 : cosn = cosn + choose(n, j)
1228 410912 : ELSE IF (MOD(j/2, 2) == 0) THEN
1229 192947 : cosn = cosn + choose(n, j)*s**j
1230 : ELSE
1231 217965 : cosn = cosn - choose(n, j)*s**j
1232 : END IF
1233 : ELSE
1234 1716964 : j = n - i
1235 1716964 : IF (j == 0) THEN
1236 849894 : cosn = cosn + choose(n, j)*c**i
1237 867070 : ELSE IF (MOD(j/2, 2) == 0) THEN
1238 237166 : cosn = cosn + choose(n, j)*c**i*s**j
1239 : ELSE
1240 629904 : cosn = cosn - choose(n, j)*c**i*s**j
1241 : END IF
1242 : END IF
1243 : END DO
1244 : END IF
1245 :
1246 938072 : END FUNCTION cosn
1247 :
1248 : ! **************************************************************************************************
1249 : !> \brief ...
1250 : !> \param n ...
1251 : !> \param c ...
1252 : !> \param s ...
1253 : !> \return ...
1254 : ! **************************************************************************************************
1255 0 : FUNCTION dcosn_dcp(n, c, s)
1256 :
1257 : INTEGER, INTENT(IN) :: n
1258 : REAL(KIND=dp), INTENT(IN) :: c, s
1259 : REAL(KIND=dp) :: dcosn_dcp
1260 :
1261 : INTEGER :: i, j
1262 :
1263 0 : dcosn_dcp = 0.0_dp
1264 :
1265 0 : IF (s < EPSILON(1.0_dp)) THEN
1266 : dcosn_dcp = 0.0_dp
1267 : ELSE
1268 0 : DO i = n, 0, -2
1269 0 : IF (i == 0) THEN
1270 : dcosn_dcp = dcosn_dcp
1271 : ELSE
1272 0 : j = n - i
1273 0 : IF (MOD(j/2, 2) == 0) THEN
1274 0 : dcosn_dcp = dcosn_dcp + choose(n, j)*REAL(i, dp)*c**(i - 1)*s**j
1275 : ELSE
1276 0 : dcosn_dcp = dcosn_dcp - choose(n, j)*REAL(i, dp)*c**(i - 1)*s**j
1277 : END IF
1278 : END IF
1279 : END DO
1280 : END IF
1281 :
1282 0 : END FUNCTION dcosn_dcp
1283 :
1284 : ! **************************************************************************************************
1285 : !> \brief ...
1286 : !> \param n ...
1287 : !> \param c ...
1288 : !> \param s ...
1289 : !> \return ...
1290 : ! **************************************************************************************************
1291 0 : FUNCTION dcosn_dsp(n, c, s)
1292 :
1293 : INTEGER, INTENT(IN) :: n
1294 : REAL(KIND=dp), INTENT(IN) :: c, s
1295 : REAL(KIND=dp) :: dcosn_dsp
1296 :
1297 : INTEGER :: i, j
1298 :
1299 0 : dcosn_dsp = 0.0_dp
1300 0 : IF (c < EPSILON(1.0_dp) .OR. s < EPSILON(1.0_dp)) THEN
1301 : dcosn_dsp = 0.0_dp
1302 : ELSE
1303 0 : DO i = n, 0, -2
1304 0 : j = n - i
1305 0 : IF (j == 0) THEN
1306 : dcosn_dsp = dcosn_dsp
1307 : ELSE
1308 0 : IF (MOD(j/2, 2) == 0) THEN
1309 0 : dcosn_dsp = dcosn_dsp + choose(n, j)*REAL(j, dp)*s**(j - 1)*c**i
1310 : ELSE
1311 0 : dcosn_dsp = dcosn_dsp - choose(n, j)*REAL(j, dp)*s**(j - 1)*c**i
1312 : END IF
1313 : END IF
1314 : END DO
1315 : END IF
1316 :
1317 0 : END FUNCTION dcosn_dsp
1318 :
1319 : ! **************************************************************************************************
1320 : !> \brief ...
1321 : !> \param n ...
1322 : !> \param c ...
1323 : !> \param s ...
1324 : !> \return ...
1325 : ! **************************************************************************************************
1326 938072 : FUNCTION sinn(n, c, s)
1327 :
1328 : INTEGER, INTENT(IN) :: n
1329 : REAL(KIND=dp), INTENT(IN) :: c, s
1330 : REAL(KIND=dp) :: sinn
1331 :
1332 : INTEGER :: i, j
1333 :
1334 938072 : sinn = 0.0_dp
1335 :
1336 938072 : IF (ABS(s) < EPSILON(1.0_dp) .OR. n == 0) THEN
1337 : sinn = 0.0_dp
1338 848594 : ELSE IF (ABS(c) < EPSILON(1.0_dp)) THEN
1339 88178 : IF (MOD(n, 2) == 0) THEN
1340 : sinn = 0.0_dp
1341 : ELSE
1342 45144 : sinn = s*(-1.0_dp)**INT((n - 1)/2)
1343 : END IF
1344 : ELSE
1345 2315996 : DO i = n - 1, 0, -2
1346 2315996 : IF (i == 0) THEN
1347 393188 : j = n
1348 393188 : IF (MOD(j/2, 2) == 0) THEN
1349 214718 : sinn = sinn + choose(n, j)*s**j
1350 : ELSE
1351 178470 : sinn = sinn - choose(n, j)*s**j
1352 : END IF
1353 : ELSE
1354 1162392 : j = n - i
1355 1162392 : IF (MOD(j/2, 2) == 0) THEN
1356 768018 : sinn = sinn + choose(n, j)*c**i*s**j
1357 : ELSE
1358 394374 : sinn = sinn - choose(n, j)*c**i*s**j
1359 : END IF
1360 : END IF
1361 : END DO
1362 : END IF
1363 :
1364 938072 : END FUNCTION sinn
1365 :
1366 : ! **************************************************************************************************
1367 : !> \brief ...
1368 : !> \param n ...
1369 : !> \param c ...
1370 : !> \param s ...
1371 : !> \return ...
1372 : ! **************************************************************************************************
1373 0 : FUNCTION dsinn_dcp(n, c, s)
1374 :
1375 : INTEGER, INTENT(IN) :: n
1376 : REAL(KIND=dp), INTENT(IN) :: c, s
1377 : REAL(KIND=dp) :: dsinn_dcp
1378 :
1379 : INTEGER :: i, j
1380 :
1381 0 : dsinn_dcp = 0.0_dp
1382 :
1383 0 : IF (c < EPSILON(1.0_dp) .OR. s < EPSILON(1.0_dp)) THEN
1384 : dsinn_dcp = 0.0_dp
1385 : ELSE
1386 0 : DO i = n - 1, 0, -2
1387 0 : IF (i == 0) THEN
1388 : dsinn_dcp = dsinn_dcp
1389 : ELSE
1390 0 : j = n - i
1391 0 : IF (MOD(j/2, 2) == 0) THEN
1392 0 : dsinn_dcp = dsinn_dcp + choose(n, j)*REAL(i, dp)*c**(i - 1)*s**j
1393 : ELSE
1394 0 : dsinn_dcp = dsinn_dcp - choose(n, j)*REAL(i, dp)*c**(i - 1)*s**j
1395 : END IF
1396 : END IF
1397 : END DO
1398 : END IF
1399 :
1400 0 : END FUNCTION dsinn_dcp
1401 :
1402 : ! **************************************************************************************************
1403 : !> \brief ...
1404 : !> \param n ...
1405 : !> \param c ...
1406 : !> \param s ...
1407 : !> \return ...
1408 : ! **************************************************************************************************
1409 0 : FUNCTION dsinn_dsp(n, c, s)
1410 :
1411 : INTEGER, INTENT(IN) :: n
1412 : REAL(KIND=dp), INTENT(IN) :: c, s
1413 : REAL(KIND=dp) :: dsinn_dsp
1414 :
1415 : INTEGER :: i, j
1416 :
1417 0 : dsinn_dsp = 0.0_dp
1418 :
1419 0 : IF (c < EPSILON(1.0_dp) .OR. s < EPSILON(1.0_dp)) THEN
1420 : dsinn_dsp = 0.0_dp
1421 : ELSE
1422 0 : DO i = n - 1, 0, -2
1423 0 : j = n - i
1424 0 : IF (j == 0) THEN
1425 : dsinn_dsp = dsinn_dsp
1426 : ELSE
1427 0 : IF (MOD(j/2, 2) == 0) THEN
1428 0 : dsinn_dsp = dsinn_dsp + choose(n, j)*c**i*REAL(j, dp)*s**(j - 1)
1429 : ELSE
1430 0 : dsinn_dsp = dsinn_dsp - choose(n, j)*c**i*REAL(j, dp)*s**(j - 1)
1431 : END IF
1432 : END IF
1433 : END DO
1434 : END IF
1435 :
1436 0 : END FUNCTION dsinn_dsp
1437 :
1438 : END MODULE spherical_harmonics
|
