Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 : !> Clebsch-Gordon Coefficients
16 : !> Tested up to l=7 (i.e. L=14)
17 : !> \todo
18 : !> 1) Check if this definition is consistent with the
19 : !> Slater-Koster module
20 : !> \par History
21 : !> JGH 28-Feb-2002 : Change of sign convention (-1^m)
22 : !> JGH 1-Mar-2002 : Clebsch-Gordon Coefficients
23 : !> - Clebsch-Gordon coefficients and Wigner 3-j symbols added as generic routines using the
24 : !> standard normalization (19.09.2022, MK)
25 : !> \author JGH 6-Oct-2000, MK
26 : ! **************************************************************************************************
27 : MODULE spherical_harmonics
28 :
29 : USE kinds, ONLY: dp
30 : USE mathconstants, ONLY: fac,&
31 : maxfac,&
32 : pi
33 : #include "../base/base_uses.f90"
34 :
35 : IMPLICIT NONE
36 :
37 : PRIVATE
38 :
39 : PUBLIC :: y_lm, dy_lm
40 : PUBLIC :: clebsch_gordon_deallocate, clebsch_gordon_init, &
41 : clebsch_gordon
42 : PUBLIC :: legendre, dlegendre
43 : PUBLIC :: Clebsch_Gordon_coefficient, Wigner_3j
44 :
45 : INTERFACE y_lm
46 : MODULE PROCEDURE rvy_lm, rry_lm, cvy_lm, ccy_lm
47 : END INTERFACE
48 :
49 : INTERFACE dy_lm
50 : MODULE PROCEDURE dry_lm, dcy_lm
51 : END INTERFACE
52 :
53 : INTERFACE clebsch_gordon
54 : MODULE PROCEDURE clebsch_gordon_real, clebsch_gordon_complex
55 : END INTERFACE
56 :
57 : INTERFACE clebsch_gordon_getm
58 : MODULE PROCEDURE getm
59 : END INTERFACE
60 :
61 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'spherical_harmonics'
62 :
63 : REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: cg_table
64 : INTEGER :: lmax = -1
65 : REAL(KIND=dp) :: osq2, sq2
66 :
67 : CONTAINS
68 :
69 : ! Clebsch-Gordon Coefficients
70 :
71 : ! **************************************************************************************************
72 : !> \brief ...
73 : !> \param l1 ...
74 : !> \param m1 ...
75 : !> \param l2 ...
76 : !> \param m2 ...
77 : !> \param clm ...
78 : ! **************************************************************************************************
79 4096 : SUBROUTINE clebsch_gordon_complex(l1, m1, l2, m2, clm)
80 : INTEGER, INTENT(IN) :: l1, m1, l2, m2
81 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: clm
82 :
83 : INTEGER :: icase, ind, l, lm, lp, n
84 :
85 4096 : l = l1 + l2
86 4096 : IF (l > lmax) CALL clebsch_gordon_init(l)
87 4096 : n = l/2 + 1
88 4096 : IF (n > SIZE(clm)) CPABORT("Array too small")
89 :
90 4096 : IF ((m1 >= 0 .AND. m2 >= 0) .OR. (m1 < 0 .AND. m2 < 0)) THEN
91 : icase = 1
92 : ELSE
93 2016 : icase = 2
94 : END IF
95 4096 : ind = order(l1, m1, l2, m2)
96 :
97 4096 : DO lp = MOD(l, 2), l, 2
98 22800 : lm = lp/2 + 1
99 22800 : clm(lm) = cg_table(ind, lm, icase)
100 : END DO
101 :
102 4096 : END SUBROUTINE clebsch_gordon_complex
103 :
104 : ! **************************************************************************************************
105 : !> \brief ...
106 : !> \param l1 ...
107 : !> \param m1 ...
108 : !> \param l2 ...
109 : !> \param m2 ...
110 : !> \param rlm ...
111 : ! **************************************************************************************************
112 113266 : SUBROUTINE clebsch_gordon_real(l1, m1, l2, m2, rlm)
113 : INTEGER, INTENT(IN) :: l1, m1, l2, m2
114 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: rlm
115 :
116 : INTEGER :: icase1, icase2, ind, l, lm, lp, mm(2), n
117 : REAL(KIND=dp) :: xsi
118 :
119 113266 : l = l1 + l2
120 113266 : IF (l > lmax) CALL clebsch_gordon_init(l)
121 113266 : n = l/2 + 1
122 113266 : IF (n > SIZE(rlm, 1)) CPABORT("Array too small")
123 :
124 113266 : ind = order(l1, m1, l2, m2)
125 113266 : mm = getm(m1, m2)
126 113266 : IF ((m1 >= 0 .AND. m2 >= 0) .OR. (m1 < 0 .AND. m2 < 0)) THEN
127 : icase1 = 1
128 : icase2 = 2
129 : ELSE
130 51512 : icase1 = 2
131 51512 : icase2 = 1
132 : END IF
133 :
134 113266 : DO lp = MOD(l, 2), l, 2
135 315274 : lm = lp/2 + 1
136 315274 : xsi = get_factor(m1, m2, mm(1))
137 315274 : rlm(lm, 1) = xsi*cg_table(ind, lm, icase1)
138 315274 : xsi = get_factor(m1, m2, mm(2))
139 315274 : rlm(lm, 2) = xsi*cg_table(ind, lm, icase2)
140 : END DO
141 :
142 113266 : END SUBROUTINE clebsch_gordon_real
143 :
144 : ! **************************************************************************************************
145 : !> \brief ...
146 : !> \param m1 ...
147 : !> \param m2 ...
148 : !> \return ...
149 : ! **************************************************************************************************
150 113266 : FUNCTION getm(m1, m2) RESULT(m)
151 : INTEGER, INTENT(IN) :: m1, m2
152 : INTEGER, DIMENSION(2) :: m
153 :
154 : INTEGER :: mm, mp
155 :
156 113266 : mp = m1 + m2
157 113266 : mm = m1 - m2
158 113266 : IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
159 51512 : mp = -ABS(mp)
160 51512 : mm = -ABS(mm)
161 : ELSE
162 61754 : mp = ABS(mp)
163 61754 : mm = ABS(mm)
164 : END IF
165 113266 : m(1) = mp
166 113266 : m(2) = mm
167 113266 : END FUNCTION getm
168 :
169 : ! **************************************************************************************************
170 : !> \brief ...
171 : !> \param m1 ...
172 : !> \param m2 ...
173 : !> \param m ...
174 : !> \return ...
175 : ! **************************************************************************************************
176 630548 : FUNCTION get_factor(m1, m2, m) RESULT(f)
177 : INTEGER, INTENT(IN) :: m1, m2, m
178 : REAL(KIND=dp) :: f
179 :
180 : INTEGER :: mx, my
181 :
182 630548 : f = 1.0_dp
183 630548 : IF (ABS(m1) >= ABS(m2)) THEN
184 : mx = m1
185 : my = m2
186 : ELSE
187 237632 : mx = m2
188 237632 : my = m1
189 : END IF
190 630548 : IF (mx*my == 0) THEN
191 : f = 1.0_dp
192 370816 : ELSE IF (m == 0) THEN
193 70128 : IF (ABS(mx) /= ABS(my)) WRITE (*, '(A,3I6)') " 1) Illegal Case ", m1, m2, m
194 70128 : IF (mx*my > 0) THEN
195 : f = 1.0_dp
196 : ELSE
197 35064 : f = 0.0_dp
198 : END IF
199 300688 : ELSE IF (ABS(mx) + ABS(my) == m) THEN
200 92704 : f = osq2
201 92704 : IF (mx < 0) f = -osq2
202 207984 : ELSE IF (ABS(mx) + ABS(my) == -m) THEN
203 92704 : f = osq2
204 115280 : ELSE IF (ABS(mx) - ABS(my) == -m) THEN
205 57640 : IF (mx*my > 0) WRITE (*, '(A,3I6)') " 2) Illegal Case ", m1, m2, m
206 57640 : IF (mx > 0) f = -osq2
207 57640 : IF (mx < 0) f = osq2
208 57640 : ELSE IF (ABS(mx) - ABS(my) == m) THEN
209 57640 : IF (mx*my < 0) WRITE (*, '(A,3I6)') " 3) Illegal Case ", m1, m2, m
210 57640 : f = osq2
211 : ELSE
212 0 : WRITE (*, '(A,3I6)') " 4) Illegal Case ", m1, m2, m
213 : END IF
214 630548 : END FUNCTION get_factor
215 :
216 : ! **************************************************************************************************
217 : !> \brief ...
218 : ! **************************************************************************************************
219 1237 : SUBROUTINE clebsch_gordon_deallocate()
220 : CHARACTER(len=*), PARAMETER :: routineN = 'clebsch_gordon_deallocate'
221 :
222 : INTEGER :: handle
223 :
224 1237 : CALL timeset(routineN, handle)
225 :
226 1237 : IF (ALLOCATED(cg_table)) THEN
227 1237 : DEALLOCATE (cg_table)
228 : END IF
229 1237 : CALL timestop(handle)
230 :
231 1237 : END SUBROUTINE clebsch_gordon_deallocate
232 :
233 : ! **************************************************************************************************
234 : !> \brief ...
235 : !> \param l ...
236 : ! **************************************************************************************************
237 1244 : SUBROUTINE clebsch_gordon_init(l)
238 : INTEGER, INTENT(IN) :: l
239 :
240 : CHARACTER(len=*), PARAMETER :: routineN = 'clebsch_gordon_init'
241 :
242 : INTEGER :: handle, i1, i2, ix, iy, l1, l2, lp, m, &
243 : m1, m2, ml, mp, n
244 :
245 1244 : CALL timeset(routineN, handle)
246 :
247 1244 : sq2 = SQRT(2.0_dp)
248 1244 : osq2 = 1.0_dp/sq2
249 :
250 1244 : IF (l < 0) CPABORT("l < 0")
251 1244 : IF (ALLOCATED(cg_table)) THEN
252 7 : DEALLOCATE (cg_table)
253 : END IF
254 : ! maximum size of table
255 1244 : n = (l**4 + 6*l**3 + 15*l**2 + 18*l + 8)/8
256 1244 : m = l + 1
257 6220 : ALLOCATE (cg_table(n, m, 2))
258 1244 : lmax = l
259 :
260 7030 : DO l1 = 0, lmax
261 24806 : DO m1 = 0, l1
262 17776 : iy = (l1*(l1 + 1))/2 + m1 + 1
263 68852 : DO l2 = l1, lmax
264 45290 : ml = 0
265 45290 : IF (l1 == l2) ml = m1
266 258904 : DO m2 = ml, l2
267 195838 : ix = (l2*(l2 + 1))/2 + m2 + 1
268 195838 : i1 = (ix*(ix - 1))/2 + iy
269 1221734 : DO lp = MOD(l1 + l2, 2), l1 + l2, 2
270 980606 : i2 = lp/2 + 1
271 980606 : mp = m2 + m1
272 980606 : cg_table(i1, i2, 1) = cgc(l1, m1, l2, m2, lp, mp)
273 980606 : mp = ABS(m2 - m1)
274 1176444 : IF (m2 >= m1) THEN
275 793190 : cg_table(i1, i2, 2) = cgc(l1, m1, lp, mp, l2, m2)
276 : ELSE
277 187416 : cg_table(i1, i2, 2) = cgc(l2, m2, lp, mp, l1, m1)
278 : END IF
279 : END DO
280 : END DO
281 : END DO
282 : END DO
283 : END DO
284 :
285 1244 : CALL timestop(handle)
286 :
287 1244 : END SUBROUTINE clebsch_gordon_init
288 :
289 : ! **************************************************************************************************
290 : !> \brief ...
291 : !> \param l1 ...
292 : !> \param m1 ...
293 : !> \param l2 ...
294 : !> \param m2 ...
295 : !> \param lp ...
296 : !> \param mp ...
297 : !> \return ...
298 : ! **************************************************************************************************
299 1961212 : FUNCTION cgc(l1, m1, l2, m2, lp, mp)
300 : INTEGER, INTENT(IN) :: l1, m1, l2, m2, lp, mp
301 : REAL(KIND=dp) :: cgc
302 :
303 : INTEGER :: la, lb, ll, ma, mb, s, t, tmax, tmin, &
304 : z1, z2, z3, z4, z5
305 : REAL(KIND=dp) :: f1, f2, pref
306 :
307 : ! m1 >= 0; m2 >= 0; mp >= 0
308 :
309 1961212 : IF (m1 < 0 .OR. m2 < 0 .OR. mp < 0) THEN
310 0 : WRITE (*, *) l1, l2, lp
311 0 : WRITE (*, *) m1, m2, mp
312 0 : CPABORT("Illegal input values")
313 : END IF
314 1961212 : IF (l2 < l1) THEN
315 346922 : la = l2
316 346922 : ma = m2
317 346922 : lb = l1
318 346922 : mb = m1
319 : ELSE
320 1614290 : la = l1
321 1614290 : ma = m1
322 1614290 : lb = l2
323 1614290 : mb = m2
324 : END IF
325 :
326 : IF (MOD(la + lb + lp, 2) == 0 .AND. la + lb >= lp .AND. lp >= lb - la &
327 1961212 : .AND. lb - mb >= 0) THEN
328 1657627 : ll = (2*lp + 1)*(2*la + 1)*(2*lb + 1)
329 : pref = 1.0_dp/SQRT(4.0_dp*pi)*0.5_dp*SQRT(REAL(ll, dp)* &
330 : (sfac(lp - mp)/sfac(lp + mp))* &
331 1657627 : (sfac(la - ma)/sfac(la + ma))*(sfac(lb - mb)/sfac(lb + mb)))
332 1657627 : s = (la + lb + lp)/2
333 1657627 : tmin = MAX(0, -lb + la - mp)
334 1657627 : tmax = MIN(lb + la - mp, lp - mp, la - ma)
335 : f1 = REAL(2*(-1)**(s - lb - ma), KIND=dp)*(sfac(lb + mb)/sfac(lb - mb))* &
336 : sfac(la + ma)/(sfac(s - lp)*sfac(s - lb))*sfac(2*s - 2*la)/sfac(s - la)* &
337 1657627 : (sfac(s)/sfac(2*s + 1))
338 1657627 : f2 = 0.0_dp
339 5116631 : DO t = tmin, tmax
340 3459004 : z1 = lp + mp + t
341 3459004 : z2 = la + lb - mp - t
342 3459004 : z3 = lp - mp - t
343 3459004 : z4 = lb - la + mp + t
344 3459004 : z5 = la - ma - t
345 5116631 : f2 = f2 + (-1)**t*(sfac(z1)/(sfac(t)*sfac(z3)))*(sfac(z2)/(sfac(z4)*sfac(z5)))
346 : END DO
347 1657627 : cgc = pref*f1*f2
348 : ELSE
349 : cgc = 0.0_dp
350 : END IF
351 :
352 1961212 : END FUNCTION cgc
353 :
354 : ! **************************************************************************************************
355 : !> \brief Calculate factorial even for integer values larger than the tabulated (pre-computed)
356 : !> values of up to 30!
357 : !> \param n Integer number for which the factorial has to be returned
358 : !> \return Factorial n!
359 : ! **************************************************************************************************
360 45618429 : FUNCTION sfac(n)
361 : INTEGER, INTENT(IN) :: n
362 : REAL(KIND=dp) :: sfac
363 :
364 : INTEGER :: i
365 :
366 45618429 : IF (n > 170) THEN
367 0 : CPABORT("Factorials greater than 170! cannot be computed with double-precision")
368 45618429 : ELSE IF (n > maxfac) THEN
369 : sfac = fac(maxfac)
370 2832511 : DO i = maxfac + 1, n
371 2832511 : sfac = REAL(i, KIND=dp)*sfac
372 : END DO
373 45204494 : ELSE IF (n >= 0) THEN
374 44822763 : sfac = fac(n)
375 : ELSE
376 : sfac = 1.0_dp
377 : END IF
378 :
379 45618429 : END FUNCTION sfac
380 :
381 : ! **************************************************************************************************
382 : !> \brief ...
383 : !> \param l1 ...
384 : !> \param m1 ...
385 : !> \param l2 ...
386 : !> \param m2 ...
387 : !> \return ...
388 : ! **************************************************************************************************
389 117362 : FUNCTION order(l1, m1, l2, m2) RESULT(ind)
390 : INTEGER, INTENT(IN) :: l1, m1, l2, m2
391 : INTEGER :: ind
392 :
393 : INTEGER :: i1, i2, ix, iy
394 :
395 117362 : i1 = (l1*(l1 + 1))/2 + ABS(m1) + 1
396 117362 : i2 = (l2*(l2 + 1))/2 + ABS(m2) + 1
397 117362 : ix = MAX(i1, i2)
398 117362 : iy = MIN(i1, i2)
399 117362 : ind = (ix*(ix - 1))/2 + iy
400 117362 : END FUNCTION order
401 :
402 : ! Calculation of Spherical Harmonics
403 :
404 : ! **************************************************************************************************
405 : !> \brief ...
406 : !> \param r ...
407 : !> \param y ...
408 : !> \param l ...
409 : !> \param m ...
410 : ! **************************************************************************************************
411 79734 : SUBROUTINE rvy_lm(r, y, l, m)
412 : !
413 : ! Real Spherical Harmonics
414 : ! _ _
415 : ! | [(2l+1)(l-|m|)!] |1/2 m cos(m p) m>=0
416 : ! Y_lm ( t, p ) = |---------------------| P_l(cos(t))
417 : ! |[2Pi(1+d_m0)(l+|m|)!]| sin(|m| p) m<0
418 : !
419 : ! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1
420 : !
421 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r
422 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: y
423 : INTEGER, INTENT(IN) :: l, m
424 :
425 : INTEGER :: i
426 : REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, sp, t, z
427 :
428 79734 : SELECT CASE (l)
429 : CASE (:-1)
430 0 : CPABORT("Negative l value")
431 : CASE (0)
432 2567 : pf = SQRT(1.0_dp/(4.0_dp*pi))
433 2567 : IF (m /= 0) CPABORT("l = 0 and m value out of bounds")
434 346689 : y(:) = pf
435 : CASE (1)
436 19675 : SELECT CASE (m)
437 : CASE DEFAULT
438 0 : CPABORT("l = 1 and m value out of bounds")
439 : CASE (1)
440 2453 : pf = SQRT(3.0_dp/(4.0_dp*pi))
441 319347 : y(:) = pf*r(1, :)
442 : CASE (0)
443 2477 : pf = SQRT(3.0_dp/(4.0_dp*pi))
444 333531 : y(:) = pf*r(3, :)
445 : CASE (-1)
446 2453 : pf = SQRT(3.0_dp/(4.0_dp*pi))
447 324277 : y(:) = pf*r(2, :)
448 : END SELECT
449 : CASE (2)
450 25939 : SELECT CASE (m)
451 : CASE DEFAULT
452 0 : CPABORT("l = 2 and m value out of bounds")
453 : CASE (2)
454 2447 : pf = SQRT(15.0_dp/(16.0_dp*pi))
455 315801 : y(:) = pf*(r(1, :)*r(1, :) - r(2, :)*r(2, :))
456 : CASE (1)
457 2453 : pf = SQRT(15.0_dp/(4.0_dp*pi))
458 319347 : y(:) = pf*r(3, :)*r(1, :)
459 : CASE (0)
460 2492 : pf = SQRT(5.0_dp/(16.0_dp*pi))
461 342396 : y(:) = pf*(3.0_dp*r(3, :)*r(3, :) - 1.0_dp)
462 : CASE (-1)
463 2453 : pf = SQRT(15.0_dp/(4.0_dp*pi))
464 319347 : y(:) = pf*r(3, :)*r(2, :)
465 : CASE (-2)
466 2447 : pf = SQRT(15.0_dp/(16.0_dp*pi))
467 325646 : y(:) = pf*2.0_dp*r(1, :)*r(2, :)
468 : END SELECT
469 : CASE (3)
470 57492 : SELECT CASE (m)
471 : CASE DEFAULT
472 0 : CPABORT("l = 3 and m value out of bounds")
473 : CASE (3)
474 1927 : pf = SQRT(35.0_dp/(32.0_dp*pi))
475 273225 : y(:) = pf*r(1, :)*(r(1, :)**2 - 3.0_dp*r(2, :)**2)
476 : CASE (2)
477 1943 : pf = SQRT(105.0_dp/(16.0_dp*pi))
478 282681 : y(:) = pf*r(3, :)*(r(1, :)**2 - r(2, :)**2)
479 : CASE (1)
480 1961 : pf = SQRT(21.0_dp/(32.0_dp*pi))
481 293319 : y(:) = pf*r(1, :)*(5.0_dp*r(3, :)*r(3, :) - 1.0_dp)
482 : CASE (0)
483 1985 : pf = SQRT(7.0_dp/(16.0_dp*pi))
484 307503 : y(:) = pf*r(3, :)*(5.0_dp*r(3, :)*r(3, :) - 3.0_dp)
485 : CASE (-1)
486 1961 : pf = SQRT(21.0_dp/(32.0_dp*pi))
487 293319 : y(:) = pf*r(2, :)*(5.0_dp*r(3, :)*r(3, :) - 1.0_dp)
488 : CASE (-2)
489 1943 : pf = SQRT(105.0_dp/(16.0_dp*pi))
490 282681 : y(:) = pf*2.0_dp*r(1, :)*r(2, :)*r(3, :)
491 : CASE (-3)
492 1927 : pf = SQRT(35.0_dp/(32.0_dp*pi))
493 284945 : y(:) = pf*r(2, :)*(3.0_dp*r(1, :)**2 - r(2, :)**2)
494 : END SELECT
495 : CASE DEFAULT
496 43845 : IF (m < -l .OR. m > l) CPABORT("m value out of bounds")
497 43845 : lpm = fac(l + ABS(m))
498 43845 : lmm = fac(l - ABS(m))
499 43845 : IF (m == 0) THEN
500 : t = 4.0_dp*pi
501 : ELSE
502 39624 : t = 2.0_dp*pi
503 : END IF
504 43845 : IF (ABS(lpm) < EPSILON(1.0_dp)) THEN
505 0 : pf = REAL(2*l + 1, KIND=dp)/t
506 : ELSE
507 43845 : pf = (REAL(2*l + 1, KIND=dp)*lmm)/(t*lpm)
508 : END IF
509 43845 : pf = SQRT(pf)
510 16540833 : DO i = 1, SIZE(r, 2)
511 16417254 : z = r(3, i)
512 : ! plm = legendre ( z, l, m )
513 16417254 : plm = legendre(z, l, ABS(m))
514 16461099 : IF (m == 0) THEN
515 1502358 : y(i) = pf*plm
516 : ELSE
517 14914896 : rxy = SQRT(r(1, i)**2 + r(2, i)**2)
518 14914896 : IF (rxy < EPSILON(1.0_dp)) THEN
519 79248 : y(i) = 0.0_dp
520 : ELSE
521 14835648 : cp = r(1, i)/rxy
522 14835648 : sp = r(2, i)/rxy
523 14835648 : IF (m > 0) THEN
524 7417824 : y(i) = pf*plm*cosn(m, cp, sp)
525 : ELSE
526 7417824 : y(i) = pf*plm*sinn(-m, cp, sp)
527 : END IF
528 : END IF
529 : END IF
530 : END DO
531 : END SELECT
532 :
533 79734 : END SUBROUTINE rvy_lm
534 :
535 : ! **************************************************************************************************
536 : !> \brief ...
537 : !> \param r ...
538 : !> \param y ...
539 : !> \param l ...
540 : !> \param m ...
541 : ! **************************************************************************************************
542 530080 : SUBROUTINE rry_lm(r, y, l, m)
543 : !
544 : ! Real Spherical Harmonics
545 : ! _ _
546 : ! | [(2l+1)(l-|m|)!] |1/2 m cos(m p) m>=0
547 : ! Y_lm ( t, p ) = |---------------------| P_l(cos(t))
548 : ! |[2Pi(1+d_m0)(l+|m|)!]| sin(|m| p) m<0
549 : !
550 : ! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1
551 : !
552 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r
553 : REAL(KIND=dp), INTENT(OUT) :: y
554 : INTEGER, INTENT(IN) :: l, m
555 :
556 : REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, sp, t, z
557 :
558 530080 : SELECT CASE (l)
559 : CASE (:-1)
560 0 : CPABORT("Negative l value")
561 : CASE (0)
562 0 : pf = SQRT(1.0_dp/(4.0_dp*pi))
563 0 : IF (m /= 0) CPABORT("l = 0 and m value out of bounds")
564 0 : y = pf
565 : CASE (1)
566 431184 : SELECT CASE (m)
567 : CASE DEFAULT
568 0 : CPABORT("l = 1 and m value out of bounds")
569 : CASE (1)
570 98416 : pf = SQRT(3.0_dp/(4.0_dp*pi))
571 98416 : y = pf*r(1)
572 : CASE (0)
573 0 : pf = SQRT(3.0_dp/(4.0_dp*pi))
574 0 : y = pf*r(3)
575 : CASE (-1)
576 98416 : pf = SQRT(3.0_dp/(4.0_dp*pi))
577 196832 : y = pf*r(2)
578 : END SELECT
579 : CASE (2)
580 283264 : SELECT CASE (m)
581 : CASE DEFAULT
582 0 : CPABORT("l = 2 and m value out of bounds")
583 : CASE (2)
584 58588 : pf = SQRT(15.0_dp/(16.0_dp*pi))
585 58588 : y = pf*(r(1)*r(1) - r(2)*r(2))
586 : CASE (1)
587 58588 : pf = SQRT(15.0_dp/(4.0_dp*pi))
588 58588 : y = pf*r(3)*r(1)
589 : CASE (0)
590 0 : pf = SQRT(5.0_dp/(16.0_dp*pi))
591 0 : y = pf*(3.0_dp*r(3)*r(3) - 1.0_dp)
592 : CASE (-1)
593 58588 : pf = SQRT(15.0_dp/(4.0_dp*pi))
594 58588 : y = pf*r(3)*r(2)
595 : CASE (-2)
596 58588 : pf = SQRT(15.0_dp/(16.0_dp*pi))
597 234352 : y = pf*2.0_dp*r(1)*r(2)
598 : END SELECT
599 : CASE (3)
600 98896 : SELECT CASE (m)
601 : CASE DEFAULT
602 0 : CPABORT("l = 3 and m value out of bounds")
603 : CASE (3)
604 8152 : pf = SQRT(35.0_dp/(32.0_dp*pi))
605 8152 : y = pf*r(1)*(r(1)**2 - 3.0_dp*r(2)**2)
606 : CASE (2)
607 8152 : pf = SQRT(105.0_dp/(16.0_dp*pi))
608 8152 : y = pf*r(3)*(r(1)**2 - r(2)**2)
609 : CASE (1)
610 8152 : pf = SQRT(21.0_dp/(32.0_dp*pi))
611 8152 : y = pf*r(1)*(5.0_dp*r(3)*r(3) - 1.0_dp)
612 : CASE (0)
613 0 : pf = SQRT(7.0_dp/(16.0_dp*pi))
614 0 : y = pf*r(3)*(5.0_dp*r(3)*r(3) - 3.0_dp)
615 : CASE (-1)
616 8152 : pf = SQRT(21.0_dp/(32.0_dp*pi))
617 8152 : y = pf*r(2)*(5.0_dp*r(3)*r(3) - 1.0_dp)
618 : CASE (-2)
619 8152 : pf = SQRT(105.0_dp/(16.0_dp*pi))
620 8152 : y = pf*2.0_dp*r(1)*r(2)*r(3)
621 : CASE (-3)
622 8152 : pf = SQRT(35.0_dp/(32.0_dp*pi))
623 48912 : y = pf*r(2)*(3.0_dp*r(1)**2 - r(2)**2)
624 : END SELECT
625 : CASE DEFAULT
626 49984 : IF (m < -l .OR. m > l) CPABORT("m value out of bounds")
627 49984 : lpm = fac(l + ABS(m))
628 49984 : lmm = fac(l - ABS(m))
629 49984 : IF (m == 0) THEN
630 : t = 4.0_dp*pi
631 : ELSE
632 49984 : t = 2.0_dp*pi
633 : END IF
634 49984 : IF (ABS(lpm) < EPSILON(1.0_dp)) THEN
635 0 : pf = REAL(2*l + 1, KIND=dp)/t
636 : ELSE
637 49984 : pf = (REAL(2*l + 1, KIND=dp)*lmm)/(t*lpm)
638 : END IF
639 49984 : pf = SQRT(pf)
640 49984 : z = r(3)
641 49984 : plm = legendre(z, l, m)
642 580064 : IF (m == 0) THEN
643 0 : y = pf*plm
644 : ELSE
645 49984 : rxy = SQRT(r(1)**2 + r(2)**2)
646 49984 : IF (rxy < EPSILON(1.0_dp)) THEN
647 496 : y = 0.0_dp
648 : ELSE
649 49488 : cp = r(1)/rxy
650 49488 : sp = r(2)/rxy
651 49488 : IF (m > 0) THEN
652 24744 : y = pf*plm*cosn(m, cp, sp)
653 : ELSE
654 24744 : y = pf*plm*sinn(-m, cp, sp)
655 : END IF
656 : END IF
657 : END IF
658 : END SELECT
659 :
660 530080 : END SUBROUTINE rry_lm
661 :
662 : ! **************************************************************************************************
663 : !> \brief ...
664 : !> \param r ...
665 : !> \param y ...
666 : !> \param l ...
667 : !> \param m ...
668 : ! **************************************************************************************************
669 19400 : SUBROUTINE cvy_lm(r, y, l, m)
670 : !
671 : ! Complex Spherical Harmonics
672 : ! _ _
673 : ! | [(2l+1)(l-|m|)!] |1/2 m
674 : ! Y_lm ( t, p ) = |------------------| P_l(cos(t)) Exp[ i m p ]
675 : ! | [4Pi(l+|m|)!]| |
676 : !
677 : ! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1
678 : !
679 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r
680 : COMPLEX(KIND=dp), DIMENSION(:), INTENT(OUT) :: y
681 : INTEGER, INTENT(IN) :: l, m
682 :
683 : INTEGER :: i, n
684 : REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, sp, t, ym, &
685 : yp, z
686 :
687 19400 : n = SIZE(r, 2)
688 19400 : SELECT CASE (l)
689 : CASE (:-1)
690 0 : CPABORT("Negative l value")
691 : CASE (0)
692 241 : pf = SQRT(1.0_dp/(4.0_dp*pi))
693 241 : IF (m /= 0) CPABORT("l = 0 and m value out of bounds")
694 142431 : y(:) = pf
695 : CASE (1)
696 1857 : SELECT CASE (m)
697 : CASE DEFAULT
698 0 : CPABORT("l = 1 and m value out of bounds")
699 : CASE (1)
700 137703 : pf = SQRT(3.0_dp/(8.0_dp*pi))
701 137703 : DO i = 1, n
702 137470 : yp = r(1, i)
703 137470 : ym = r(2, i)
704 137703 : y(i) = pf*CMPLX(yp, ym, KIND=dp)
705 : END DO
706 : CASE (0)
707 233 : pf = SQRT(3.0_dp/(4.0_dp*pi))
708 137703 : y(:) = pf*r(3, :)
709 : CASE (-1)
710 137703 : pf = SQRT(3.0_dp/(8.0_dp*pi))
711 138169 : DO i = 1, n
712 137470 : yp = r(1, i)
713 137470 : ym = r(2, i)
714 137703 : y(i) = pf*CMPLX(yp, -ym, KIND=dp)
715 : END DO
716 : END SELECT
717 : CASE (2)
718 2703 : SELECT CASE (m)
719 : CASE DEFAULT
720 0 : CPABORT("l = 2 and m value out of bounds")
721 : CASE (2)
722 133566 : pf = SQRT(15.0_dp/(32.0_dp*pi))
723 133566 : DO i = 1, n
724 133340 : yp = (r(1, i)*r(1, i) - r(2, i)*r(2, i))
725 133340 : ym = 2.0_dp*r(1, i)*r(2, i)
726 133566 : y(i) = pf*CMPLX(yp, ym, KIND=dp)
727 : END DO
728 : CASE (1)
729 137703 : pf = SQRT(15.0_dp/(8.0_dp*pi))
730 137703 : DO i = 1, n
731 137470 : yp = r(3, i)*r(1, i)
732 137470 : ym = r(3, i)*r(2, i)
733 137703 : y(i) = pf*CMPLX(yp, ym, KIND=dp)
734 : END DO
735 : CASE (0)
736 240 : pf = SQRT(5.0_dp/(16.0_dp*pi))
737 141840 : y(:) = pf*(3.0_dp*r(3, :)*r(3, :) - 1.0_dp)
738 : CASE (-1)
739 137703 : pf = SQRT(15.0_dp/(8.0_dp*pi))
740 137703 : DO i = 1, n
741 137470 : yp = r(3, i)*r(1, i)
742 137470 : ym = r(3, i)*r(2, i)
743 137703 : y(i) = pf*CMPLX(yp, -ym, KIND=dp)
744 : END DO
745 : CASE (-2)
746 133566 : pf = SQRT(15.0_dp/(32.0_dp*pi))
747 134498 : DO i = 1, n
748 133340 : yp = (r(1, i)*r(1, i) - r(2, i)*r(2, i))
749 133340 : ym = 2.0_dp*r(1, i)*r(2, i)
750 133566 : y(i) = pf*CMPLX(yp, -ym, KIND=dp)
751 : END DO
752 : END SELECT
753 : CASE (3)
754 17302 : SELECT CASE (m)
755 : CASE DEFAULT
756 0 : CPABORT("l = 3 and m value out of bounds")
757 : CASE (3)
758 122337 : pf = SQRT(35.0_dp/(64.0_dp*pi))
759 122337 : DO i = 1, n
760 122130 : yp = r(1, i)*(r(1, i)**2 - 3.0_dp*r(2, i)**2)
761 122130 : ym = r(2, i)*(3.0_dp*r(1, i)**2 - r(2, i)**2)
762 122337 : y(i) = pf*CMPLX(yp, ym, KIND=dp)
763 : END DO
764 : CASE (2)
765 129429 : pf = SQRT(105.0_dp/(32.0_dp*pi))
766 129429 : DO i = 1, n
767 129210 : yp = r(3, i)*(r(1, i)**2 - r(2, i)**2)
768 129210 : ym = 2.0_dp*r(1, i)*r(2, i)*r(3, i)
769 129429 : y(i) = pf*CMPLX(yp, ym, KIND=dp)
770 : END DO
771 : CASE (1)
772 136521 : pf = SQRT(21.0_dp/(64.0_dp*pi))
773 136521 : DO i = 1, n
774 136290 : yp = r(1, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
775 136290 : ym = r(2, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
776 136521 : y(i) = pf*CMPLX(yp, ym, KIND=dp)
777 : END DO
778 : CASE (0)
779 231 : pf = SQRT(7.0_dp/(16.0_dp*pi))
780 136521 : y(:) = pf*r(3, :)*(5.0_dp*r(3, :)*r(3, :) - 3.0_dp)
781 : CASE (-1)
782 136521 : pf = SQRT(21.0_dp/(64.0_dp*pi))
783 136521 : DO i = 1, n
784 136290 : yp = r(1, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
785 136290 : ym = r(2, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
786 136521 : y(i) = pf*CMPLX(yp, -ym, KIND=dp)
787 : END DO
788 : CASE (-2)
789 129429 : pf = SQRT(105.0_dp/(32.0_dp*pi))
790 129429 : DO i = 1, n
791 129210 : yp = r(3, i)*(r(1, i)**2 - r(2, i)**2)
792 129210 : ym = 2.0_dp*r(1, i)*r(2, i)*r(3, i)
793 129429 : y(i) = pf*CMPLX(yp, -ym, KIND=dp)
794 : END DO
795 : CASE (-3)
796 122337 : pf = SQRT(35.0_dp/(64.0_dp*pi))
797 123675 : DO i = 1, n
798 122130 : yp = r(1, i)*(r(1, i)**2 - 3.0_dp*r(2, i)**2)
799 122130 : ym = r(2, i)*(3.0_dp*r(1, i)**2 - r(2, i)**2)
800 122337 : y(i) = pf*CMPLX(yp, -ym, KIND=dp)
801 : END DO
802 : END SELECT
803 : CASE DEFAULT
804 15757 : IF (m < -l .OR. m > l) CPABORT("m value out of bounds")
805 15757 : lpm = fac(l + ABS(m))
806 15757 : lmm = fac(l - ABS(m))
807 15757 : t = 4.0_dp*pi
808 15757 : IF (ABS(lpm) < EPSILON(1.0_dp)) THEN
809 0 : pf = REAL(2*l + 1, KIND=dp)/t
810 : ELSE
811 15757 : pf = (REAL(2*l + 1, KIND=dp)*lmm)/(t*lpm)
812 : END IF
813 15757 : pf = SQRT(pf)
814 9331787 : DO i = 1, n
815 9296630 : z = r(3, i)
816 9296630 : plm = legendre(z, l, m)
817 9312387 : IF (m == 0) THEN
818 811250 : y(i) = pf*plm
819 : ELSE
820 8485380 : rxy = SQRT(r(1, i)**2 + r(2, i)**2)
821 8485380 : IF (rxy < EPSILON(1.0_dp)) THEN
822 28764 : y(i) = 0.0_dp
823 : ELSE
824 8456616 : cp = r(1, i)/rxy
825 8456616 : sp = r(2, i)/rxy
826 8456616 : IF (m > 0) THEN
827 4228308 : yp = cosn(m, cp, sp)
828 4228308 : ym = sinn(m, cp, sp)
829 4228308 : y(i) = pf*plm*CMPLX(yp, ym, KIND=dp)
830 : ELSE
831 4228308 : yp = cosn(-m, cp, sp)
832 4228308 : ym = sinn(-m, cp, sp)
833 4228308 : y(i) = pf*plm*CMPLX(yp, -ym, KIND=dp)
834 : END IF
835 : END IF
836 : END IF
837 : END DO
838 : END SELECT
839 :
840 19400 : END SUBROUTINE cvy_lm
841 :
842 : ! **************************************************************************************************
843 : !> \brief ...
844 : !> \param r ...
845 : !> \param y ...
846 : !> \param l ...
847 : !> \param m ...
848 : ! **************************************************************************************************
849 0 : SUBROUTINE ccy_lm(r, y, l, m)
850 : !
851 : ! Complex Spherical Harmonics
852 : ! _ _
853 : ! | [(2l+1)(l-|m|)!] |1/2 m
854 : ! Y_lm ( t, p ) = |------------------| P_l(cos(t)) Exp[ i m p ]
855 : ! | [4Pi(l+|m|)!]| |
856 : !
857 : ! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1
858 : !
859 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r
860 : COMPLEX(KIND=dp), INTENT(OUT) :: y
861 : INTEGER, INTENT(IN) :: l, m
862 :
863 : REAL(KIND=dp) :: cp, lmm, lpm, pf, plm, rxy, sp, t, ym, &
864 : yp, z
865 :
866 0 : SELECT CASE (l)
867 : CASE (:-1)
868 0 : CPABORT("Negative l value")
869 : CASE (0)
870 0 : pf = SQRT(1.0_dp/(4.0_dp*pi))
871 0 : IF (m /= 0) CPABORT("l = 0 and m value out of bounds")
872 0 : y = pf
873 : CASE (1)
874 0 : SELECT CASE (m)
875 : CASE DEFAULT
876 0 : CPABORT("l = 1 and m value out of bounds")
877 : CASE (1)
878 0 : pf = SQRT(3.0_dp/(8.0_dp*pi))
879 0 : yp = r(1)
880 0 : ym = r(2)
881 0 : y = pf*CMPLX(yp, ym, KIND=dp)
882 : CASE (0)
883 0 : pf = SQRT(3.0_dp/(4.0_dp*pi))
884 0 : y = pf*r(3)
885 : CASE (-1)
886 0 : pf = SQRT(3.0_dp/(8.0_dp*pi))
887 0 : yp = r(1)
888 0 : ym = r(2)
889 0 : y = pf*CMPLX(yp, -ym, KIND=dp)
890 : END SELECT
891 : CASE (2)
892 0 : SELECT CASE (m)
893 : CASE DEFAULT
894 0 : CPABORT("l = 2 and m value out of bounds")
895 : CASE (2)
896 0 : pf = SQRT(15.0_dp/(32.0_dp*pi))
897 0 : yp = (r(1)*r(1) - r(2)*r(2))
898 0 : ym = 2.0_dp*r(1)*r(2)
899 0 : y = pf*CMPLX(yp, ym, KIND=dp)
900 : CASE (1)
901 0 : pf = SQRT(15.0_dp/(8.0_dp*pi))
902 0 : yp = r(3)*r(1)
903 0 : ym = r(3)*r(2)
904 0 : y = pf*CMPLX(yp, ym, KIND=dp)
905 : CASE (0)
906 0 : pf = SQRT(5.0_dp/(16.0_dp*pi))
907 0 : y = pf*(3.0_dp*r(3)*r(3) - 1.0_dp)
908 : CASE (-1)
909 0 : pf = SQRT(15.0_dp/(8.0_dp*pi))
910 0 : yp = r(3)*r(1)
911 0 : ym = r(3)*r(2)
912 0 : y = pf*CMPLX(yp, -ym, KIND=dp)
913 : CASE (-2)
914 0 : pf = SQRT(15.0_dp/(32.0_dp*pi))
915 0 : yp = (r(1)*r(1) - r(2)*r(2))
916 0 : ym = 2.0_dp*r(1)*r(2)
917 0 : y = pf*CMPLX(yp, -ym, KIND=dp)
918 : END SELECT
919 : CASE (3)
920 0 : SELECT CASE (m)
921 : CASE DEFAULT
922 0 : CPABORT("l = 3 and m value out of bounds")
923 : CASE (3)
924 0 : pf = SQRT(35.0_dp/(64.0_dp*pi))
925 0 : yp = r(1)*(r(1)**2 - 3.0_dp*r(2)**2)
926 0 : ym = r(2)*(3.0_dp*r(1)**2 - r(2)**2)
927 0 : y = pf*CMPLX(yp, ym, KIND=dp)
928 : CASE (2)
929 0 : pf = SQRT(105.0_dp/(32.0_dp*pi))
930 0 : yp = r(3)*(r(1)**2 - r(2)**2)
931 0 : ym = 2.0_dp*r(1)*r(2)*r(3)
932 0 : y = pf*CMPLX(yp, ym, KIND=dp)
933 : CASE (1)
934 0 : pf = SQRT(21.0_dp/(64.0_dp*pi))
935 0 : yp = r(1)*(5.0_dp*r(3)*r(3) - 1.0_dp)
936 0 : ym = r(2)*(5.0_dp*r(3)*r(3) - 1.0_dp)
937 0 : y = pf*CMPLX(yp, ym, KIND=dp)
938 : CASE (0)
939 0 : pf = SQRT(7.0_dp/(16.0_dp*pi))
940 0 : y = pf*r(3)*(5.0_dp*r(3)*r(3) - 3.0_dp)
941 : CASE (-1)
942 0 : pf = SQRT(21.0_dp/(64.0_dp*pi))
943 0 : yp = r(1)*(5.0_dp*r(3)*r(3) - 1.0_dp)
944 0 : ym = r(2)*(5.0_dp*r(3)*r(3) - 1.0_dp)
945 0 : y = pf*CMPLX(yp, -ym, KIND=dp)
946 : CASE (-2)
947 0 : pf = SQRT(105.0_dp/(32.0_dp*pi))
948 0 : yp = r(3)*(r(1)**2 - r(2)**2)
949 0 : ym = 2.0_dp*r(1)*r(2)*r(3)
950 0 : y = pf*CMPLX(yp, -ym, KIND=dp)
951 : CASE (-3)
952 0 : pf = SQRT(35.0_dp/(64.0_dp*pi))
953 0 : yp = r(1)*(r(1)**2 - 3.0_dp*r(2)**2)
954 0 : ym = r(2)*(3.0_dp*r(1)**2 - r(2)**2)
955 0 : y = pf*CMPLX(yp, -ym, KIND=dp)
956 : END SELECT
957 : CASE DEFAULT
958 0 : IF (m < -l .OR. m > l) CPABORT("m value out of bounds")
959 0 : lpm = fac(l + ABS(m))
960 0 : lmm = fac(l - ABS(m))
961 0 : t = 4.0_dp*pi
962 0 : IF (ABS(lpm) < EPSILON(1.0_dp)) THEN
963 0 : pf = REAL(2*l + 1, KIND=dp)/t
964 : ELSE
965 0 : pf = (REAL(2*l + 1, KIND=dp)*lmm)/(t*lpm)
966 : END IF
967 0 : pf = SQRT(pf)
968 0 : z = r(3)
969 0 : plm = legendre(z, l, m)
970 0 : IF (m == 0) THEN
971 0 : y = pf*plm
972 : ELSE
973 0 : rxy = SQRT(r(1)**2 + r(2)**2)
974 0 : IF (rxy < EPSILON(1.0_dp)) THEN
975 0 : y = 0.0_dp
976 : ELSE
977 0 : cp = r(1)/rxy
978 0 : sp = r(2)/rxy
979 0 : IF (m > 0) THEN
980 0 : yp = cosn(m, cp, sp)
981 0 : ym = sinn(m, cp, sp)
982 0 : y = pf*plm*CMPLX(yp, ym, KIND=dp)
983 : ELSE
984 0 : yp = cosn(-m, cp, sp)
985 0 : ym = sinn(-m, cp, sp)
986 0 : y = pf*plm*CMPLX(yp, -ym, KIND=dp)
987 : END IF
988 : END IF
989 : END IF
990 : END SELECT
991 :
992 0 : END SUBROUTINE ccy_lm
993 :
994 : ! Calculation of derivatives of Spherical Harmonics
995 :
996 : ! **************************************************************************************************
997 : !> \brief ...
998 : !> \param c ...
999 : !> \param dy ...
1000 : !> \param l ...
1001 : !> \param m ...
1002 : ! **************************************************************************************************
1003 817784 : SUBROUTINE dry_lm(c, dy, l, m)
1004 : !
1005 : ! Real Spherical Harmonics
1006 : ! _ _
1007 : ! | [(2l+1)(l-|m|)!] |1/2 m cos(m p) m>=0
1008 : ! Y_lm ( t, p ) = |---------------------| P_l(cos(t))
1009 : ! |[2Pi(1+d_m0)(l+|m|)!]| sin(|m| p) m<0
1010 : !
1011 : ! Input: c == (t,p)
1012 : ! Output: dy == (dy/dt, dy/dp)
1013 : !
1014 : ! x == sin(t)*cos(p)
1015 : ! y == sin(t)*sin(p)
1016 : ! z == cos(t)
1017 : !
1018 :
1019 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: c
1020 : REAL(KIND=dp), DIMENSION(2), INTENT(OUT) :: dy
1021 : INTEGER, INTENT(IN) :: l, m
1022 :
1023 : REAL(KIND=dp) :: cp, ct, dplm, lmm, lpm, p, pf, rxy, sp, &
1024 : st, t, tt, y, z
1025 : REAL(KIND=dp), DIMENSION(3) :: r
1026 :
1027 817784 : t = c(1)
1028 817784 : ct = COS(t)
1029 817784 : st = SIN(t)
1030 817784 : p = c(2)
1031 817784 : cp = COS(p)
1032 817784 : sp = SIN(p)
1033 817784 : r(1) = st*cp
1034 817784 : r(2) = st*sp
1035 817784 : r(3) = ct
1036 :
1037 : ! dY/dp
1038 817784 : IF (m == 0) THEN
1039 287704 : dy(2) = 0.0_dp
1040 : ELSE
1041 530080 : CALL y_lm(r, y, l, -m)
1042 530080 : dy(2) = -REAL(m, KIND=dp)*y
1043 : END IF
1044 :
1045 : ! dY/dt
1046 817784 : SELECT CASE (l)
1047 : CASE (:-1)
1048 0 : CPABORT("Negative l value")
1049 : CASE (0)
1050 116972 : IF (m /= 0) CPABORT("l = 0 and m value out of bounds")
1051 116972 : dy(1) = 0.0_dp
1052 : CASE (1)
1053 292940 : SELECT CASE (m)
1054 : CASE DEFAULT
1055 0 : CPABORT("l = 1 and m value out of bounds")
1056 : CASE (1)
1057 98416 : pf = SQRT(3.0_dp/(4.0_dp*pi))
1058 98416 : dy(1) = pf*ct*cp
1059 : CASE (0)
1060 98416 : pf = SQRT(3.0_dp/(4.0_dp*pi))
1061 98416 : dy(1) = -pf*st
1062 : CASE (-1)
1063 98416 : pf = SQRT(3.0_dp/(4.0_dp*pi))
1064 295248 : dy(1) = pf*ct*sp
1065 : END SELECT
1066 : CASE (2)
1067 57064 : SELECT CASE (m)
1068 : CASE DEFAULT
1069 0 : CPABORT("l = 2 and m value out of bounds")
1070 : CASE (2)
1071 58588 : pf = SQRT(15.0_dp/(16.0_dp*pi))
1072 58588 : dy(1) = pf*2.0_dp*st*ct*COS(2._dp*p)
1073 : CASE (1)
1074 58588 : pf = SQRT(15.0_dp/(4.0_dp*pi))
1075 58588 : dy(1) = pf*cp*(ct*ct - st*st)
1076 : CASE (0)
1077 58588 : pf = SQRT(5.0_dp/(16.0_dp*pi))
1078 58588 : dy(1) = -pf*6.0_dp*ct*st
1079 : CASE (-1)
1080 58588 : pf = SQRT(15.0_dp/(4.0_dp*pi))
1081 58588 : dy(1) = pf*sp*(ct*ct - st*st)
1082 : CASE (-2)
1083 58588 : pf = SQRT(15.0_dp/(16.0_dp*pi))
1084 292940 : dy(1) = pf*2.0_dp*st*ct*SIN(2._dp*p)
1085 : END SELECT
1086 : CASE (3)
1087 55560 : SELECT CASE (m)
1088 : CASE DEFAULT
1089 0 : CPABORT("l = 3 and m value out of bounds")
1090 : CASE (3)
1091 8152 : pf = SQRT(35.0_dp/(32.0_dp*pi))
1092 8152 : dy(1) = pf*3.0_dp*COS(3._dp*p)*ct*st*st
1093 : CASE (2)
1094 8152 : pf = SQRT(105.0_dp/(16.0_dp*pi))
1095 8152 : dy(1) = pf*2.0_dp*COS(2._dp*p)*ct*st
1096 : CASE (1)
1097 8152 : pf = SQRT(21.0_dp/(32.0_dp*pi))
1098 8152 : dy(1) = pf*cp*(ct*(5.0_dp*ct - 1.0_dp) - 5.0_dp*st*st)
1099 : CASE (0)
1100 8152 : pf = SQRT(7.0_dp/(16.0_dp*pi))
1101 8152 : dy(1) = pf*r(3)*(3.0_dp - 15.0_dp*ct*ct)*st
1102 : CASE (-1)
1103 8152 : pf = SQRT(21.0_dp/(32.0_dp*pi))
1104 8152 : dy(1) = pf*sp*(ct*(5.0_dp*ct - 1.0_dp) - 5.0_dp*st*st)
1105 : CASE (-2)
1106 8152 : pf = SQRT(105.0_dp/(16.0_dp*pi))
1107 8152 : dy(1) = pf*2.0_dp*SIN(2._dp*p)*ct*st
1108 : CASE (-3)
1109 8152 : pf = SQRT(35.0_dp/(32.0_dp*pi))
1110 57064 : dy(1) = pf*3.0_dp*SIN(3._dp*p)*ct*st*st
1111 : END SELECT
1112 : CASE DEFAULT
1113 55560 : IF (m < -l .OR. m > l) CPABORT("m value out of bounds")
1114 55560 : lpm = fac(l + ABS(m))
1115 55560 : lmm = fac(l - ABS(m))
1116 : IF (m == 0) THEN
1117 : tt = 4.0_dp*pi
1118 : ELSE
1119 : tt = 2.0_dp*pi
1120 : END IF
1121 : IF (ABS(lpm) < EPSILON(1.0_dp)) THEN
1122 : pf = REAL(2*l + 1, KIND=dp)/tt
1123 : ELSE
1124 55560 : pf = (REAL(2*l + 1, KIND=dp)*lmm)/(tt*lpm)
1125 : END IF
1126 55560 : pf = SQRT(pf)
1127 55560 : z = ct
1128 55560 : dplm = dlegendre(z, l, m)
1129 873344 : IF (m == 0) THEN
1130 : y = pf*dplm
1131 : ELSE
1132 49984 : rxy = SQRT(r(1)**2 + r(2)**2)
1133 : IF (rxy < EPSILON(1.0_dp)) THEN
1134 : y = 0.0_dp
1135 : ELSE
1136 : IF (m > 0) THEN
1137 : y = pf*dplm*cosn(m, cp, sp)
1138 : ELSE
1139 : y = pf*dplm*sinn(-m, cp, sp)
1140 : END IF
1141 : END IF
1142 : END IF
1143 : END SELECT
1144 :
1145 817784 : END SUBROUTINE dry_lm
1146 :
1147 : ! **************************************************************************************************
1148 : !> \brief ...
1149 : !> \param c ...
1150 : !> \param dy ...
1151 : !> \param l ...
1152 : !> \param m ...
1153 : ! **************************************************************************************************
1154 0 : SUBROUTINE dcy_lm(c, dy, l, m)
1155 : !
1156 : ! Complex Spherical Harmonics
1157 : ! _ _
1158 : ! | [(2l+1)(l-|m|)!] |1/2 m
1159 : ! Y_lm ( t, p ) = |------------------| P_l(cos(t)) Exp[ i m p ]
1160 : ! | [4Pi(l+|m|)!]| |
1161 : !
1162 : ! Input: r == (x,y,z) : normalised x^2 + y^2 + z^2 = 1
1163 : !
1164 : ! Input: c == (t,p)
1165 : ! Output: dy == (dy/dt, dy/dp)
1166 : !
1167 : ! x == sin(t)*cos(p)
1168 : ! y == sin(t)*sin(p)
1169 : ! z == cos(t)
1170 : !
1171 :
1172 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: c
1173 : COMPLEX(KIND=dp), DIMENSION(2), INTENT(OUT) :: dy
1174 : INTEGER, INTENT(IN) :: l, m
1175 :
1176 : REAL(KIND=dp), DIMENSION(2) :: dd
1177 :
1178 0 : CPABORT("Not implemented")
1179 0 : CALL dry_lm(c, dd, l, m)
1180 0 : dy(1) = CMPLX(dd(1), 0.0_dp, KIND=dp)
1181 0 : dy(2) = CMPLX(dd(2), 0.0_dp, KIND=dp)
1182 :
1183 0 : END SUBROUTINE dcy_lm
1184 :
1185 : ! **************************************************************************************************
1186 : !> \brief ...
1187 : !> \param x ...
1188 : !> \param l ...
1189 : !> \param m ...
1190 : !> \return ...
1191 : ! **************************************************************************************************
1192 26256040 : FUNCTION legendre(x, l, m) RESULT(plm)
1193 :
1194 : REAL(KIND=dp), INTENT(IN) :: x
1195 : INTEGER, INTENT(IN) :: l, m
1196 : REAL(KIND=dp) :: plm
1197 :
1198 : INTEGER :: il, im, mm
1199 : REAL(KIND=dp) :: fact, pll, pmm, pmmp1, somx2
1200 :
1201 26256040 : IF (ABS(x) > 1.0_dp) CPABORT("x value > 1")
1202 26256040 : SELECT CASE (l)
1203 : CASE (:-1)
1204 0 : CPABORT("Negative l value")
1205 : CASE (0)
1206 198 : plm = 1.0_dp
1207 : CASE (1)
1208 396 : SELECT CASE (ABS(m))
1209 : CASE DEFAULT
1210 0 : CPABORT("l = 1 and m value out of bounds")
1211 : CASE (1)
1212 0 : plm = SQRT(1.0_dp - x*x)
1213 : CASE (0)
1214 198 : plm = x
1215 : END SELECT
1216 : CASE (2)
1217 396 : SELECT CASE (ABS(m))
1218 : CASE DEFAULT
1219 0 : CPABORT("l = 2 and m value out of bounds")
1220 : CASE (2)
1221 0 : plm = 3.0_dp*(1.0_dp - x*x)
1222 : CASE (1)
1223 0 : plm = 3.0_dp*x*SQRT(1.0_dp - x*x)
1224 : CASE (0)
1225 198 : plm = 1.5_dp*x*x - 0.5_dp
1226 : END SELECT
1227 : CASE (3)
1228 3686874 : SELECT CASE (ABS(m))
1229 : CASE DEFAULT
1230 0 : CPABORT("l = 3 and m value out of bounds")
1231 : CASE (3)
1232 0 : plm = 15.0_dp*(1.0_dp - x*x)**1.5_dp
1233 : CASE (2)
1234 0 : plm = 15.0_dp*x*(1.0_dp - x*x)
1235 : CASE (1)
1236 0 : plm = (7.5_dp*x*x - 1.5_dp)*SQRT(1.0_dp - x*x)
1237 : CASE (0)
1238 198 : plm = 2.5_dp*x**3 - 1.5_dp*x
1239 : END SELECT
1240 : CASE (4)
1241 7063400 : SELECT CASE (ABS(m))
1242 : CASE DEFAULT
1243 0 : CPABORT("l = 4 and m value out of bounds")
1244 : CASE (4)
1245 745432 : plm = 105.0_dp*(1.0_dp - x*x)**2
1246 : CASE (3)
1247 792632 : plm = 105.0_dp*x*(1.0_dp - x*x)**1.5_dp
1248 : CASE (2)
1249 842192 : plm = (52.5_dp*x*x - 7.5_dp)*(1.0_dp - x*x)
1250 : CASE (1)
1251 856352 : plm = (17.5_dp*x**3 - 7.5_dp*x)*SQRT(1.0_dp - x*x)
1252 : CASE (0)
1253 3686676 : plm = 4.375_dp*x**4 - 3.75_dp*x**2 + 0.375_dp
1254 : END SELECT
1255 : CASE (5)
1256 7644020 : SELECT CASE (ABS(m))
1257 : CASE DEFAULT
1258 0 : CPABORT("l = 5 and m value out of bounds")
1259 : CASE (5)
1260 510384 : plm = 945.0_dp*(1.0_dp - x*x)**2.5_dp
1261 : CASE (4)
1262 564664 : plm = 945.0_dp*x*(1.0_dp - x*x)**2
1263 : CASE (3)
1264 621304 : plm = -(-472.5_dp*x*x + 52.5_dp)*(1.0_dp - x*x)**1.5_dp
1265 : CASE (2)
1266 649624 : plm = (157.5_dp*x**3 - 52.5_dp*x)*(1.0_dp - x*x)
1267 : CASE (1)
1268 : plm = -(-39.375_dp*x**4 + 26.25_dp*x*x - &
1269 680304 : 1.875_dp)*SQRT(1.0_dp - x*x)
1270 : CASE (0)
1271 3376724 : plm = 7.875_dp*x**5 - 8.75_dp*x**3 + 1.875_dp*x
1272 : END SELECT
1273 : CASE (6)
1274 19191972 : SELECT CASE (ABS(m))
1275 : CASE DEFAULT
1276 0 : CPABORT("l = 6 and m value out of bounds")
1277 : CASE (6)
1278 522804 : plm = 10395.0_dp*(1.0_dp - x*x)**3
1279 : CASE (5)
1280 580624 : plm = 10395.0_dp*x*(1.0_dp - x*x)**2.5_dp
1281 : CASE (4)
1282 640804 : plm = (5197.5_dp*x*x - 472.5_dp)*(1.0_dp - x*x)**2
1283 : CASE (3)
1284 : plm = -(-1732.5_dp*x**3 + 472.5_dp*x)* &
1285 679744 : (1.0_dp - x*x)**1.5_dp
1286 : CASE (2)
1287 : plm = (433.125_dp*x**4 - 236.25_dp*x**2 + &
1288 721044 : 13.125_dp)*(1.0_dp - x*x)
1289 : CASE (1)
1290 : plm = -(-86.625_dp*x**5 + 78.75_dp*x**3 - &
1291 734024 : 13.125_dp*x)*SQRT(1.0_dp - x*x)
1292 : CASE (0)
1293 : plm = 14.4375_dp*x**6 - 19.6875_dp*x**4 + &
1294 4267296 : 6.5625_dp*x**2 - 0.3125_dp
1295 : END SELECT
1296 : CASE DEFAULT
1297 14924676 : mm = ABS(m)
1298 14924676 : IF (mm > l) CPABORT("m out of bounds")
1299 : ! use recurence from numerical recipies
1300 14924676 : pmm = 1.0_dp
1301 14924676 : IF (mm > 0) THEN
1302 13762024 : somx2 = SQRT((1.0_dp - x)*(1.0_dp + x))
1303 13762024 : fact = 1.0_dp
1304 71397248 : DO im = 1, mm
1305 57635224 : pmm = pmm*fact*somx2
1306 71397248 : fact = fact + 2.0_dp
1307 : END DO
1308 : END IF
1309 41180716 : IF (l == mm) THEN
1310 : plm = pmm
1311 : ELSE
1312 14154612 : pmmp1 = x*REAL(2*mm + 1, KIND=dp)*pmm
1313 14154612 : IF (l == mm + 1) THEN
1314 : plm = pmmp1
1315 : ELSE
1316 78758448 : DO il = mm + 2, l
1317 : pll = (x*REAL(2*il - 1, KIND=dp)*pmmp1 - &
1318 65595740 : REAL(il + mm - 1, KIND=dp)*pmm)/REAL(il - mm, KIND=dp)
1319 65595740 : pmm = pmmp1
1320 78758448 : pmmp1 = pll
1321 : END DO
1322 : plm = pll
1323 : END IF
1324 : END IF
1325 : END SELECT
1326 :
1327 26256040 : END FUNCTION legendre
1328 :
1329 : ! **************************************************************************************************
1330 : !> \brief ...
1331 : !> \param x ...
1332 : !> \param l ...
1333 : !> \param m ...
1334 : !> \return ...
1335 : ! **************************************************************************************************
1336 547436 : FUNCTION dlegendre(x, l, m) RESULT(dplm)
1337 : REAL(KIND=dp), INTENT(IN) :: x
1338 : INTEGER, INTENT(IN) :: l, m
1339 : REAL(KIND=dp) :: dplm
1340 :
1341 : INTEGER :: mm
1342 :
1343 547436 : IF (ABS(x) > 1.0_dp) CPABORT("x value > 1")
1344 547436 : SELECT CASE (l)
1345 : CASE (0)
1346 124 : dplm = 0.0_dp
1347 : CASE (1)
1348 248 : SELECT CASE (ABS(m))
1349 : CASE DEFAULT
1350 0 : CPABORT("l = 1 and m value out of bounds")
1351 : CASE (1)
1352 0 : dplm = -x/SQRT(1.0_dp - x*x)
1353 : CASE (0)
1354 124 : dplm = 1.0_dp
1355 : END SELECT
1356 : CASE (2)
1357 248 : SELECT CASE (ABS(m))
1358 : CASE DEFAULT
1359 0 : CPABORT("l = 2 and m value out of bounds")
1360 : CASE (2)
1361 0 : dplm = -6.0_dp*x
1362 : CASE (1)
1363 0 : dplm = 3.0_dp*SQRT(1.0_dp - x*x) - 3.0_dp*x*x/SQRT(1.0_dp - x*x)
1364 : CASE (0)
1365 124 : dplm = 3.0_dp*x
1366 : END SELECT
1367 : CASE (3)
1368 26116 : SELECT CASE (ABS(m))
1369 : CASE DEFAULT
1370 0 : CPABORT("l = 3 and m value out of bounds")
1371 : CASE (3)
1372 0 : dplm = -45.0_dp*SQRT(1.0_dp - x*x)*x
1373 : CASE (2)
1374 0 : dplm = 15.0_dp*(1.0_dp - x*x) - 30.0_dp*x*x
1375 : CASE (1)
1376 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)
1377 : CASE (0)
1378 124 : dplm = 7.5_dp*x*x - 1.5_dp
1379 : END SELECT
1380 : CASE (4)
1381 55560 : SELECT CASE (ABS(m))
1382 : CASE DEFAULT
1383 0 : CPABORT("l = 4 and m value out of bounds")
1384 : CASE (4)
1385 5776 : dplm = -420*x*(1 - x*x)
1386 : CASE (3)
1387 5776 : 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)
1388 : CASE (2)
1389 5776 : dplm = 105.0_dp*x*(1.0_dp - x*x) - 2.0_dp*x*(52.5_dp*x*x - 7.5_dp)
1390 : CASE (1)
1391 5776 : IF (ABS(x) - 1.0_dp < EPSILON(1.0_dp)) THEN
1392 : dplm = 0.0_dp
1393 : ELSE
1394 : dplm = (17.5_dp*3.0_dp*x**2 - 7.5_dp)*SQRT(1.0_dp - x*x) - &
1395 0 : x*(17.5_dp*x**3 - 7.5_dp*x)/SQRT(1.0_dp - x*x)
1396 : END IF
1397 : CASE (0)
1398 25992 : dplm = 4.375_dp*4.0_dp*x**3 - 3.75_dp*2.0_dp*x
1399 : END SELECT
1400 : CASE (5)
1401 521072 : SELECT CASE (ABS(m))
1402 : CASE DEFAULT
1403 0 : CPABORT("l = 5 and m value out of bounds")
1404 : CASE (5)
1405 5376 : dplm = -945.0_dp*5.0_dp*x*(1.0_dp - x*x)**1.5_dp
1406 : CASE (4)
1407 5376 : dplm = 945.0_dp*((1.0_dp - x*x)**2 - 2.0_dp*x*x*(1.0_dp - x*x))
1408 : CASE (3)
1409 : dplm = 945.0_dp*x*(1.0_dp - x*x)**1.5_dp - &
1410 5376 : 3.0_dp*x*(472.5_dp*x*x - 52.5_dp)*(1.0_dp - x*x)**0.5_dp
1411 : CASE (2)
1412 : dplm = (3.0_dp*157.5_dp*x**2 - 52.5_dp)*(1.0_dp - x*x) - &
1413 5376 : (157.5_dp*x**3 - 52.5_dp*x)*(-2.0_dp*x)
1414 : CASE (1)
1415 5376 : IF (ABS(x) - 1.0_dp < EPSILON(1.0_dp)) THEN
1416 : dplm = 0.0_dp
1417 : ELSE
1418 : dplm = -(-39.375_dp*4.0_dp*x*x*x + 2.0_dp*26.25_dp*x)*SQRT(1.0_dp - x*x) + &
1419 0 : x*(-39.375_dp*x**4 + 26.25_dp*x*x - 1.875_dp)/SQRT(1.0_dp - x*x)
1420 : END IF
1421 : CASE (0)
1422 29568 : dplm = 5.0_dp*7.875_dp*x**4 - 3.0_dp*8.75_dp*x**2 + 1.875_dp
1423 : END SELECT
1424 : CASE (6)
1425 491504 : SELECT CASE (ABS(m))
1426 : CASE DEFAULT
1427 0 : CPABORT("l = 6 and m value out of bounds")
1428 : CASE (6)
1429 75616 : dplm = -10395.0_dp*6.0_dp*x*(1.0_dp - x*x)**2
1430 : CASE (5)
1431 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)
1432 : CASE (4)
1433 : dplm = 2.0_dp*5197.5_dp*x*(1.0_dp - x*x)**2 - &
1434 75616 : 4.0_dp*x*(5197.5_dp*x*x - 472.5_dp)*(1.0_dp - x*x)
1435 : CASE (3)
1436 : dplm = -(-3.0_dp*1732.5_dp*x*x + 472.5_dp)*(1.0_dp - x*x)**1.5_dp + &
1437 75616 : (-1732.5_dp*x**3 + 472.5_dp*x)*3.0_dp*x*(1.0_dp - x*x)**0.5_dp
1438 : CASE (2)
1439 : dplm = (433.125_dp*4.0_dp*x**3 - 2.0_dp*236.25_dp*x)*(1.0_dp - x*x) - &
1440 75616 : 2.0_dp*x*(433.125_dp*x**4 - 236.25_dp*x**2 + 13.125_dp)
1441 : CASE (1)
1442 75616 : IF (ABS(x) - 1.0_dp < EPSILON(1.0_dp)) THEN
1443 : dplm = 0.0_dp
1444 : ELSE
1445 : 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) + &
1446 0 : x*(-86.625_dp*x**5 + 78.75_dp*x**3 - 13.125_dp*x)/SQRT(1.0_dp - x*x)
1447 : END IF
1448 : CASE (0)
1449 : dplm = 14.4375_dp*6.0_dp*x**5 - 19.6875_dp*4.0_dp*x**3 + &
1450 491504 : 6.5625_dp*2.0_dp*x
1451 : END SELECT
1452 : CASE DEFAULT
1453 0 : mm = ABS(m)
1454 0 : IF (mm > l) CPABORT("m out of bounds")
1455 :
1456 : !From Wikipedia: dPlm(x) = -(l+1)x/(x^2-1)*Plm(x) + (l-m+1)/(x^2-1)Pl+1m(x)
1457 547436 : IF (ABS(x) - 1.0_dp < EPSILON(1.0_dp)) THEN
1458 : dplm = 0.0_dp
1459 : ELSE
1460 : dplm = -REAL(l + 1, dp)*x/(x**2 - 1.0_dp)*legendre(x, l, m) &
1461 0 : + REAL(l - m + 1, dp)/(x**2 - 1.0_dp)*legendre(x, l + 1, m)
1462 : END IF
1463 : END SELECT
1464 :
1465 547436 : END FUNCTION dlegendre
1466 :
1467 : ! **************************************************************************************************
1468 : !> \brief ...
1469 : !> \param x ...
1470 : !> \param l ...
1471 : !> \return ...
1472 : ! **************************************************************************************************
1473 0 : FUNCTION dPof1(x, l)
1474 :
1475 : REAL(KIND=dp), INTENT(IN) :: x
1476 : INTEGER, INTENT(IN) :: l
1477 : REAL(KIND=dp) :: dPof1
1478 :
1479 0 : IF (ABS(x) - 1.0_dp > EPSILON(1.0_dp)) THEN
1480 0 : CPABORT("|x| is not 1")
1481 : END IF
1482 0 : IF (x > 0.0_dp) THEN
1483 0 : SELECT CASE (l)
1484 : CASE (:-1)
1485 0 : CPABORT("Negative l value")
1486 : CASE (0)
1487 0 : dPof1 = 0.0_dp
1488 : CASE (1)
1489 0 : dPof1 = 1.0_dp
1490 : CASE (2)
1491 0 : dPof1 = 3.0_dp
1492 : CASE (3)
1493 0 : dPof1 = 6.0_dp
1494 : CASE (4)
1495 0 : dPof1 = 10.0_dp
1496 : CASE (5)
1497 0 : dPof1 = 15.0_dp
1498 : CASE (6)
1499 0 : dPof1 = 21.0_dp
1500 : CASE (7:)
1501 0 : CPABORT("Not implemented")
1502 : END SELECT
1503 : ELSE
1504 0 : SELECT CASE (l)
1505 : CASE (:-1)
1506 0 : CPABORT("Negative l value")
1507 : CASE (0)
1508 0 : dPof1 = 0.0_dp
1509 : CASE (1)
1510 0 : dPof1 = 1.0_dp
1511 : CASE (2)
1512 0 : dPof1 = -3.0_dp
1513 : CASE (3)
1514 0 : dPof1 = 6.0_dp
1515 : CASE (4)
1516 0 : dPof1 = -10.0_dp
1517 : CASE (5)
1518 0 : dPof1 = 15.0_dp
1519 : CASE (6)
1520 0 : dPof1 = -21.0_dp
1521 : CASE (7:)
1522 0 : CPABORT("Not implemented")
1523 : END SELECT
1524 : END IF
1525 :
1526 0 : END FUNCTION dPof1
1527 :
1528 : ! **************************************************************************************************
1529 : !> \brief ...
1530 : !> \param n ...
1531 : !> \param k ...
1532 : !> \return ...
1533 : ! **************************************************************************************************
1534 68808434 : FUNCTION choose(n, k)
1535 :
1536 : INTEGER, INTENT(IN) :: n, k
1537 : REAL(KIND=dp) :: choose
1538 :
1539 68808434 : IF (n >= k) THEN
1540 68808434 : choose = REAL(NINT(fac(n)/(fac(k)*fac(n - k))), KIND=dp)
1541 : ELSE
1542 : choose = 0.0_dp
1543 : END IF
1544 :
1545 68808434 : END FUNCTION choose
1546 :
1547 : ! **************************************************************************************************
1548 : !> \brief ...
1549 : !> \param n ...
1550 : !> \param c ...
1551 : !> \param s ...
1552 : !> \return ...
1553 : ! **************************************************************************************************
1554 15899184 : FUNCTION cosn(n, c, s)
1555 :
1556 : INTEGER, INTENT(IN) :: n
1557 : REAL(KIND=dp), INTENT(IN) :: c, s
1558 : REAL(KIND=dp) :: cosn
1559 :
1560 : INTEGER :: i, j
1561 :
1562 15899184 : cosn = 0.0_dp
1563 15899184 : IF (ABS(c) < EPSILON(1.0_dp) .OR. n == 0) THEN
1564 736772 : IF (MOD(n, 2) == 0) THEN
1565 337504 : cosn = (-1.0_dp)**INT(n/2)
1566 : ELSE
1567 : cosn = 0.0_dp
1568 : END IF
1569 : ELSE
1570 53789862 : DO i = n, 0, -2
1571 53789862 : IF (i == 0) THEN
1572 6921056 : j = n
1573 : IF (j == 0) THEN
1574 : cosn = cosn + choose(n, j)
1575 6921056 : ELSE IF (MOD(j/2, 2) == 0) THEN
1576 2830284 : cosn = cosn + choose(n, j)*s**j
1577 : ELSE
1578 4090772 : cosn = cosn - choose(n, j)*s**j
1579 : END IF
1580 : ELSE
1581 31706394 : j = n - i
1582 31706394 : IF (j == 0) THEN
1583 15162412 : cosn = cosn + choose(n, j)*c**i
1584 16543982 : ELSE IF (MOD(j/2, 2) == 0) THEN
1585 5149460 : cosn = cosn + choose(n, j)*c**i*s**j
1586 : ELSE
1587 11394522 : cosn = cosn - choose(n, j)*c**i*s**j
1588 : END IF
1589 : END IF
1590 : END DO
1591 : END IF
1592 :
1593 15899184 : END FUNCTION cosn
1594 :
1595 : ! **************************************************************************************************
1596 : !> \brief ...
1597 : !> \param n ...
1598 : !> \param c ...
1599 : !> \param s ...
1600 : !> \return ...
1601 : ! **************************************************************************************************
1602 0 : FUNCTION dcosn_dcp(n, c, s)
1603 :
1604 : INTEGER, INTENT(IN) :: n
1605 : REAL(KIND=dp), INTENT(IN) :: c, s
1606 : REAL(KIND=dp) :: dcosn_dcp
1607 :
1608 : INTEGER :: i, j
1609 :
1610 0 : dcosn_dcp = 0.0_dp
1611 :
1612 0 : IF (s < EPSILON(1.0_dp)) THEN
1613 : dcosn_dcp = 0.0_dp
1614 : ELSE
1615 0 : DO i = n, 0, -2
1616 0 : IF (i == 0) THEN
1617 : dcosn_dcp = dcosn_dcp
1618 : ELSE
1619 0 : j = n - i
1620 0 : IF (MOD(j/2, 2) == 0) THEN
1621 0 : dcosn_dcp = dcosn_dcp + choose(n, j)*REAL(i, dp)*c**(i - 1)*s**j
1622 : ELSE
1623 0 : dcosn_dcp = dcosn_dcp - choose(n, j)*REAL(i, dp)*c**(i - 1)*s**j
1624 : END IF
1625 : END IF
1626 : END DO
1627 : END IF
1628 :
1629 0 : END FUNCTION dcosn_dcp
1630 :
1631 : ! **************************************************************************************************
1632 : !> \brief ...
1633 : !> \param n ...
1634 : !> \param c ...
1635 : !> \param s ...
1636 : !> \return ...
1637 : ! **************************************************************************************************
1638 0 : FUNCTION dcosn_dsp(n, c, s)
1639 :
1640 : INTEGER, INTENT(IN) :: n
1641 : REAL(KIND=dp), INTENT(IN) :: c, s
1642 : REAL(KIND=dp) :: dcosn_dsp
1643 :
1644 : INTEGER :: i, j
1645 :
1646 0 : dcosn_dsp = 0.0_dp
1647 0 : IF (c < EPSILON(1.0_dp) .OR. s < EPSILON(1.0_dp)) THEN
1648 : dcosn_dsp = 0.0_dp
1649 : ELSE
1650 0 : DO i = n, 0, -2
1651 0 : j = n - i
1652 0 : IF (j == 0) THEN
1653 : dcosn_dsp = dcosn_dsp
1654 : ELSE
1655 0 : IF (MOD(j/2, 2) == 0) THEN
1656 0 : dcosn_dsp = dcosn_dsp + choose(n, j)*REAL(j, dp)*s**(j - 1)*c**i
1657 : ELSE
1658 0 : dcosn_dsp = dcosn_dsp - choose(n, j)*REAL(j, dp)*s**(j - 1)*c**i
1659 : END IF
1660 : END IF
1661 : END DO
1662 : END IF
1663 :
1664 0 : END FUNCTION dcosn_dsp
1665 :
1666 : ! **************************************************************************************************
1667 : !> \brief ...
1668 : !> \param n ...
1669 : !> \param c ...
1670 : !> \param s ...
1671 : !> \return ...
1672 : ! **************************************************************************************************
1673 15899184 : FUNCTION sinn(n, c, s)
1674 :
1675 : INTEGER, INTENT(IN) :: n
1676 : REAL(KIND=dp), INTENT(IN) :: c, s
1677 : REAL(KIND=dp) :: sinn
1678 :
1679 : INTEGER :: i, j
1680 :
1681 15899184 : sinn = 0.0_dp
1682 :
1683 15899184 : IF (ABS(s) < EPSILON(1.0_dp) .OR. n == 0) THEN
1684 : sinn = 0.0_dp
1685 15162412 : ELSE IF (ABS(c) < EPSILON(1.0_dp)) THEN
1686 736772 : IF (MOD(n, 2) == 0) THEN
1687 : sinn = 0.0_dp
1688 : ELSE
1689 399268 : sinn = s*(-1.0_dp)**INT((n - 1)/2)
1690 : END IF
1691 : ELSE
1692 44606624 : DO i = n - 1, 0, -2
1693 44606624 : IF (i == 0) THEN
1694 7842088 : j = n
1695 7842088 : IF (MOD(j/2, 2) == 0) THEN
1696 4593728 : sinn = sinn + choose(n, j)*s**j
1697 : ELSE
1698 3248360 : sinn = sinn - choose(n, j)*s**j
1699 : END IF
1700 : ELSE
1701 22338896 : j = n - i
1702 22338896 : IF (MOD(j/2, 2) == 0) THEN
1703 14739904 : sinn = sinn + choose(n, j)*c**i*s**j
1704 : ELSE
1705 7598992 : sinn = sinn - choose(n, j)*c**i*s**j
1706 : END IF
1707 : END IF
1708 : END DO
1709 : END IF
1710 :
1711 15899184 : END FUNCTION sinn
1712 :
1713 : ! **************************************************************************************************
1714 : !> \brief ...
1715 : !> \param n ...
1716 : !> \param c ...
1717 : !> \param s ...
1718 : !> \return ...
1719 : ! **************************************************************************************************
1720 0 : FUNCTION dsinn_dcp(n, c, s)
1721 :
1722 : INTEGER, INTENT(IN) :: n
1723 : REAL(KIND=dp), INTENT(IN) :: c, s
1724 : REAL(KIND=dp) :: dsinn_dcp
1725 :
1726 : INTEGER :: i, j
1727 :
1728 0 : dsinn_dcp = 0.0_dp
1729 :
1730 0 : IF (c < EPSILON(1.0_dp) .OR. s < EPSILON(1.0_dp)) THEN
1731 : dsinn_dcp = 0.0_dp
1732 : ELSE
1733 0 : DO i = n - 1, 0, -2
1734 0 : IF (i == 0) THEN
1735 : dsinn_dcp = dsinn_dcp
1736 : ELSE
1737 0 : j = n - i
1738 0 : IF (MOD(j/2, 2) == 0) THEN
1739 0 : dsinn_dcp = dsinn_dcp + choose(n, j)*REAL(i, dp)*c**(i - 1)*s**j
1740 : ELSE
1741 0 : dsinn_dcp = dsinn_dcp - choose(n, j)*REAL(i, dp)*c**(i - 1)*s**j
1742 : END IF
1743 : END IF
1744 : END DO
1745 : END IF
1746 :
1747 0 : END FUNCTION dsinn_dcp
1748 :
1749 : ! **************************************************************************************************
1750 : !> \brief ...
1751 : !> \param n ...
1752 : !> \param c ...
1753 : !> \param s ...
1754 : !> \return ...
1755 : ! **************************************************************************************************
1756 0 : FUNCTION dsinn_dsp(n, c, s)
1757 :
1758 : INTEGER, INTENT(IN) :: n
1759 : REAL(KIND=dp), INTENT(IN) :: c, s
1760 : REAL(KIND=dp) :: dsinn_dsp
1761 :
1762 : INTEGER :: i, j
1763 :
1764 0 : dsinn_dsp = 0.0_dp
1765 :
1766 0 : IF (c < EPSILON(1.0_dp) .OR. s < EPSILON(1.0_dp)) THEN
1767 : dsinn_dsp = 0.0_dp
1768 : ELSE
1769 0 : DO i = n - 1, 0, -2
1770 0 : j = n - i
1771 0 : IF (j == 0) THEN
1772 : dsinn_dsp = dsinn_dsp
1773 : ELSE
1774 0 : IF (MOD(j/2, 2) == 0) THEN
1775 0 : dsinn_dsp = dsinn_dsp + choose(n, j)*c**i*REAL(j, dp)*s**(j - 1)
1776 : ELSE
1777 0 : dsinn_dsp = dsinn_dsp - choose(n, j)*c**i*REAL(j, dp)*s**(j - 1)
1778 : END IF
1779 : END IF
1780 : END DO
1781 : END IF
1782 :
1783 0 : END FUNCTION dsinn_dsp
1784 :
1785 : ! **************************************************************************************************
1786 : !> \brief Compute the Clebsch-Gordon coefficient C = < j1 m1 j2 m2 | J M > = < j1 j2; m1 m2 | J M >
1787 : !> \param j1 Angular momentum quantum number of the first state | j1 m1 >
1788 : !> \param m1 Magnetic quantum number of the first first state | j1 m1 >
1789 : !> \param j2 Angular momentum quantum number of the second state | j2 m2 >
1790 : !> \param m2 Magnetic quantum number of the second state | j2 m2 >
1791 : !> \param J Angular momentum quantum number of the coupled state | J M >
1792 : !> \param M Magnetic quantum number of the coupled state | J M >
1793 : !> \param CG_coeff Clebsch-Gordon coefficient C^{JM}_{j1 m1 j2 m2}
1794 : !> \author Matthias Krack (16.09.2022, based on a program by D. G. Simpson)
1795 : !> \note Generic routine allowing also for fractional arguments. It should return CG coefficients
1796 : !> consistent with the standard definition and normalization, e.g. of Wolfram Mathematica.
1797 : !> The input parameters have to be integer or half-integer.
1798 : ! **************************************************************************************************
1799 0 : SUBROUTINE Clebsch_Gordon_coefficient(j1, m1, j2, m2, J, M, CG_coeff)
1800 :
1801 : REAL(KIND=dp), INTENT(IN) :: j1, m1, j2, m2, J, M
1802 : REAL(KIND=dp), INTENT(OUT) :: CG_coeff
1803 :
1804 : INTEGER :: k, kmax
1805 : REAL(KIND=dp) :: sumk, t
1806 :
1807 : ! Check validity of the input parameters
1808 0 : IF (j1 < 0.0_dp) &
1809 0 : CPABORT("The angular momentum quantum number j1 has to be nonnegative")
1810 0 : IF (.NOT. (is_integer(j1) .OR. is_integer(2.0_dp*j1))) &
1811 0 : CPABORT("The angular momentum quantum number j1 has to be integer or half-integer")
1812 0 : IF (j2 < 0.0_dp) &
1813 0 : CPABORT("The angular momentum quantum number j2 has to be nonnegative")
1814 0 : IF (.NOT. (is_integer(j2) .OR. is_integer(2.0_dp*j2))) &
1815 0 : CPABORT("The angular momentum quantum number j2 has to be integer or half-integer")
1816 0 : IF (J < 0.0_dp) &
1817 0 : CPABORT("The angular momentum quantum number J has to be nonnegative")
1818 0 : IF (.NOT. (is_integer(J) .OR. is_integer(2.0_dp*J))) &
1819 0 : CPABORT("The angular momentum quantum number J has to be integer or half-integer")
1820 0 : IF ((ABS(m1) - j1) > EPSILON(m1)) &
1821 0 : CPABORT("The angular momentum quantum number m1 has to satisfy -j1 <= m1 <= j1")
1822 0 : IF (.NOT. (is_integer(m1) .OR. is_integer(2.0_dp*m1))) &
1823 0 : CPABORT("The angular momentum quantum number m1 has to be integer or half-integer")
1824 0 : IF ((ABS(m2) - j2) > EPSILON(m2)) &
1825 0 : CPABORT("The angular momentum quantum number m2 has to satisfy -j2 <= m1 <= j2")
1826 0 : IF (.NOT. (is_integer(m2) .OR. is_integer(2.0_dp*m2))) &
1827 0 : CPABORT("The angular momentum quantum number m2 has to be integer or half-integer")
1828 0 : IF ((ABS(M) - J) > EPSILON(M)) &
1829 0 : CPABORT("The angular momentum quantum number M has to satisfy -J <= M <= J")
1830 0 : IF (.NOT. (is_integer(M) .OR. is_integer(2.0_dp*M))) &
1831 0 : CPABORT("The angular momentum quantum number M has to be integer or half-integer")
1832 :
1833 : IF (is_integer(j1 + j2 + J) .AND. &
1834 : is_integer(j1 + m1) .AND. &
1835 : is_integer(j2 + m2) .AND. &
1836 : is_integer(J + M) .AND. &
1837 0 : is_integer(J - j1 - m2) .AND. &
1838 : is_integer(J - j2 + m1)) THEN
1839 : IF ((J < ABS(j1 - j2)) .OR. &
1840 : (J > j1 + j2) .OR. &
1841 : (ABS(m1) > j1) .OR. &
1842 0 : (ABS(m2) > j2) .OR. &
1843 : (ABS(M) > J)) THEN
1844 0 : CG_coeff = 0.0_dp
1845 : ELSE
1846 : ! Compute Clebsch-Gordan coefficient
1847 0 : sumk = 0.0_dp
1848 0 : kmax = NINT(MAX(j1 + j2 - J, j1 - J + m2, j2 - J - m1, j1 - m1, j2 + m2))
1849 0 : DO k = 0, kmax
1850 0 : IF (j1 + j2 - J - k < 0.0_dp) CYCLE
1851 0 : IF (J - j1 - m2 + k < 0.0_dp) CYCLE
1852 0 : IF (J - j2 + m1 + k < 0.0_dp) CYCLE
1853 0 : IF (j1 - m1 - k < 0.0_dp) CYCLE
1854 0 : IF (j2 + m2 - k < 0.0_dp) CYCLE
1855 : t = sfac(NINT(j1 + j2 - J - k))*sfac(NINT(J - j1 - m2 + k))* &
1856 : sfac(NINT(J - j2 + m1 + k))*sfac(NINT(j1 - m1 - k))* &
1857 0 : sfac(NINT(j2 + m2 - k))*sfac(k)
1858 0 : IF (MOD(k, 2) == 1) t = -t
1859 0 : sumk = sumk + 1.0_dp/t
1860 : END DO
1861 : CG_coeff = SQRT((2.0_dp*J + 1)/sfac(NINT(j1 + j2 + J + 1.0_dp)))* &
1862 : SQRT(sfac(NINT(j1 + j2 - J))*sfac(NINT(j2 + J - j1))*sfac(NINT(J + j1 - j2)))* &
1863 : SQRT(sfac(NINT(j1 + m1))*sfac(NINT(j1 - m1))* &
1864 : sfac(NINT(j2 + m2))*sfac(NINT(j2 - m2))* &
1865 0 : sfac(NINT(J + M))*sfac(NINT(J - M)))*sumk
1866 : END IF
1867 : ELSE
1868 0 : CPABORT("Invalid set of input parameters < j1 m1 j2 m2 | J M > specified")
1869 : END IF
1870 :
1871 0 : END SUBROUTINE Clebsch_Gordon_coefficient
1872 :
1873 : ! **************************************************************************************************
1874 : !> \brief Compute the Wigner 3-j symbol
1875 : !> / j1 j2 j3 \
1876 : !> \ m1 m2 m3 /
1877 : !> using the Clebsch-Gordon coefficients
1878 : !> \param j1 Angular momentum quantum number of the first state | j1 m1 >
1879 : !> \param m1 Magnetic quantum number of the first first state | j1 m1 >
1880 : !> \param j2 Angular momentum quantum number of the second state | j2 m2 >
1881 : !> \param m2 Magnetic quantum number of the second state | j2 m2 >
1882 : !> \param j3 Angular momentum quantum number of the third state | j3 m3 >
1883 : !> \param m3 Magnetic quantum number of the third state | j3 m3 >
1884 : !> \param W_3j Wigner 3-j symbol
1885 : !> \author Matthias Krack (16.09.2022, MK)
1886 : ! **************************************************************************************************
1887 0 : SUBROUTINE Wigner_3j(j1, m1, j2, m2, j3, m3, W_3j)
1888 :
1889 : REAL(KIND=dp), INTENT(IN) :: j1, m1, j2, m2, j3, m3
1890 : REAL(KIND=dp), INTENT(OUT) :: W_3j
1891 :
1892 : REAL(KIND=dp) :: CG_coeff
1893 :
1894 : IF ((ABS(m1 + m2 + m3) < EPSILON(m1)) .AND. &
1895 0 : (ABS(j1 - j2) <= j3) .AND. (j3 <= ABS(j1 + j2)) .AND. &
1896 : is_integer(j1 + j2 + j3)) THEN
1897 0 : IF (is_integer(j1 - j2 - m3)) THEN
1898 0 : CALL Clebsch_Gordon_coefficient(j1, m1, j2, m2, j3, -m3, CG_coeff)
1899 0 : W_3j = (-1.0_dp)**INT(j1 - j2 - m3)*CG_coeff/SQRT(2.0_dp*j3 + 1.0_dp)
1900 : ELSE
1901 0 : CPABORT("j1 - j2 - m3 results in an invalid non-integer exponent")
1902 : END IF
1903 : ELSE
1904 0 : W_3j = 0.0_dp
1905 : END IF
1906 :
1907 0 : END SUBROUTINE Wigner_3j
1908 :
1909 : ! **************************************************************************************************
1910 : !> \brief Check if the input value is an integer number within double-precision accuracy
1911 : !> \param x Input value to be checked
1912 : !> \return True if the input value is integer otherwise false
1913 : !> \author Matthias Krack (16.09.2022, MK)
1914 : ! **************************************************************************************************
1915 0 : FUNCTION is_integer(x)
1916 :
1917 : REAL(KIND=dp), INTENT(IN) :: x
1918 : LOGICAL :: is_integer
1919 :
1920 0 : IF ((ABS(x) - AINT(ABS(x), KIND=dp)) > EPSILON(x)) THEN
1921 : is_integer = .FALSE.
1922 : ELSE
1923 0 : is_integer = .TRUE.
1924 : END IF
1925 :
1926 0 : END FUNCTION is_integer
1927 :
1928 : END MODULE spherical_harmonics
|
