72 cout <<
"Constructing ylm" << endl;
79 for (
int l=0 ; l<nz ; l++) {
85 for (
int k=0 ; k<np ; k++) {
86 for (
int j=0 ; j<nt ; j++) {
87 for (
int i=0 ; i<nr ; i++) {
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);
97 ylmvec[0]->
set(l,k,j,i)=1.0*
sqrt(1.0/4.0/M_PI);
101 if (nylm <4) {abort();}
else {
103 ylmvec[1]->
set(l,k,j,i)=zval*
sqrt(3.0/4.0/M_PI);
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);
111 if (nylm <9) {abort();}
else {
113 ylmvec[4]->
set(l,k,j,i)=(3.0*zval*zval-rval*rval)*
sqrt(5.0/16.0/M_PI);
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);
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);
124 if (nylm <16) {abort();}
else {
126 ylmvec[9]->
set(l,k,j,i)=(5.0*
pow(zval,3)-3.0*zval*rval*rval)*
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);
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);
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);
147 if (nylm <25) {abort();}
else {
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);
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);
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);
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);
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);
175 if (nylm <36) {abort();}
else {
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);
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);
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);
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);
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);
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);
208 if (nylm <49) {abort();}
else {
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);
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);
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);
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);
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);
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);
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);
245 cout <<
"l>6 not implemented!!!!!!!"<< endl;