63 #include "et_bin_bhns_extr.h" 75 cout <<
"Constructing ylm" << endl;
82 for (
int l=0 ; l<nz ; l++) {
88 for (
int k=0 ; k<np ; k++) {
89 for (
int j=0 ; j<nt ; j++) {
90 for (
int i=0 ; i<nr ; i++) {
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);
100 ylmvec[0]->
set(l,k,j,i)=1.0*
sqrt(1.0/4.0/M_PI);
104 if (nylm <4) {abort();}
else {
106 ylmvec[1]->
set(l,k,j,i)=zval*
sqrt(3.0/4.0/M_PI);
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);
114 if (nylm <9) {abort();}
else {
116 ylmvec[4]->
set(l,k,j,i)=(3.0*zval*zval-rval*rval)*
sqrt(5.0/16.0/M_PI);
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);
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);
127 if (nylm <16) {abort();}
else {
129 ylmvec[9]->
set(l,k,j,i)=(5.0*
pow(zval,3)-3.0*zval*rval*rval)*
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);
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);
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);
150 if (nylm <25) {abort();}
else {
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);
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);
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);
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);
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);
178 if (nylm <36) {abort();}
else {
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);
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);
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);
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);
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);
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);
211 if (nylm <49) {abort();}
else {
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);
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);
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);
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);
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);
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);
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);
248 cout <<
"l>6 not implemented!!!!!!!"<< endl;
273 const double* b1 = mapping.
get_beta() ;
275 double rlim=a1[nz-1]+b1[nz-1];
280 for (
int n=0; n<nylm; n++) {
282 Cmp ylmsource=(*ylmvec[n]*source);
284 for (
int l=0; l<nz; l++) {
286 if(symc!=2304 && symc!=1280)symcheck=0;
293 if(mm>=1)intvec[n]*=2.0;
304 if(mm>0&&ncount==0) {
309 if(mm>0&&ncount==1) {
321 if(mm>0&&ncount==0) {
326 if(mm>0&&ncount==1) {
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
const double * get_alpha() const
Returns the pointer on the array alpha.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Cmp sqrt(const Cmp &)
Square root.
int get_base_t(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
void get_integrals(int nylm, double *intvec, Cmp **ylmvec, Cmp source) const
Computes multipole moments.
void get_ylm(int nylm, Cmp **ylmvec) const
Constructs spherical harmonics.
Base_val base
Bases on which the spectral expansion is performed.
const double * get_beta() const
Returns the pointer on the array beta.
Map & mp
Mapping associated with the star.
int get_nzone() const
Returns the number of domains.
double integrale() const
Computes the integral over all space of *this .
Cmp pow(const Cmp &, int)
Power .
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Tbl & set(int l)
Read/write of the value in a given domain.
Coord y
y coordinate centered on the grid
Coord x
x coordinate centered on the grid
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Coord z
z coordinate centered on the grid
Valeur va
The numerical value of the Cmp.