(git:b279b6b)
spherical_harmonics.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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  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  l = l1 + l2
86  IF (l > lmax) CALL clebsch_gordon_init(l)
87  n = l/2 + 1
88  IF (n > SIZE(clm)) cpabort("Array too small")
89 
90  IF ((m1 >= 0 .AND. m2 >= 0) .OR. (m1 < 0 .AND. m2 < 0)) THEN
91  icase = 1
92  ELSE
93  icase = 2
94  END IF
95  ind = order(l1, m1, l2, m2)
96 
97  DO lp = mod(l, 2), l, 2
98  lm = lp/2 + 1
99  clm(lm) = cg_table(ind, lm, icase)
100  END DO
101 
102  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  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  l = l1 + l2
120  IF (l > lmax) CALL clebsch_gordon_init(l)
121  n = l/2 + 1
122  IF (n > SIZE(rlm, 1)) cpabort("Array too small")
123 
124  ind = order(l1, m1, l2, m2)
125  mm = getm(m1, m2)
126  IF ((m1 >= 0 .AND. m2 >= 0) .OR. (m1 < 0 .AND. m2 < 0)) THEN
127  icase1 = 1
128  icase2 = 2
129  ELSE
130  icase1 = 2
131  icase2 = 1
132  END IF
133 
134  DO lp = mod(l, 2), l, 2
135  lm = lp/2 + 1
136  xsi = get_factor(m1, m2, mm(1))
137  rlm(lm, 1) = xsi*cg_table(ind, lm, icase1)
138  xsi = get_factor(m1, m2, mm(2))
139  rlm(lm, 2) = xsi*cg_table(ind, lm, icase2)
140  END DO
141 
142  END SUBROUTINE clebsch_gordon_real
143 
144 ! **************************************************************************************************
145 !> \brief ...
146 !> \param m1 ...
147 !> \param m2 ...
148 !> \return ...
149 ! **************************************************************************************************
150  FUNCTION getm(m1, m2) RESULT(m)
151  INTEGER, INTENT(IN) :: m1, m2
152  INTEGER, DIMENSION(2) :: m
153 
154  INTEGER :: mm, mp
155 
156  mp = m1 + m2
157  mm = m1 - m2
158  IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
159  mp = -abs(mp)
160  mm = -abs(mm)
161  ELSE
162  mp = abs(mp)
163  mm = abs(mm)
164  END IF
165  m(1) = mp
166  m(2) = mm
167  END FUNCTION getm
168 
169 ! **************************************************************************************************
170 !> \brief ...
171 !> \param m1 ...
172 !> \param m2 ...
173 !> \param m ...
174 !> \return ...
175 ! **************************************************************************************************
176  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  f = 1.0_dp
183  IF (abs(m1) >= abs(m2)) THEN
184  mx = m1
185  my = m2
186  ELSE
187  mx = m2
188  my = m1
189  END IF
190  IF (mx*my == 0) THEN
191  f = 1.0_dp
192  ELSE IF (m == 0) THEN
193  IF (abs(mx) /= abs(my)) WRITE (*, '(A,3I6)') " 1) Illegal Case ", m1, m2, m
194  IF (mx*my > 0) THEN
195  f = 1.0_dp
196  ELSE
197  f = 0.0_dp
198  END IF
199  ELSE IF (abs(mx) + abs(my) == m) THEN
200  f = osq2
201  IF (mx < 0) f = -osq2
202  ELSE IF (abs(mx) + abs(my) == -m) THEN
203  f = osq2
204  ELSE IF (abs(mx) - abs(my) == -m) THEN
205  IF (mx*my > 0) WRITE (*, '(A,3I6)') " 2) Illegal Case ", m1, m2, m
206  IF (mx > 0) f = -osq2
207  IF (mx < 0) f = osq2
208  ELSE IF (abs(mx) - abs(my) == m) THEN
209  IF (mx*my < 0) WRITE (*, '(A,3I6)') " 3) Illegal Case ", m1, m2, m
210  f = osq2
211  ELSE
212  WRITE (*, '(A,3I6)') " 4) Illegal Case ", m1, m2, m
213  END IF
214  END FUNCTION get_factor
215 
216 ! **************************************************************************************************
217 !> \brief ...
218 ! **************************************************************************************************
219  SUBROUTINE clebsch_gordon_deallocate()
220  CHARACTER(len=*), PARAMETER :: routinen = 'clebsch_gordon_deallocate'
221 
222  INTEGER :: handle
223 
224  CALL timeset(routinen, handle)
225 
226  IF (ALLOCATED(cg_table)) THEN
227  DEALLOCATE (cg_table)
228  END IF
229  CALL timestop(handle)
230 
231  END SUBROUTINE clebsch_gordon_deallocate
232 
233 ! **************************************************************************************************
234 !> \brief ...
235 !> \param l ...
236 ! **************************************************************************************************
237  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  CALL timeset(routinen, handle)
246 
247  sq2 = sqrt(2.0_dp)
248  osq2 = 1.0_dp/sq2
249 
250  IF (l < 0) cpabort("l < 0")
251  IF (ALLOCATED(cg_table)) THEN
252  DEALLOCATE (cg_table)
253  END IF
254  ! maximum size of table
255  n = (l**4 + 6*l**3 + 15*l**2 + 18*l + 8)/8
256  m = l + 1
257  ALLOCATE (cg_table(n, m, 2))
258  lmax = l
259 
260  DO l1 = 0, lmax
261  DO m1 = 0, l1
262  iy = (l1*(l1 + 1))/2 + m1 + 1
263  DO l2 = l1, lmax
264  ml = 0
265  IF (l1 == l2) ml = m1
266  DO m2 = ml, l2
267  ix = (l2*(l2 + 1))/2 + m2 + 1
268  i1 = (ix*(ix - 1))/2 + iy
269  DO lp = mod(l1 + l2, 2), l1 + l2, 2
270  i2 = lp/2 + 1
271  mp = m2 + m1
272  cg_table(i1, i2, 1) = cgc(l1, m1, l2, m2, lp, mp)
273  mp = abs(m2 - m1)
274  IF (m2 >= m1) THEN
275  cg_table(i1, i2, 2) = cgc(l1, m1, lp, mp, l2, m2)
276  ELSE
277  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  CALL timestop(handle)
286 
287  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  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  IF (m1 < 0 .OR. m2 < 0 .OR. mp < 0) THEN
310  WRITE (*, *) l1, l2, lp
311  WRITE (*, *) m1, m2, mp
312  cpabort("Illegal input values")
313  END IF
314  IF (l2 < l1) THEN
315  la = l2
316  ma = m2
317  lb = l1
318  mb = m1
319  ELSE
320  la = l1
321  ma = m1
322  lb = l2
323  mb = m2
324  END IF
325 
326  IF (mod(la + lb + lp, 2) == 0 .AND. la + lb >= lp .AND. lp >= lb - la &
327  .AND. lb - mb >= 0) THEN
328  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  (sfac(la - ma)/sfac(la + ma))*(sfac(lb - mb)/sfac(lb + mb)))
332  s = (la + lb + lp)/2
333  tmin = max(0, -lb + la - mp)
334  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  (sfac(s)/sfac(2*s + 1))
338  f2 = 0.0_dp
339  DO t = tmin, tmax
340  z1 = lp + mp + t
341  z2 = la + lb - mp - t
342  z3 = lp - mp - t
343  z4 = lb - la + mp + t
344  z5 = la - ma - t
345  f2 = f2 + (-1)**t*(sfac(z1)/(sfac(t)*sfac(z3)))*(sfac(z2)/(sfac(z4)*sfac(z5)))
346  END DO
347  cgc = pref*f1*f2
348  ELSE
349  cgc = 0.0_dp
350  END IF
351 
352  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  FUNCTION sfac(n)
361  INTEGER, INTENT(IN) :: n
362  REAL(kind=dp) :: sfac
363 
364  INTEGER :: i
365 
366  IF (n > 170) THEN
367  cpabort("Factorials greater than 170! cannot be computed with double-precision")
368  ELSE IF (n > maxfac) THEN
369  sfac = fac(maxfac)
370  DO i = maxfac + 1, n
371  sfac = real(i, kind=dp)*sfac
372  END DO
373  ELSE IF (n >= 0) THEN
374  sfac = fac(n)
375  ELSE
376  sfac = 1.0_dp
377  END IF
378 
379  END FUNCTION sfac
380 
381 ! **************************************************************************************************
382 !> \brief ...
383 !> \param l1 ...
384 !> \param m1 ...
385 !> \param l2 ...
386 !> \param m2 ...
387 !> \return ...
388 ! **************************************************************************************************
389  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  i1 = (l1*(l1 + 1))/2 + abs(m1) + 1
396  i2 = (l2*(l2 + 1))/2 + abs(m2) + 1
397  ix = max(i1, i2)
398  iy = min(i1, i2)
399  ind = (ix*(ix - 1))/2 + iy
400  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  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  SELECT CASE (l)
429  CASE (:-1)
430  cpabort("Negative l value")
431  CASE (0)
432  pf = sqrt(1.0_dp/(4.0_dp*pi))
433  IF (m /= 0) cpabort("l = 0 and m value out of bounds")
434  y(:) = pf
435  CASE (1)
436  SELECT CASE (m)
437  CASE DEFAULT
438  cpabort("l = 1 and m value out of bounds")
439  CASE (1)
440  pf = sqrt(3.0_dp/(4.0_dp*pi))
441  y(:) = pf*r(1, :)
442  CASE (0)
443  pf = sqrt(3.0_dp/(4.0_dp*pi))
444  y(:) = pf*r(3, :)
445  CASE (-1)
446  pf = sqrt(3.0_dp/(4.0_dp*pi))
447  y(:) = pf*r(2, :)
448  END SELECT
449  CASE (2)
450  SELECT CASE (m)
451  CASE DEFAULT
452  cpabort("l = 2 and m value out of bounds")
453  CASE (2)
454  pf = sqrt(15.0_dp/(16.0_dp*pi))
455  y(:) = pf*(r(1, :)*r(1, :) - r(2, :)*r(2, :))
456  CASE (1)
457  pf = sqrt(15.0_dp/(4.0_dp*pi))
458  y(:) = pf*r(3, :)*r(1, :)
459  CASE (0)
460  pf = sqrt(5.0_dp/(16.0_dp*pi))
461  y(:) = pf*(3.0_dp*r(3, :)*r(3, :) - 1.0_dp)
462  CASE (-1)
463  pf = sqrt(15.0_dp/(4.0_dp*pi))
464  y(:) = pf*r(3, :)*r(2, :)
465  CASE (-2)
466  pf = sqrt(15.0_dp/(16.0_dp*pi))
467  y(:) = pf*2.0_dp*r(1, :)*r(2, :)
468  END SELECT
469  CASE (3)
470  SELECT CASE (m)
471  CASE DEFAULT
472  cpabort("l = 3 and m value out of bounds")
473  CASE (3)
474  pf = sqrt(35.0_dp/(32.0_dp*pi))
475  y(:) = pf*r(1, :)*(r(1, :)**2 - 3.0_dp*r(2, :)**2)
476  CASE (2)
477  pf = sqrt(105.0_dp/(16.0_dp*pi))
478  y(:) = pf*r(3, :)*(r(1, :)**2 - r(2, :)**2)
479  CASE (1)
480  pf = sqrt(21.0_dp/(32.0_dp*pi))
481  y(:) = pf*r(1, :)*(5.0_dp*r(3, :)*r(3, :) - 1.0_dp)
482  CASE (0)
483  pf = sqrt(7.0_dp/(16.0_dp*pi))
484  y(:) = pf*r(3, :)*(5.0_dp*r(3, :)*r(3, :) - 3.0_dp)
485  CASE (-1)
486  pf = sqrt(21.0_dp/(32.0_dp*pi))
487  y(:) = pf*r(2, :)*(5.0_dp*r(3, :)*r(3, :) - 1.0_dp)
488  CASE (-2)
489  pf = sqrt(105.0_dp/(16.0_dp*pi))
490  y(:) = pf*2.0_dp*r(1, :)*r(2, :)*r(3, :)
491  CASE (-3)
492  pf = sqrt(35.0_dp/(32.0_dp*pi))
493  y(:) = pf*r(2, :)*(3.0_dp*r(1, :)**2 - r(2, :)**2)
494  END SELECT
495  CASE DEFAULT
496  IF (m < -l .OR. m > l) cpabort("m value out of bounds")
497  lpm = fac(l + abs(m))
498  lmm = fac(l - abs(m))
499  IF (m == 0) THEN
500  t = 4.0_dp*pi
501  ELSE
502  t = 2.0_dp*pi
503  END IF
504  IF (abs(lpm) < epsilon(1.0_dp)) THEN
505  pf = real(2*l + 1, kind=dp)/t
506  ELSE
507  pf = (real(2*l + 1, kind=dp)*lmm)/(t*lpm)
508  END IF
509  pf = sqrt(pf)
510  DO i = 1, SIZE(r, 2)
511  z = r(3, i)
512 ! plm = legendre ( z, l, m )
513  plm = legendre(z, l, abs(m))
514  IF (m == 0) THEN
515  y(i) = pf*plm
516  ELSE
517  rxy = sqrt(r(1, i)**2 + r(2, i)**2)
518  IF (rxy < epsilon(1.0_dp)) THEN
519  y(i) = 0.0_dp
520  ELSE
521  cp = r(1, i)/rxy
522  sp = r(2, i)/rxy
523  IF (m > 0) THEN
524  y(i) = pf*plm*cosn(m, cp, sp)
525  ELSE
526  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  END SUBROUTINE rvy_lm
534 
535 ! **************************************************************************************************
536 !> \brief ...
537 !> \param r ...
538 !> \param y ...
539 !> \param l ...
540 !> \param m ...
541 ! **************************************************************************************************
542  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  SELECT CASE (l)
559  CASE (:-1)
560  cpabort("Negative l value")
561  CASE (0)
562  pf = sqrt(1.0_dp/(4.0_dp*pi))
563  IF (m /= 0) cpabort("l = 0 and m value out of bounds")
564  y = pf
565  CASE (1)
566  SELECT CASE (m)
567  CASE DEFAULT
568  cpabort("l = 1 and m value out of bounds")
569  CASE (1)
570  pf = sqrt(3.0_dp/(4.0_dp*pi))
571  y = pf*r(1)
572  CASE (0)
573  pf = sqrt(3.0_dp/(4.0_dp*pi))
574  y = pf*r(3)
575  CASE (-1)
576  pf = sqrt(3.0_dp/(4.0_dp*pi))
577  y = pf*r(2)
578  END SELECT
579  CASE (2)
580  SELECT CASE (m)
581  CASE DEFAULT
582  cpabort("l = 2 and m value out of bounds")
583  CASE (2)
584  pf = sqrt(15.0_dp/(16.0_dp*pi))
585  y = pf*(r(1)*r(1) - r(2)*r(2))
586  CASE (1)
587  pf = sqrt(15.0_dp/(4.0_dp*pi))
588  y = pf*r(3)*r(1)
589  CASE (0)
590  pf = sqrt(5.0_dp/(16.0_dp*pi))
591  y = pf*(3.0_dp*r(3)*r(3) - 1.0_dp)
592  CASE (-1)
593  pf = sqrt(15.0_dp/(4.0_dp*pi))
594  y = pf*r(3)*r(2)
595  CASE (-2)
596  pf = sqrt(15.0_dp/(16.0_dp*pi))
597  y = pf*2.0_dp*r(1)*r(2)
598  END SELECT
599  CASE (3)
600  SELECT CASE (m)
601  CASE DEFAULT
602  cpabort("l = 3 and m value out of bounds")
603  CASE (3)
604  pf = sqrt(35.0_dp/(32.0_dp*pi))
605  y = pf*r(1)*(r(1)**2 - 3.0_dp*r(2)**2)
606  CASE (2)
607  pf = sqrt(105.0_dp/(16.0_dp*pi))
608  y = pf*r(3)*(r(1)**2 - r(2)**2)
609  CASE (1)
610  pf = sqrt(21.0_dp/(32.0_dp*pi))
611  y = pf*r(1)*(5.0_dp*r(3)*r(3) - 1.0_dp)
612  CASE (0)
613  pf = sqrt(7.0_dp/(16.0_dp*pi))
614  y = pf*r(3)*(5.0_dp*r(3)*r(3) - 3.0_dp)
615  CASE (-1)
616  pf = sqrt(21.0_dp/(32.0_dp*pi))
617  y = pf*r(2)*(5.0_dp*r(3)*r(3) - 1.0_dp)
618  CASE (-2)
619  pf = sqrt(105.0_dp/(16.0_dp*pi))
620  y = pf*2.0_dp*r(1)*r(2)*r(3)
621  CASE (-3)
622  pf = sqrt(35.0_dp/(32.0_dp*pi))
623  y = pf*r(2)*(3.0_dp*r(1)**2 - r(2)**2)
624  END SELECT
625  CASE DEFAULT
626  IF (m < -l .OR. m > l) cpabort("m value out of bounds")
627  lpm = fac(l + abs(m))
628  lmm = fac(l - abs(m))
629  IF (m == 0) THEN
630  t = 4.0_dp*pi
631  ELSE
632  t = 2.0_dp*pi
633  END IF
634  IF (abs(lpm) < epsilon(1.0_dp)) THEN
635  pf = real(2*l + 1, kind=dp)/t
636  ELSE
637  pf = (real(2*l + 1, kind=dp)*lmm)/(t*lpm)
638  END IF
639  pf = sqrt(pf)
640  z = r(3)
641  plm = legendre(z, l, m)
642  IF (m == 0) THEN
643  y = pf*plm
644  ELSE
645  rxy = sqrt(r(1)**2 + r(2)**2)
646  IF (rxy < epsilon(1.0_dp)) THEN
647  y = 0.0_dp
648  ELSE
649  cp = r(1)/rxy
650  sp = r(2)/rxy
651  IF (m > 0) THEN
652  y = pf*plm*cosn(m, cp, sp)
653  ELSE
654  y = pf*plm*sinn(-m, cp, sp)
655  END IF
656  END IF
657  END IF
658  END SELECT
659 
660  END SUBROUTINE rry_lm
661 
662 ! **************************************************************************************************
663 !> \brief ...
664 !> \param r ...
665 !> \param y ...
666 !> \param l ...
667 !> \param m ...
668 ! **************************************************************************************************
669  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  n = SIZE(r, 2)
688  SELECT CASE (l)
689  CASE (:-1)
690  cpabort("Negative l value")
691  CASE (0)
692  pf = sqrt(1.0_dp/(4.0_dp*pi))
693  IF (m /= 0) cpabort("l = 0 and m value out of bounds")
694  y(:) = pf
695  CASE (1)
696  SELECT CASE (m)
697  CASE DEFAULT
698  cpabort("l = 1 and m value out of bounds")
699  CASE (1)
700  pf = sqrt(3.0_dp/(8.0_dp*pi))
701  DO i = 1, n
702  yp = r(1, i)
703  ym = r(2, i)
704  y(i) = pf*cmplx(yp, ym, kind=dp)
705  END DO
706  CASE (0)
707  pf = sqrt(3.0_dp/(4.0_dp*pi))
708  y(:) = pf*r(3, :)
709  CASE (-1)
710  pf = sqrt(3.0_dp/(8.0_dp*pi))
711  DO i = 1, n
712  yp = r(1, i)
713  ym = r(2, i)
714  y(i) = pf*cmplx(yp, -ym, kind=dp)
715  END DO
716  END SELECT
717  CASE (2)
718  SELECT CASE (m)
719  CASE DEFAULT
720  cpabort("l = 2 and m value out of bounds")
721  CASE (2)
722  pf = sqrt(15.0_dp/(32.0_dp*pi))
723  DO i = 1, n
724  yp = (r(1, i)*r(1, i) - r(2, i)*r(2, i))
725  ym = 2.0_dp*r(1, i)*r(2, i)
726  y(i) = pf*cmplx(yp, ym, kind=dp)
727  END DO
728  CASE (1)
729  pf = sqrt(15.0_dp/(8.0_dp*pi))
730  DO i = 1, n
731  yp = r(3, i)*r(1, i)
732  ym = r(3, i)*r(2, i)
733  y(i) = pf*cmplx(yp, ym, kind=dp)
734  END DO
735  CASE (0)
736  pf = sqrt(5.0_dp/(16.0_dp*pi))
737  y(:) = pf*(3.0_dp*r(3, :)*r(3, :) - 1.0_dp)
738  CASE (-1)
739  pf = sqrt(15.0_dp/(8.0_dp*pi))
740  DO i = 1, n
741  yp = r(3, i)*r(1, i)
742  ym = r(3, i)*r(2, i)
743  y(i) = pf*cmplx(yp, -ym, kind=dp)
744  END DO
745  CASE (-2)
746  pf = sqrt(15.0_dp/(32.0_dp*pi))
747  DO i = 1, n
748  yp = (r(1, i)*r(1, i) - r(2, i)*r(2, i))
749  ym = 2.0_dp*r(1, i)*r(2, i)
750  y(i) = pf*cmplx(yp, -ym, kind=dp)
751  END DO
752  END SELECT
753  CASE (3)
754  SELECT CASE (m)
755  CASE DEFAULT
756  cpabort("l = 3 and m value out of bounds")
757  CASE (3)
758  pf = sqrt(35.0_dp/(64.0_dp*pi))
759  DO i = 1, n
760  yp = r(1, i)*(r(1, i)**2 - 3.0_dp*r(2, i)**2)
761  ym = r(2, i)*(3.0_dp*r(1, i)**2 - r(2, i)**2)
762  y(i) = pf*cmplx(yp, ym, kind=dp)
763  END DO
764  CASE (2)
765  pf = sqrt(105.0_dp/(32.0_dp*pi))
766  DO i = 1, n
767  yp = r(3, i)*(r(1, i)**2 - r(2, i)**2)
768  ym = 2.0_dp*r(1, i)*r(2, i)*r(3, i)
769  y(i) = pf*cmplx(yp, ym, kind=dp)
770  END DO
771  CASE (1)
772  pf = sqrt(21.0_dp/(64.0_dp*pi))
773  DO i = 1, n
774  yp = r(1, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
775  ym = r(2, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
776  y(i) = pf*cmplx(yp, ym, kind=dp)
777  END DO
778  CASE (0)
779  pf = sqrt(7.0_dp/(16.0_dp*pi))
780  y(:) = pf*r(3, :)*(5.0_dp*r(3, :)*r(3, :) - 3.0_dp)
781  CASE (-1)
782  pf = sqrt(21.0_dp/(64.0_dp*pi))
783  DO i = 1, n
784  yp = r(1, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
785  ym = r(2, i)*(5.0_dp*r(3, i)*r(3, i) - 1.0_dp)
786  y(i) = pf*cmplx(yp, -ym, kind=dp)
787  END DO
788  CASE (-2)
789  pf = sqrt(105.0_dp/(32.0_dp*pi))
790  DO i = 1, n
791  yp = r(3, i)*(r(1, i)**2 - r(2, i)**2)
792  ym = 2.0_dp*r(1, i)*r(2, i)*r(3, i)
793  y(i) = pf*cmplx(yp, -ym, kind=dp)
794  END DO
795  CASE (-3)
796  pf = sqrt(35.0_dp/(64.0_dp*pi))
797  DO i = 1, n
798  yp = r(1, i)*(r(1, i)**2 - 3.0_dp*r(2, i)**2)
799  ym = r(2, i)*(3.0_dp*r(1, i)**2 - r(2, i)**2)
800  y(i) = pf*cmplx(yp, -ym, kind=dp)
801  END DO
802  END SELECT
803  CASE DEFAULT
804  IF (m < -l .OR. m > l) cpabort("m value out of bounds")
805  lpm = fac(l + abs(m))
806  lmm = fac(l - abs(m))
807  t = 4.0_dp*pi
808  IF (abs(lpm) < epsilon(1.0_dp)) THEN
809  pf = real(2*l + 1, kind=dp)/t
810  ELSE
811  pf = (real(2*l + 1, kind=dp)*lmm)/(t*lpm)
812  END IF
813  pf = sqrt(pf)
814  DO i = 1, n
815  z = r(3, i)
816  plm = legendre(z, l, m)
817  IF (m == 0) THEN
818  y(i) = pf*plm
819  ELSE
820  rxy = sqrt(r(1, i)**2 + r(2, i)**2)
821  IF (rxy < epsilon(1.0_dp)) THEN
822  y(i) = 0.0_dp
823  ELSE
824  cp = r(1, i)/rxy
825  sp = r(2, i)/rxy
826  IF (m > 0) THEN
827  yp = cosn(m, cp, sp)
828  ym = sinn(m, cp, sp)
829  y(i) = pf*plm*cmplx(yp, ym, kind=dp)
830  ELSE
831  yp = cosn(-m, cp, sp)
832  ym = sinn(-m, cp, sp)
833  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  END SUBROUTINE cvy_lm
841 
842 ! **************************************************************************************************
843 !> \brief ...
844 !> \param r ...
845 !> \param y ...
846 !> \param l ...
847 !> \param m ...
848 ! **************************************************************************************************
849  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  SELECT CASE (l)
867  CASE (:-1)
868  cpabort("Negative l value")
869  CASE (0)
870  pf = sqrt(1.0_dp/(4.0_dp*pi))
871  IF (m /= 0) cpabort("l = 0 and m value out of bounds")
872  y = pf
873  CASE (1)
874  SELECT CASE (m)
875  CASE DEFAULT
876  cpabort("l = 1 and m value out of bounds")
877  CASE (1)
878  pf = sqrt(3.0_dp/(8.0_dp*pi))
879  yp = r(1)
880  ym = r(2)
881  y = pf*cmplx(yp, ym, kind=dp)
882  CASE (0)
883  pf = sqrt(3.0_dp/(4.0_dp*pi))
884  y = pf*r(3)
885  CASE (-1)
886  pf = sqrt(3.0_dp/(8.0_dp*pi))
887  yp = r(1)
888  ym = r(2)
889  y = pf*cmplx(yp, -ym, kind=dp)
890  END SELECT
891  CASE (2)
892  SELECT CASE (m)
893  CASE DEFAULT
894  cpabort("l = 2 and m value out of bounds")
895  CASE (2)
896  pf = sqrt(15.0_dp/(32.0_dp*pi))
897  yp = (r(1)*r(1) - r(2)*r(2))
898  ym = 2.0_dp*r(1)*r(2)
899  y = pf*cmplx(yp, ym, kind=dp)
900  CASE (1)
901  pf = sqrt(15.0_dp/(8.0_dp*pi))
902  yp = r(3)*r(1)
903  ym = r(3)*r(2)
904  y = pf*cmplx(yp, ym, kind=dp)
905  CASE (0)
906  pf = sqrt(5.0_dp/(16.0_dp*pi))
907  y = pf*(3.0_dp*r(3)*r(3) - 1.0_dp)
908  CASE (-1)
909  pf = sqrt(15.0_dp/(8.0_dp*pi))
910  yp = r(3)*r(1)
911  ym = r(3)*r(2)
912  y = pf*cmplx(yp, -ym, kind=dp)
913  CASE (-2)
914  pf = sqrt(15.0_dp/(32.0_dp*pi))
915  yp = (r(1)*r(1) - r(2)*r(2))
916  ym = 2.0_dp*r(1)*r(2)
917  y = pf*cmplx(yp, -ym, kind=dp)
918  END SELECT
919  CASE (3)
920  SELECT CASE (m)
921  CASE DEFAULT
922  cpabort("l = 3 and m value out of bounds")
923  CASE (3)
924  pf = sqrt(35.0_dp/(64.0_dp*pi))
925  yp = r(1)*(r(1)**2 - 3.0_dp*r(2)**2)
926  ym = r(2)*(3.0_dp*r(1)**2 - r(2)**2)
927  y = pf*cmplx(yp, ym, kind=dp)
928  CASE (2)
929  pf = sqrt(105.0_dp/(32.0_dp*pi))
930  yp = r(3)*(r(1)**2 - r(2)**2)
931  ym = 2.0_dp*r(1)*r(2)*r(3)
932  y = pf*cmplx(yp, ym, kind=dp)
933  CASE (1)
934  pf = sqrt(21.0_dp/(64.0_dp*pi))
935  yp = r(1)*(5.0_dp*r(3)*r(3) - 1.0_dp)
936  ym = r(2)*(5.0_dp*r(3)*r(3) - 1.0_dp)
937  y = pf*cmplx(yp, ym, kind=dp)
938  CASE (0)
939  pf = sqrt(7.0_dp/(16.0_dp*pi))
940  y = pf*r(3)*(5.0_dp*r(3)*r(3) - 3.0_dp)
941  CASE (-1)
942  pf = sqrt(21.0_dp/(64.0_dp*pi))
943  yp = r(1)*(5.0_dp*r(3)*r(3) - 1.0_dp)
944  ym = r(2)*(5.0_dp*r(3)*r(3) - 1.0_dp)
945  y = pf*cmplx(yp, -ym, kind=dp)
946  CASE (-2)
947  pf = sqrt(105.0_dp/(32.0_dp*pi))
948  yp = r(3)*(r(1)**2 - r(2)**2)
949  ym = 2.0_dp*r(1)*r(2)*r(3)
950  y = pf*cmplx(yp, -ym, kind=dp)
951  CASE (-3)
952  pf = sqrt(35.0_dp/(64.0_dp*pi))
953  yp = r(1)*(r(1)**2 - 3.0_dp*r(2)**2)
954  ym = r(2)*(3.0_dp*r(1)**2 - r(2)**2)
955  y = pf*cmplx(yp, -ym, kind=dp)
956  END SELECT
957  CASE DEFAULT
958  IF (m < -l .OR. m > l) cpabort("m value out of bounds")
959  lpm = fac(l + abs(m))
960  lmm = fac(l - abs(m))
961  t = 4.0_dp*pi
962  IF (abs(lpm) < epsilon(1.0_dp)) THEN
963  pf = real(2*l + 1, kind=dp)/t
964  ELSE
965  pf = (real(2*l + 1, kind=dp)*lmm)/(t*lpm)
966  END IF
967  pf = sqrt(pf)
968  z = r(3)
969  plm = legendre(z, l, m)
970  IF (m == 0) THEN
971  y = pf*plm
972  ELSE
973  rxy = sqrt(r(1)**2 + r(2)**2)
974  IF (rxy < epsilon(1.0_dp)) THEN
975  y = 0.0_dp
976  ELSE
977  cp = r(1)/rxy
978  sp = r(2)/rxy
979  IF (m > 0) THEN
980  yp = cosn(m, cp, sp)
981  ym = sinn(m, cp, sp)
982  y = pf*plm*cmplx(yp, ym, kind=dp)
983  ELSE
984  yp = cosn(-m, cp, sp)
985  ym = sinn(-m, cp, sp)
986  y = pf*plm*cmplx(yp, -ym, kind=dp)
987  END IF
988  END IF
989  END IF
990  END SELECT
991 
992  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  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  t = c(1)
1028  ct = cos(t)
1029  st = sin(t)
1030  p = c(2)
1031  cp = cos(p)
1032  sp = sin(p)
1033  r(1) = st*cp
1034  r(2) = st*sp
1035  r(3) = ct
1036 
1037 ! dY/dp
1038  IF (m == 0) THEN
1039  dy(2) = 0.0_dp
1040  ELSE
1041  CALL y_lm(r, y, l, -m)
1042  dy(2) = -real(m, kind=dp)*y
1043  END IF
1044 
1045 ! dY/dt
1046  SELECT CASE (l)
1047  CASE (:-1)
1048  cpabort("Negative l value")
1049  CASE (0)
1050  IF (m /= 0) cpabort("l = 0 and m value out of bounds")
1051  dy(1) = 0.0_dp
1052  CASE (1)
1053  SELECT CASE (m)
1054  CASE DEFAULT
1055  cpabort("l = 1 and m value out of bounds")
1056  CASE (1)
1057  pf = sqrt(3.0_dp/(4.0_dp*pi))
1058  dy(1) = pf*ct*cp
1059  CASE (0)
1060  pf = sqrt(3.0_dp/(4.0_dp*pi))
1061  dy(1) = -pf*st
1062  CASE (-1)
1063  pf = sqrt(3.0_dp/(4.0_dp*pi))
1064  dy(1) = pf*ct*sp
1065  END SELECT
1066  CASE (2)
1067  SELECT CASE (m)
1068  CASE DEFAULT
1069  cpabort("l = 2 and m value out of bounds")
1070  CASE (2)
1071  pf = sqrt(15.0_dp/(16.0_dp*pi))
1072  dy(1) = pf*2.0_dp*st*ct*cos(2._dp*p)
1073  CASE (1)
1074  pf = sqrt(15.0_dp/(4.0_dp*pi))
1075  dy(1) = pf*cp*(ct*ct - st*st)
1076  CASE (0)
1077  pf = sqrt(5.0_dp/(16.0_dp*pi))
1078  dy(1) = -pf*6.0_dp*ct*st
1079  CASE (-1)
1080  pf = sqrt(15.0_dp/(4.0_dp*pi))
1081  dy(1) = pf*sp*(ct*ct - st*st)
1082  CASE (-2)
1083  pf = sqrt(15.0_dp/(16.0_dp*pi))
1084  dy(1) = pf*2.0_dp*st*ct*sin(2._dp*p)
1085  END SELECT
1086  CASE (3)
1087  SELECT CASE (m)
1088  CASE DEFAULT
1089  cpabort("l = 3 and m value out of bounds")
1090  CASE (3)
1091  pf = sqrt(35.0_dp/(32.0_dp*pi))
1092  dy(1) = pf*3.0_dp*cos(3._dp*p)*ct*st*st
1093  CASE (2)
1094  pf = sqrt(105.0_dp/(16.0_dp*pi))
1095  dy(1) = pf*2.0_dp*cos(2._dp*p)*ct*st
1096  CASE (1)
1097  pf = sqrt(21.0_dp/(32.0_dp*pi))
1098  dy(1) = pf*cp*(ct*(5.0_dp*ct - 1.0_dp) - 5.0_dp*st*st)
1099  CASE (0)
1100  pf = sqrt(7.0_dp/(16.0_dp*pi))
1101  dy(1) = pf*r(3)*(3.0_dp - 15.0_dp*ct*ct)*st
1102  CASE (-1)
1103  pf = sqrt(21.0_dp/(32.0_dp*pi))
1104  dy(1) = pf*sp*(ct*(5.0_dp*ct - 1.0_dp) - 5.0_dp*st*st)
1105  CASE (-2)
1106  pf = sqrt(105.0_dp/(16.0_dp*pi))
1107  dy(1) = pf*2.0_dp*sin(2._dp*p)*ct*st
1108  CASE (-3)
1109  pf = sqrt(35.0_dp/(32.0_dp*pi))
1110  dy(1) = pf*3.0_dp*sin(3._dp*p)*ct*st*st
1111  END SELECT
1112  CASE DEFAULT
1113  IF (m < -l .OR. m > l) cpabort("m value out of bounds")
1114  lpm = fac(l + abs(m))
1115  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  pf = (real(2*l + 1, kind=dp)*lmm)/(tt*lpm)
1125  END IF
1126  pf = sqrt(pf)
1127  z = ct
1128  dplm = dlegendre(z, l, m)
1129  IF (m == 0) THEN
1130  y = pf*dplm
1131  ELSE
1132  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  END SUBROUTINE dry_lm
1146 
1147 ! **************************************************************************************************
1148 !> \brief ...
1149 !> \param c ...
1150 !> \param dy ...
1151 !> \param l ...
1152 !> \param m ...
1153 ! **************************************************************************************************
1154  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  cpabort("Not implemented")
1179  CALL dry_lm(c, dd, l, m)
1180  dy(1) = cmplx(dd(1), 0.0_dp, kind=dp)
1181  dy(2) = cmplx(dd(2), 0.0_dp, kind=dp)
1182 
1183  END SUBROUTINE dcy_lm
1184 
1185 ! **************************************************************************************************
1186 !> \brief ...
1187 !> \param x ...
1188 !> \param l ...
1189 !> \param m ...
1190 !> \return ...
1191 ! **************************************************************************************************
1192  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  IF (abs(x) > 1.0_dp) cpabort("x value > 1")
1202  SELECT CASE (l)
1203  CASE (:-1)
1204  cpabort("Negative l value")
1205  CASE (0)
1206  plm = 1.0_dp
1207  CASE (1)
1208  SELECT CASE (abs(m))
1209  CASE DEFAULT
1210  cpabort("l = 1 and m value out of bounds")
1211  CASE (1)
1212  plm = sqrt(1.0_dp - x*x)
1213  CASE (0)
1214  plm = x
1215  END SELECT
1216  CASE (2)
1217  SELECT CASE (abs(m))
1218  CASE DEFAULT
1219  cpabort("l = 2 and m value out of bounds")
1220  CASE (2)
1221  plm = 3.0_dp*(1.0_dp - x*x)
1222  CASE (1)
1223  plm = 3.0_dp*x*sqrt(1.0_dp - x*x)
1224  CASE (0)
1225  plm = 1.5_dp*x*x - 0.5_dp
1226  END SELECT
1227  CASE (3)
1228  SELECT CASE (abs(m))
1229  CASE DEFAULT
1230  cpabort("l = 3 and m value out of bounds")
1231  CASE (3)
1232  plm = 15.0_dp*(1.0_dp - x*x)**1.5_dp
1233  CASE (2)
1234  plm = 15.0_dp*x*(1.0_dp - x*x)
1235  CASE (1)
1236  plm = (7.5_dp*x*x - 1.5_dp)*sqrt(1.0_dp - x*x)
1237  CASE (0)
1238  plm = 2.5_dp*x**3 - 1.5_dp*x
1239  END SELECT
1240  CASE (4)
1241  SELECT CASE (abs(m))
1242  CASE DEFAULT
1243  cpabort("l = 4 and m value out of bounds")
1244  CASE (4)
1245  plm = 105.0_dp*(1.0_dp - x*x)**2
1246  CASE (3)
1247  plm = 105.0_dp*x*(1.0_dp - x*x)**1.5_dp
1248  CASE (2)
1249  plm = (52.5_dp*x*x - 7.5_dp)*(1.0_dp - x*x)
1250  CASE (1)
1251  plm = (17.5_dp*x**3 - 7.5_dp*x)*sqrt(1.0_dp - x*x)
1252  CASE (0)
1253  plm = 4.375_dp*x**4 - 3.75_dp*x**2 + 0.375_dp
1254  END SELECT
1255  CASE (5)
1256  SELECT CASE (abs(m))
1257  CASE DEFAULT
1258  cpabort("l = 5 and m value out of bounds")
1259  CASE (5)
1260  plm = 945.0_dp*(1.0_dp - x*x)**2.5_dp
1261  CASE (4)
1262  plm = 945.0_dp*x*(1.0_dp - x*x)**2
1263  CASE (3)
1264  plm = -(-472.5_dp*x*x + 52.5_dp)*(1.0_dp - x*x)**1.5_dp
1265  CASE (2)
1266  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  1.875_dp)*sqrt(1.0_dp - x*x)
1270  CASE (0)
1271  plm = 7.875_dp*x**5 - 8.75_dp*x**3 + 1.875_dp*x
1272  END SELECT
1273  CASE (6)
1274  SELECT CASE (abs(m))
1275  CASE DEFAULT
1276  cpabort("l = 6 and m value out of bounds")
1277  CASE (6)
1278  plm = 10395.0_dp*(1.0_dp - x*x)**3
1279  CASE (5)
1280  plm = 10395.0_dp*x*(1.0_dp - x*x)**2.5_dp
1281  CASE (4)
1282  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  (1.0_dp - x*x)**1.5_dp
1286  CASE (2)
1287  plm = (433.125_dp*x**4 - 236.25_dp*x**2 + &
1288  13.125_dp)*(1.0_dp - x*x)
1289  CASE (1)
1290  plm = -(-86.625_dp*x**5 + 78.75_dp*x**3 - &
1291  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  6.5625_dp*x**2 - 0.3125_dp
1295  END SELECT
1296  CASE DEFAULT
1297  mm = abs(m)
1298  IF (mm > l) cpabort("m out of bounds")
1299 ! use recurence from numerical recipies
1300  pmm = 1.0_dp
1301  IF (mm > 0) THEN
1302  somx2 = sqrt((1.0_dp - x)*(1.0_dp + x))
1303  fact = 1.0_dp
1304  DO im = 1, mm
1305  pmm = pmm*fact*somx2
1306  fact = fact + 2.0_dp
1307  END DO
1308  END IF
1309  IF (l == mm) THEN
1310  plm = pmm
1311  ELSE
1312  pmmp1 = x*real(2*mm + 1, kind=dp)*pmm
1313  IF (l == mm + 1) THEN
1314  plm = pmmp1
1315  ELSE
1316  DO il = mm + 2, l
1317  pll = (x*real(2*il - 1, kind=dp)*pmmp1 - &
1318  REAL(il + mm - 1, kind=dp)*pmm)/real(il - mm, kind=dp)
1319  pmm = pmmp1
1320  pmmp1 = pll
1321  END DO
1322  plm = pll
1323  END IF
1324  END IF
1325  END SELECT
1326 
1327  END FUNCTION legendre
1328 
1329 ! **************************************************************************************************
1330 !> \brief ...
1331 !> \param x ...
1332 !> \param l ...
1333 !> \param m ...
1334 !> \return ...
1335 ! **************************************************************************************************
1336  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  IF (abs(x) > 1.0_dp) cpabort("x value > 1")
1344  SELECT CASE (l)
1345  CASE (0)
1346  dplm = 0.0_dp
1347  CASE (1)
1348  SELECT CASE (abs(m))
1349  CASE DEFAULT
1350  cpabort("l = 1 and m value out of bounds")
1351  CASE (1)
1352  dplm = -x/sqrt(1.0_dp - x*x)
1353  CASE (0)
1354  dplm = 1.0_dp
1355  END SELECT
1356  CASE (2)
1357  SELECT CASE (abs(m))
1358  CASE DEFAULT
1359  cpabort("l = 2 and m value out of bounds")
1360  CASE (2)
1361  dplm = -6.0_dp*x
1362  CASE (1)
1363  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  dplm = 3.0_dp*x
1366  END SELECT
1367  CASE (3)
1368  SELECT CASE (abs(m))
1369  CASE DEFAULT
1370  cpabort("l = 3 and m value out of bounds")
1371  CASE (3)
1372  dplm = -45.0_dp*sqrt(1.0_dp - x*x)*x
1373  CASE (2)
1374  dplm = 15.0_dp*(1.0_dp - x*x) - 30.0_dp*x*x
1375  CASE (1)
1376  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  dplm = 7.5_dp*x*x - 1.5_dp
1379  END SELECT
1380  CASE (4)
1381  SELECT CASE (abs(m))
1382  CASE DEFAULT
1383  cpabort("l = 4 and m value out of bounds")
1384  CASE (4)
1385  dplm = -420*x*(1 - x*x)
1386  CASE (3)
1387  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  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  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  x*(17.5_dp*x**3 - 7.5_dp*x)/sqrt(1.0_dp - x*x)
1396  END IF
1397  CASE (0)
1398  dplm = 4.375_dp*4.0_dp*x**3 - 3.75_dp*2.0_dp*x
1399  END SELECT
1400  CASE (5)
1401  SELECT CASE (abs(m))
1402  CASE DEFAULT
1403  cpabort("l = 5 and m value out of bounds")
1404  CASE (5)
1405  dplm = -945.0_dp*5.0_dp*x*(1.0_dp - x*x)**1.5_dp
1406  CASE (4)
1407  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  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  (157.5_dp*x**3 - 52.5_dp*x)*(-2.0_dp*x)
1414  CASE (1)
1415  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  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  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  SELECT CASE (abs(m))
1426  CASE DEFAULT
1427  cpabort("l = 6 and m value out of bounds")
1428  CASE (6)
1429  dplm = -10395.0_dp*6.0_dp*x*(1.0_dp - x*x)**2
1430  CASE (5)
1431  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  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  (-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  2.0_dp*x*(433.125_dp*x**4 - 236.25_dp*x**2 + 13.125_dp)
1441  CASE (1)
1442  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  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  6.5625_dp*2.0_dp*x
1451  END SELECT
1452  CASE DEFAULT
1453  mm = abs(m)
1454  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  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  + real(l - m + 1, dp)/(x**2 - 1.0_dp)*legendre(x, l + 1, m)
1462  END IF
1463  END SELECT
1464 
1465  END FUNCTION dlegendre
1466 
1467 ! **************************************************************************************************
1468 !> \brief ...
1469 !> \param x ...
1470 !> \param l ...
1471 !> \return ...
1472 ! **************************************************************************************************
1473  FUNCTION dpof1(x, l)
1474 
1475  REAL(kind=dp), INTENT(IN) :: x
1476  INTEGER, INTENT(IN) :: l
1477  REAL(kind=dp) :: dpof1
1478 
1479  IF (abs(x) - 1.0_dp > epsilon(1.0_dp)) THEN
1480  cpabort("|x| is not 1")
1481  END IF
1482  IF (x > 0.0_dp) THEN
1483  SELECT CASE (l)
1484  CASE (:-1)
1485  cpabort("Negative l value")
1486  CASE (0)
1487  dpof1 = 0.0_dp
1488  CASE (1)
1489  dpof1 = 1.0_dp
1490  CASE (2)
1491  dpof1 = 3.0_dp
1492  CASE (3)
1493  dpof1 = 6.0_dp
1494  CASE (4)
1495  dpof1 = 10.0_dp
1496  CASE (5)
1497  dpof1 = 15.0_dp
1498  CASE (6)
1499  dpof1 = 21.0_dp
1500  CASE (7:)
1501  cpabort("Not implemented")
1502  END SELECT
1503  ELSE
1504  SELECT CASE (l)
1505  CASE (:-1)
1506  cpabort("Negative l value")
1507  CASE (0)
1508  dpof1 = 0.0_dp
1509  CASE (1)
1510  dpof1 = 1.0_dp
1511  CASE (2)
1512  dpof1 = -3.0_dp
1513  CASE (3)
1514  dpof1 = 6.0_dp
1515  CASE (4)
1516  dpof1 = -10.0_dp
1517  CASE (5)
1518  dpof1 = 15.0_dp
1519  CASE (6)
1520  dpof1 = -21.0_dp
1521  CASE (7:)
1522  cpabort("Not implemented")
1523  END SELECT
1524  END IF
1525 
1526  END FUNCTION dpof1
1527 
1528 ! **************************************************************************************************
1529 !> \brief ...
1530 !> \param n ...
1531 !> \param k ...
1532 !> \return ...
1533 ! **************************************************************************************************
1534  FUNCTION choose(n, k)
1535 
1536  INTEGER, INTENT(IN) :: n, k
1537  REAL(kind=dp) :: choose
1538 
1539  IF (n >= k) THEN
1540  choose = real(nint(fac(n)/(fac(k)*fac(n - k))), kind=dp)
1541  ELSE
1542  choose = 0.0_dp
1543  END IF
1544 
1545  END FUNCTION choose
1546 
1547 ! **************************************************************************************************
1548 !> \brief ...
1549 !> \param n ...
1550 !> \param c ...
1551 !> \param s ...
1552 !> \return ...
1553 ! **************************************************************************************************
1554  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  cosn = 0.0_dp
1563  IF (abs(c) < epsilon(1.0_dp) .OR. n == 0) THEN
1564  IF (mod(n, 2) == 0) THEN
1565  cosn = (-1.0_dp)**int(n/2)
1566  ELSE
1567  cosn = 0.0_dp
1568  END IF
1569  ELSE
1570  DO i = n, 0, -2
1571  IF (i == 0) THEN
1572  j = n
1573  IF (j == 0) THEN
1574  cosn = cosn + choose(n, j)
1575  ELSE IF (mod(j/2, 2) == 0) THEN
1576  cosn = cosn + choose(n, j)*s**j
1577  ELSE
1578  cosn = cosn - choose(n, j)*s**j
1579  END IF
1580  ELSE
1581  j = n - i
1582  IF (j == 0) THEN
1583  cosn = cosn + choose(n, j)*c**i
1584  ELSE IF (mod(j/2, 2) == 0) THEN
1585  cosn = cosn + choose(n, j)*c**i*s**j
1586  ELSE
1587  cosn = cosn - choose(n, j)*c**i*s**j
1588  END IF
1589  END IF
1590  END DO
1591  END IF
1592 
1593  END FUNCTION cosn
1594 
1595 ! **************************************************************************************************
1596 !> \brief ...
1597 !> \param n ...
1598 !> \param c ...
1599 !> \param s ...
1600 !> \return ...
1601 ! **************************************************************************************************
1602  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  dcosn_dcp = 0.0_dp
1611 
1612  IF (s < epsilon(1.0_dp)) THEN
1613  dcosn_dcp = 0.0_dp
1614  ELSE
1615  DO i = n, 0, -2
1616  IF (i == 0) THEN
1617  dcosn_dcp = dcosn_dcp
1618  ELSE
1619  j = n - i
1620  IF (mod(j/2, 2) == 0) THEN
1621  dcosn_dcp = dcosn_dcp + choose(n, j)*real(i, dp)*c**(i - 1)*s**j
1622  ELSE
1623  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  END FUNCTION dcosn_dcp
1630 
1631 ! **************************************************************************************************
1632 !> \brief ...
1633 !> \param n ...
1634 !> \param c ...
1635 !> \param s ...
1636 !> \return ...
1637 ! **************************************************************************************************
1638  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  dcosn_dsp = 0.0_dp
1647  IF (c < epsilon(1.0_dp) .OR. s < epsilon(1.0_dp)) THEN
1648  dcosn_dsp = 0.0_dp
1649  ELSE
1650  DO i = n, 0, -2
1651  j = n - i
1652  IF (j == 0) THEN
1653  dcosn_dsp = dcosn_dsp
1654  ELSE
1655  IF (mod(j/2, 2) == 0) THEN
1656  dcosn_dsp = dcosn_dsp + choose(n, j)*real(j, dp)*s**(j - 1)*c**i
1657  ELSE
1658  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  END FUNCTION dcosn_dsp
1665 
1666 ! **************************************************************************************************
1667 !> \brief ...
1668 !> \param n ...
1669 !> \param c ...
1670 !> \param s ...
1671 !> \return ...
1672 ! **************************************************************************************************
1673  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  sinn = 0.0_dp
1682 
1683  IF (abs(s) < epsilon(1.0_dp) .OR. n == 0) THEN
1684  sinn = 0.0_dp
1685  ELSE IF (abs(c) < epsilon(1.0_dp)) THEN
1686  IF (mod(n, 2) == 0) THEN
1687  sinn = 0.0_dp
1688  ELSE
1689  sinn = s*(-1.0_dp)**int((n - 1)/2)
1690  END IF
1691  ELSE
1692  DO i = n - 1, 0, -2
1693  IF (i == 0) THEN
1694  j = n
1695  IF (mod(j/2, 2) == 0) THEN
1696  sinn = sinn + choose(n, j)*s**j
1697  ELSE
1698  sinn = sinn - choose(n, j)*s**j
1699  END IF
1700  ELSE
1701  j = n - i
1702  IF (mod(j/2, 2) == 0) THEN
1703  sinn = sinn + choose(n, j)*c**i*s**j
1704  ELSE
1705  sinn = sinn - choose(n, j)*c**i*s**j
1706  END IF
1707  END IF
1708  END DO
1709  END IF
1710 
1711  END FUNCTION sinn
1712 
1713 ! **************************************************************************************************
1714 !> \brief ...
1715 !> \param n ...
1716 !> \param c ...
1717 !> \param s ...
1718 !> \return ...
1719 ! **************************************************************************************************
1720  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  dsinn_dcp = 0.0_dp
1729 
1730  IF (c < epsilon(1.0_dp) .OR. s < epsilon(1.0_dp)) THEN
1731  dsinn_dcp = 0.0_dp
1732  ELSE
1733  DO i = n - 1, 0, -2
1734  IF (i == 0) THEN
1735  dsinn_dcp = dsinn_dcp
1736  ELSE
1737  j = n - i
1738  IF (mod(j/2, 2) == 0) THEN
1739  dsinn_dcp = dsinn_dcp + choose(n, j)*real(i, dp)*c**(i - 1)*s**j
1740  ELSE
1741  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  END FUNCTION dsinn_dcp
1748 
1749 ! **************************************************************************************************
1750 !> \brief ...
1751 !> \param n ...
1752 !> \param c ...
1753 !> \param s ...
1754 !> \return ...
1755 ! **************************************************************************************************
1756  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  dsinn_dsp = 0.0_dp
1765 
1766  IF (c < epsilon(1.0_dp) .OR. s < epsilon(1.0_dp)) THEN
1767  dsinn_dsp = 0.0_dp
1768  ELSE
1769  DO i = n - 1, 0, -2
1770  j = n - i
1771  IF (j == 0) THEN
1772  dsinn_dsp = dsinn_dsp
1773  ELSE
1774  IF (mod(j/2, 2) == 0) THEN
1775  dsinn_dsp = dsinn_dsp + choose(n, j)*c**i*real(j, dp)*s**(j - 1)
1776  ELSE
1777  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  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  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  IF (j1 < 0.0_dp) &
1809  cpabort("The angular momentum quantum number j1 has to be nonnegative")
1810  IF (.NOT. (is_integer(j1) .OR. is_integer(2.0_dp*j1))) &
1811  cpabort("The angular momentum quantum number j1 has to be integer or half-integer")
1812  IF (j2 < 0.0_dp) &
1813  cpabort("The angular momentum quantum number j2 has to be nonnegative")
1814  IF (.NOT. (is_integer(j2) .OR. is_integer(2.0_dp*j2))) &
1815  cpabort("The angular momentum quantum number j2 has to be integer or half-integer")
1816  IF (j < 0.0_dp) &
1817  cpabort("The angular momentum quantum number J has to be nonnegative")
1818  IF (.NOT. (is_integer(j) .OR. is_integer(2.0_dp*j))) &
1819  cpabort("The angular momentum quantum number J has to be integer or half-integer")
1820  IF ((abs(m1) - j1) > epsilon(m1)) &
1821  cpabort("The angular momentum quantum number m1 has to satisfy -j1 <= m1 <= j1")
1822  IF (.NOT. (is_integer(m1) .OR. is_integer(2.0_dp*m1))) &
1823  cpabort("The angular momentum quantum number m1 has to be integer or half-integer")
1824  IF ((abs(m2) - j2) > epsilon(m2)) &
1825  cpabort("The angular momentum quantum number m2 has to satisfy -j2 <= m1 <= j2")
1826  IF (.NOT. (is_integer(m2) .OR. is_integer(2.0_dp*m2))) &
1827  cpabort("The angular momentum quantum number m2 has to be integer or half-integer")
1828  IF ((abs(m) - j) > epsilon(m)) &
1829  cpabort("The angular momentum quantum number M has to satisfy -J <= M <= J")
1830  IF (.NOT. (is_integer(m) .OR. is_integer(2.0_dp*m))) &
1831  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  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  (abs(m2) > j2) .OR. &
1843  (abs(m) > j)) THEN
1844  cg_coeff = 0.0_dp
1845  ELSE
1846  ! Compute Clebsch-Gordan coefficient
1847  sumk = 0.0_dp
1848  kmax = nint(max(j1 + j2 - j, j1 - j + m2, j2 - j - m1, j1 - m1, j2 + m2))
1849  DO k = 0, kmax
1850  IF (j1 + j2 - j - k < 0.0_dp) cycle
1851  IF (j - j1 - m2 + k < 0.0_dp) cycle
1852  IF (j - j2 + m1 + k < 0.0_dp) cycle
1853  IF (j1 - m1 - k < 0.0_dp) cycle
1854  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  sfac(nint(j2 + m2 - k))*sfac(k)
1858  IF (mod(k, 2) == 1) t = -t
1859  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  sfac(nint(j + m))*sfac(nint(j - m)))*sumk
1866  END IF
1867  ELSE
1868  cpabort("Invalid set of input parameters < j1 m1 j2 m2 | J M > specified")
1869  END IF
1870 
1871  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  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  (abs(j1 - j2) <= j3) .AND. (j3 <= abs(j1 + j2)) .AND. &
1896  is_integer(j1 + j2 + j3)) THEN
1897  IF (is_integer(j1 - j2 - m3)) THEN
1898  CALL clebsch_gordon_coefficient(j1, m1, j2, m2, j3, -m3, cg_coeff)
1899  w_3j = (-1.0_dp)**int(j1 - j2 - m3)*cg_coeff/sqrt(2.0_dp*j3 + 1.0_dp)
1900  ELSE
1901  cpabort("j1 - j2 - m3 results in an invalid non-integer exponent")
1902  END IF
1903  ELSE
1904  w_3j = 0.0_dp
1905  END IF
1906 
1907  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  FUNCTION is_integer(x)
1916 
1917  REAL(kind=dp), INTENT(IN) :: x
1918  LOGICAL :: is_integer
1919 
1920  IF ((abs(x) - aint(abs(x), kind=dp)) > epsilon(x)) THEN
1921  is_integer = .false.
1922  ELSE
1923  is_integer = .true.
1924  END IF
1925 
1926  END FUNCTION is_integer
1927 
1928 END MODULE spherical_harmonics
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
mathconstants::maxfac
integer, parameter, public maxfac
Definition: mathconstants.F:36
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition: mathconstants.F:37
spherical_harmonics
Calculate spherical harmonics.
Definition: spherical_harmonics.F:27
spherical_harmonics::legendre
real(kind=dp) function, public legendre(x, l, m)
...
Definition: spherical_harmonics.F:1193
spherical_harmonics::clebsch_gordon_coefficient
subroutine, public clebsch_gordon_coefficient(j1, m1, j2, m2, J, M, CG_coeff)
Compute the Clebsch-Gordon coefficient C = < j1 m1 j2 m2 | J M > = < j1 j2; m1 m2 | J M >
Definition: spherical_harmonics.F:1800
spherical_harmonics::wigner_3j
subroutine, public wigner_3j(j1, m1, j2, m2, j3, m3, W_3j)
Compute the Wigner 3-j symbol / j1 j2 j3 \ \ m1 m2 m3 / using the Clebsch-Gordon coefficients.
Definition: spherical_harmonics.F:1888
spherical_harmonics::clebsch_gordon_init
subroutine, public clebsch_gordon_init(l)
...
Definition: spherical_harmonics.F:238
spherical_harmonics::clebsch_gordon_deallocate
subroutine, public clebsch_gordon_deallocate()
...
Definition: spherical_harmonics.F:220
spherical_harmonics::dlegendre
real(kind=dp) function, public dlegendre(x, l, m)
...
Definition: spherical_harmonics.F:1337