1 /*
2 Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17 02110-1301, USA.
18
19 */
20
21 #include <config.h>
22 #include <stdlib.h>
23 #include <math.h>
24 #include <float.h>
25 #include <fortran_types.h>
26
27 #include <gsl/gsl_sf_legendre.h>
28
29 /* normalization constants for the spherical harmonics */
30 static double sph_cnsts[16] = {
31 0.282094791773878, /* l = 0 */
32 0.488602511902920, 0.488602511902920, 0.488602511902920, /* l = 1 */
33 /* l = 2 */
34 0.182091405098680, 0.364182810197360, 0.630783130505040, 0.364182810197360, 0.182091405098680,
35 /* l = 3 */
36 0.59004359, 2.890611444, 0.4570458, 0.373176333, 0.4570458, 1.445305722, 0.59004359
37 };
38
39 /* Computes real spherical harmonics ylm in the direction of vector r:
40 ylm = c * plm( cos(theta) ) * sin(m*phi) for m < 0
41 ylm = c * plm( cos(theta) ) * cos(m*phi) for m >= 0
42 with (theta,phi) the polar angles of r, c a positive normalization
43 constant and plm associated legendre polynomials. */
44
FC_FUNC_(oct_ylm,OCT_YLM)45 void FC_FUNC_(oct_ylm, OCT_YLM)
46 (const fint * np, const double *x, const double *y, const double *z, const fint *l, const fint *m, double * restrict ylm)
47 {
48 double r, r2, r3, rr, rx, ry, rz, cosphi, sinphi, cosm, sinm, phase;
49 int i, ip;
50
51 if(l[0] == 0) {
52 for (ip = 0; ip < np[0]; ip++) ylm[ip] = sph_cnsts[0];
53 return;
54 }
55
56 for (ip = 0; ip < np[0]; ip++){
57
58 r2 = x[ip]*x[ip] + y[ip]*y[ip] + z[ip]*z[ip];
59
60 /* if r=0, direction is undefined => make ylm=0 except for l=0 */
61 if(r2 < 1.0e-15){
62 ylm[ip] = 0.0;
63 continue;
64 }
65
66 switch(l[0]){
67 case 1:
68 rr = 1.0/sqrt(r2);
69 switch(m[0]){
70 case -1: ylm[ip] = -sph_cnsts[1]*rr*y[ip]; continue;
71 case 0: ylm[ip] = sph_cnsts[2]*rr*z[ip]; continue;
72 case 1: ylm[ip] = -sph_cnsts[3]*rr*x[ip]; continue;
73 }
74 case 2:
75 switch(m[0]){
76 case -2: ylm[ip] = sph_cnsts[4]*6.0*x[ip]*y[ip]/r2;
77 continue;
78 case -1: ylm[ip] = -sph_cnsts[5]*3.0*y[ip]*z[ip]/r2;
79 continue;
80 case 0: ylm[ip] = sph_cnsts[6]*0.5*(3.0*z[ip]*z[ip]/r2 - 1.0);
81 continue;
82 case 1: ylm[ip] = -sph_cnsts[7]*3.0*x[ip]*z[ip]/r2;
83 continue;
84 case 2: ylm[ip] = sph_cnsts[8]*3.0*(x[ip]*x[ip] - y[ip]*y[ip])/r2;
85 continue;
86 }
87 case 3:
88 r3 = r2 * sqrt(r2);
89 switch(m[0]){
90 case -3: ylm[ip] = sph_cnsts[9] *((3.0*x[ip]*x[ip]-y[ip]*y[ip]) * y[ip])/r3;
91 continue;
92 case -2: ylm[ip] = sph_cnsts[10]*(x[ip]*y[ip]*z[ip])/r3;
93 continue;
94 case -1: ylm[ip] = sph_cnsts[11]*(y[ip]*(5.0*z[ip]*z[ip]-r2))/r3;
95 continue;
96 case 0: ylm[ip] = sph_cnsts[12]*(z[ip]*(5.0*z[ip]*z[ip]-3.0*r2))/r3;
97 continue;
98 case 1: ylm[ip] = sph_cnsts[13]*(x[ip]*(5.0*z[ip]*z[ip]-r2))/r3;
99 continue;
100 case 2: ylm[ip] = sph_cnsts[14]*((x[ip]*x[ip]-y[ip]*y[ip])*z[ip])/r3;
101 continue;
102 case 3: ylm[ip] = sph_cnsts[15]*((x[ip]*x[ip]-3.0*y[ip]*y[ip])*x[ip])/r3;
103 continue;
104 }
105 }
106
107 /* get phase */
108 rr = 1.0/sqrt(r2);
109
110 rx = x[ip]*rr;
111 ry = y[ip]*rr;
112 rz = z[ip]*rr;
113
114 r = hypot(rx, ry);
115 if(r < 1e-20) r = 1e-20; /* one never knows... */
116
117 cosphi = rx/r;
118 sinphi = ry/r;
119
120 /* compute sin(mphi) and cos(mphi) by adding cos/sin */
121 cosm = 1.; sinm=0.;
122 for(i = 0; i < abs(m[0]); i++){
123 double a = cosm, b = sinm;
124 cosm = a*cosphi - b*sinphi;
125 sinm = a*sinphi + b*cosphi;
126 }
127 phase = m[0] < 0 ? sinm : cosm;
128 phase = m[0] == 0 ? phase : sqrt(2.0)*phase;
129
130 /* adding small number (~= 10^-308) to avoid floating invalids */
131 rz = rz + DBL_MIN;
132
133 r = gsl_sf_legendre_sphPlm(l[0], abs(m[0]), rz);
134
135 /* I am not sure whether we are including the Condon-Shortley factor (-1)^m */
136 ylm[ip] = r*phase;
137 }
138 }
139
140