LORENE
et_bin_bhns_extr_ylm.C
1/*
2 * Method of class Et_bin_bhns_extr to construct spherical harmonics
3 *
4 * (see file et_bin_bhns_extr.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2004 Joshua A. Faber
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28char et_bin_bhns_extr_ylm_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_ylm.C,v 1.5 2014/10/13 08:52:55 j_novak Exp $" ;
29
30/*
31 * $Id: et_bin_bhns_extr_ylm.C,v 1.5 2014/10/13 08:52:55 j_novak Exp $
32 * $Log: et_bin_bhns_extr_ylm.C,v $
33 * Revision 1.5 2014/10/13 08:52:55 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.4 2014/10/06 15:13:08 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.3 2005/02/28 23:18:07 k_taniguchi
40 * Change the functions to constant ones
41 *
42 * Revision 1.2 2005/01/03 19:52:56 k_taniguchi
43 * Change a factor multiplied/divided by sqrt(2).
44 *
45 * Revision 1.1 2004/12/29 16:30:46 k_taniguchi
46 * *** empty log message ***
47 *
48 *
49 * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_ylm.C,v 1.5 2014/10/13 08:52:55 j_novak Exp $
50 *
51 */
52
53// C headers
54#include <cstdlib>
55#include <cmath>
56
57// Lorene headers
58#include "map.h"
59#include "tenseur.h"
60#include "et_bin_bhns_extr.h"
61
62namespace Lorene {
63void Et_bin_bhns_extr::get_ylm(int nylm, Cmp** ylmvec) const {
64
65 // IMPORTANT NOTE:
66 // For Y_lm with m>=1, we have the real and imaginary parts,
67 // not Y_{l,m} and Y_{l,-m}. This changes the normalization
68 // properties. In order to normalize properly, we multiply
69 // all fields in get_integrals below by a factor of 2.0 when
70 // m>=1.
71
72 cout << "Constructing ylm" << endl;
73
74 int nz = mp.get_mg()->get_nzone() ;
75 int nr = mp.get_mg()->get_nr(0) ;
76 int np = mp.get_mg()->get_np(0) ;
77 int nt = mp.get_mg()->get_nt(0) ;
78
79 for (int l=0 ; l<nz ; l++) {
80
81 Mtbl Xabs (mp.x) ;
82 Mtbl Yabs (mp.y) ;
83 Mtbl Zabs (mp.z) ;
84
85 for (int k=0 ; k<np ; k++) {
86 for (int j=0 ; j<nt ; j++) {
87 for (int i=0 ; i<nr ; i++) {
88
89 double xval=Xabs(l,k,j,i);
90 double yval=Yabs(l,k,j,i);
91 double zval=Zabs(l,k,j,i);
92 double rval=sqrt(xval*xval+yval*yval+zval*zval);
93
94 // cout <<l<<" "<<k<<" "<<j<<" "<<i<<endl;
95
96 //l=0,m=0
97 ylmvec[0]->set(l,k,j,i)=1.0*sqrt(1.0/4.0/M_PI);
98 // cout << " 0 " << endl;
99 // l=1 included?
100 if (nylm>1 ) {
101 if (nylm <4) {abort();} else {
102 //l=1,m=0
103 ylmvec[1]->set(l,k,j,i)=zval*sqrt(3.0/4.0/M_PI);
104 //l=1,m=1
105 ylmvec[2]->set(l,k,j,i)=-1.0*xval*sqrt(3.0/8.0/M_PI);
106 ylmvec[3]->set(l,k,j,i)=-1.0*yval*sqrt(3.0/8.0/M_PI);
107 }
108 }
109 // l=2 included?
110 if (nylm>4 ) {
111 if (nylm <9) {abort();} else {
112 //l=2,m=0
113 ylmvec[4]->set(l,k,j,i)=(3.0*zval*zval-rval*rval)*sqrt(5.0/16.0/M_PI);
114 //l=2,m=1
115 ylmvec[5]->set(l,k,j,i)=-1.0*zval*xval*sqrt(15.0/8.0/M_PI);
116 ylmvec[6]->set(l,k,j,i)=-1.0*zval*yval*sqrt(15.0/8.0/M_PI);
117 //l=2,m=2
118 ylmvec[7]->set(l,k,j,i)=(xval*xval-yval*yval)*sqrt(15.0/32.0/M_PI);
119 ylmvec[8]->set(l,k,j,i)=2.0*xval*yval*sqrt(15.0/32.0/M_PI);
120 }
121 }
122 // l=3 included?
123 if (nylm>9 ) {
124 if (nylm <16) {abort();} else {
125 //l=3,m=0
126 ylmvec[9]->set(l,k,j,i)=(5.0*pow(zval,3)-3.0*zval*rval*rval)*
127 sqrt(7.0/16.0/M_PI);
128 //l=3,m=1
129 ylmvec[10]->set(l,k,j,i)=-1.0*(5.0*zval*zval-rval*rval)*xval*
130 sqrt(21.0/64.0/M_PI);
131 ylmvec[11]->set(l,k,j,i)=-1.0*(5.0*zval*zval-rval*rval)*yval*
132 sqrt(21.0/64.0/M_PI);
133 //l=3,m=2
134 ylmvec[12]->set(l,k,j,i)=zval*(xval*xval-yval*yval)*
135 sqrt(105./32.0/M_PI);
136 ylmvec[13]->set(l,k,j,i)=zval*(2.0*xval*yval)*
137 sqrt(105./32.0/M_PI);
138 //l=3,m=3
139 ylmvec[14]->set(l,k,j,i)=-1.0*(pow(xval,3)-3.0*xval*yval*yval)*
140 sqrt(35.0/64.0/M_PI);
141 ylmvec[15]->set(l,k,j,i)=-1.0*(3.0*xval*xval*yval-pow(yval,3))*
142 sqrt(35.0/64.0/M_PI);
143 }
144 }
145 // l=4 included?
146 if (nylm>16 ) {
147 if (nylm <25) {abort();} else {
148 //l=4,m=0
149 ylmvec[16]->set(l,k,j,i)=(35.0*pow(zval,4)-30.0*zval*zval*rval*rval+3*pow(rval,4))*
150 sqrt(9.0/256.0/M_PI);
151 //l=4,m=1
152 ylmvec[17]->set(l,k,j,i)=-1.0*(7.0*pow(zval,3)-3*zval*rval*rval)*xval*
153 sqrt(45.0/64.0/M_PI);
154 ylmvec[18]->set(l,k,j,i)=-1.0*(7.0*pow(zval,3)-3*zval*rval*rval)*yval*
155 sqrt(45.0/64.0/M_PI);
156 //l=4,m=2
157 ylmvec[19]->set(l,k,j,i)=(7.0*zval*zval-rval*rval)*(xval*xval-yval*yval)*
158 sqrt(45./128.0/M_PI);
159 ylmvec[20]->set(l,k,j,i)=(7.0*zval*zval-rval*rval)*(2.0*xval*yval)*
160 sqrt(45./128.0/M_PI);
161 //l=4,m=3
162 ylmvec[21]->set(l,k,j,i)=-1.0*zval*(pow(xval,3)-3.0*xval*yval*yval)*
163 sqrt(315.0/64.0/M_PI);
164 ylmvec[22]->set(l,k,j,i)=-1.0*zval*(3.0*xval*xval*yval-pow(yval,3))*
165 sqrt(315.0/64.0/M_PI);
166 //l=4,m=4
167 ylmvec[23]->set(l,k,j,i)=(pow(xval,4)-6*xval*xval*yval*yval+pow(yval,4))*
168 sqrt(315.0/512.0/M_PI);
169 ylmvec[24]->set(l,k,j,i)=4.0*xval*yval*(xval*xval-yval*yval)*
170 sqrt(315.0/512.0/M_PI);
171 }
172 }
173 // l=5 included?
174 if (nylm>25 ) {
175 if (nylm <36) {abort();} else {
176 //l=5,m=0
177 ylmvec[25]->set(l,k,j,i)=(63.0*pow(zval,5)-70.0*pow(zval,3)*rval*rval+15*zval*pow(rval,4))*
178 sqrt(11.0/256.0/M_PI);
179 //l=5,m=1
180 ylmvec[26]->set(l,k,j,i)=-1.0*(21.0*pow(zval,4)-14*zval*zval*rval*rval+pow(rval,4))*xval*
181 sqrt(165.0/512.0/M_PI);
182 ylmvec[27]->set(l,k,j,i)=-1.0*(21.0*pow(zval,4)-14*zval*zval*rval*rval+pow(rval,4))*yval*
183 sqrt(165.0/512.0/M_PI);
184 //l=5,m=2
185 ylmvec[28]->set(l,k,j,i)=(3.0*pow(zval,3)-zval*rval*rval)*(xval*xval-yval*yval)*
186 sqrt(1155./128.0/M_PI);
187 ylmvec[29]->set(l,k,j,i)=(3.0*pow(zval,3)-zval*rval*rval)*(2.0*xval*yval)*
188 sqrt(1155./128.0/M_PI);
189 //l=5,m=3
190 ylmvec[30]->set(l,k,j,i)=-1.0*(9.0*zval*zval-rval*rval)*(pow(xval,3)-3.0*xval*yval*yval)*
191 sqrt(385.0/1024.0/M_PI);
192 ylmvec[31]->set(l,k,j,i)=-1.0*(9.0*zval*zval-rval*rval)*(3.0*xval*xval*yval-pow(yval,3))*
193 sqrt(385.0/1024.0/M_PI);
194 //l=5,m=4
195 ylmvec[32]->set(l,k,j,i)=zval*(pow(xval,4)-6*xval*xval*yval*yval+pow(yval,4))*
196 sqrt(3465.0/512.0/M_PI);
197 ylmvec[33]->set(l,k,j,i)=zval*4.0*xval*yval*(xval*xval-yval*yval)*
198 sqrt(3465.0/512.0/M_PI);
199 //l=5,m=5
200 ylmvec[34]->set(l,k,j,i)=-1.0*(pow(xval,5)-10.0*pow(xval,3)*yval*yval+5.0*xval*pow(yval,4))*
201 sqrt(693.0/1024.0/M_PI);
202 ylmvec[35]->set(l,k,j,i)=-1.0*(5.0*pow(xval,4)*yval-10.0*xval*xval*pow(yval,3)+pow(yval,5))*
203 sqrt(693.0/1024.0/M_PI);
204 }
205 }
206 // l=6 included?
207 if (nylm>36 ) {
208 if (nylm <49) {abort();} else {
209 //l=6,m=0
210 ylmvec[36]->set(l,k,j,i)=(231.0*pow(zval,6)-315.0*pow(zval,4)*rval*rval+105.0*zval*zval*pow(rval,4)-5.0*pow(rval,6))*
211 sqrt(13.0/1024.0/M_PI);
212 //l=6,m=1
213 ylmvec[37]->set(l,k,j,i)=-1.0*(33.0*pow(zval,5)-30.0*pow(zval,3)*rval*rval+5.0*zval*pow(rval,4))*xval*
214 sqrt(273.0/512.0/M_PI);
215 ylmvec[38]->set(l,k,j,i)=-1.0*(33.0*pow(zval,5)-30.0*pow(zval,3)*rval*rval+5.0*zval*pow(rval,4))*yval*
216 sqrt(273.0/512.0/M_PI);
217 //l=6,m=2
218 ylmvec[39]->set(l,k,j,i)=(33.0*pow(zval,4)-18.0*zval*zval*rval*rval+pow(rval,4))*(xval*xval-yval*yval)*
219 sqrt(1365./4096.0/M_PI);
220 ylmvec[40]->set(l,k,j,i)=(33.0*pow(zval,4)-18.0*zval*zval*rval*rval+pow(rval,4))*(2.0*xval*yval)*
221 sqrt(1365./4096.0/M_PI);
222 //l=6,m=3
223 ylmvec[41]->set(l,k,j,i)=-1.0*(11.0*pow(zval,3)-3.0*zval*rval*rval)*(pow(xval,3)-3.0*xval*yval*yval)*
224 sqrt(1365.0/1024.0/M_PI);
225 ylmvec[42]->set(l,k,j,i)=-1.0*(11.0*pow(zval,3)-3.0*zval*rval*rval)*(3.0*xval*xval*yval-pow(yval,3))*
226 sqrt(1365.0/1024.0/M_PI);
227 //l=6,m=4
228 ylmvec[43]->set(l,k,j,i)=(11.0*zval*zval-rval*rval)*(pow(xval,4)-6*xval*xval*yval*yval+pow(yval,4))*
229 sqrt(819.0/2048.0/M_PI);
230 ylmvec[44]->set(l,k,j,i)=(11.0*zval*zval-rval*rval)*4.0*xval*yval*(xval*xval-yval*yval)*
231 sqrt(819.0/2048.0/M_PI);
232 //l=6,m=5
233 ylmvec[45]->set(l,k,j,i)=-1.0*zval*(pow(xval,5)-10.0*pow(xval,3)*yval*yval+5.0*xval*pow(yval,4))*
234 sqrt(9009.0/1024.0/M_PI);
235 ylmvec[46]->set(l,k,j,i)=-1.0*zval*(5.0*pow(xval,4)*yval-10.0*xval*xval*pow(yval,3)+pow(yval,5))*
236 sqrt(9009.0/1024.0/M_PI);
237 //l=6,m=6
238 ylmvec[47]->set(l,k,j,i)=(pow(xval,6)-15.0*pow(xval,4)*yval*yval+15.0*xval*xval*pow(yval,4)-pow(yval,6))*
239 sqrt(3003.0/4096.0/M_PI);
240 ylmvec[48]->set(l,k,j,i)=(6.0*pow(xval,5)*yval-20.0*pow(xval,3)*pow(yval,3)+6.0*xval*pow(yval,5))*
241 sqrt(3003.0/4096.0/M_PI);
242 }
243 }
244 if(nylm >49) {
245 cout << "l>6 not implemented!!!!!!!"<< endl;
246 abort();
247 }
248 }
249 }
250 }
251 }
252
253}
254
255void Et_bin_bhns_extr::get_integrals(int nylm, double* intvec, Cmp** ylmvec,
256 Cmp source) const {
257
258 // As mentioned in the comment before get_ylm, our real/imaginary
259 // representation of the Y_lm Cmp's does not agree with the normalization
260 // used to define
261 // the spherical harmonic decomposition of Cmp arrays. Thus, we multiply
262 // all terms with m>=1 by a factor of 2.0 in order to
263 // produce the correct decomposition.
264
265 int nz=mp.get_mg()->get_nzone() ;
266
267 Map_af mapping (mp);
268
269 const double* a1 = mapping.get_alpha() ;
270 const double* b1 = mapping.get_beta() ;
271
272 double rlim=a1[nz-1]+b1[nz-1];
273
274 int ll=0;
275 int mm=0;
276 int ncount=0;
277 for (int n=0; n<nylm; n++) {
278
279 Cmp ylmsource=(*ylmvec[n]*source);
280 int symcheck=1;
281 for (int l=0; l<nz; l++) {
282 int symc=ylmsource.va.base.get_base_t(l);
283 if(symc!=2304 && symc!=1280)symcheck=0;
284 }
285 if(symcheck==1) {
286 intvec[n]=ylmsource.integrale()/(2.0*ll+1.0)/sqrt(2.0*M_PI)/pow(rlim,ll+1);
287 } else {
288 intvec[n]=0;
289 }
290 if(mm>=1)intvec[n]*=2.0;
291
292 int lnew=0;
293 int mnew=0;
294 int nnew=0;
295 if(mm<ll) {
296 if(mm==0) {
297 lnew=ll;
298 mnew=mm+1;
299 nnew=0;
300 }
301 if(mm>0&&ncount==0) {
302 lnew=ll;
303 mnew=mm;
304 nnew=1;
305 }
306 if(mm>0&&ncount==1) {
307 lnew=ll;
308 mnew=mm+1;
309 nnew=0;
310 }
311 }
312 if(mm==ll) {
313 if(mm==0) {
314 lnew=ll+1;
315 mnew=0;
316 nnew=0;
317 }
318 if(mm>0&&ncount==0) {
319 lnew=ll;
320 mnew=mm;
321 nnew=1;
322 }
323 if(mm>0&&ncount==1) {
324 lnew=ll+1;
325 mnew=0;
326 nnew=0;
327 }
328 }
329 ll=lnew;
330 mm=mnew;
331 ncount=nnew;
332 }
333}
334}
int get_base_t(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
Definition base_val.h:411
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
Tbl & set(int l)
Read/write of the value in a given domain.
Definition cmp.h:724
double integrale() const
Computes the integral over all space of *this .
Definition cmp_integ.C:55
void get_ylm(int nylm, Cmp **ylmvec) const
Constructs spherical harmonics.
void get_integrals(int nylm, double *intvec, Cmp **ylmvec, Cmp source) const
Computes multipole moments.
Map & mp
Mapping associated with the star.
Definition etoile.h:429
Affine radial mapping.
Definition map.h:2027
const double * get_beta() const
Returns the pointer on the array beta.
Definition map_af.C:481
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition map_af.C:477
Coord y
y coordinate centered on the grid
Definition map.h:727
Coord z
z coordinate centered on the grid
Definition map.h:728
Coord x
x coordinate centered on the grid
Definition map.h:726
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Multi-domain array.
Definition mtbl.h:118
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene prototypes.
Definition app_hor.h:64