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