9 #ifdef TIDAL_TENSOR_3RD
17 T2{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
18 #ifdef TIDAL_TENSOR_3RD
19 T3{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
30 std::cerr<<
"Error: Data reading fails! requiring data number is 1, only obtain "<<rcount<<
".\n";
36 T1[0] = T1[1] = T1[2] = 0;
37 for(
PS::S32 i=0; i<9; i++) T2[i] = 0;
38 #ifdef TIDAL_TENSOR_3RD
39 for(
PS::S32 i=0; i<10; i++) T3[i] = 0;
52 #ifdef TIDAL_TENSOR_3RD
60 _ptcl_tt[0].pos =
PS::F64vec(lscale, 0, -lscale) + _ptcl_cm.pos;
61 _ptcl_tt[1].pos =
PS::F64vec(0, lscale, -lscale) + _ptcl_cm.pos;
62 _ptcl_tt[2].pos =
PS::F64vec(-lscale, 0, -lscale) + _ptcl_cm.pos;
63 _ptcl_tt[3].pos =
PS::F64vec(0, -lscale, -lscale) + _ptcl_cm.pos;
64 _ptcl_tt[4].pos =
PS::F64vec(lscale, 0, lscale) + _ptcl_cm.pos;
65 _ptcl_tt[5].pos =
PS::F64vec(0, lscale, lscale) + _ptcl_cm.pos;
66 _ptcl_tt[6].pos =
PS::F64vec(-lscale, 0, lscale) + _ptcl_cm.pos;
67 _ptcl_tt[7].pos =
PS::F64vec(0, -lscale, lscale) + _ptcl_cm.pos;
77 _ptcl_tt[0].pos =
PS::F64vec( lscale, 0, -lscale*0.707106781186548) + _ptcl_cm.pos;
78 _ptcl_tt[1].pos =
PS::F64vec(-lscale, 0, -lscale*0.707106781186548) + _ptcl_cm.pos;
79 _ptcl_tt[2].pos =
PS::F64vec(0, lscale, lscale*0.707106781186548) + _ptcl_cm.pos;
80 _ptcl_tt[3].pos =
PS::F64vec(0, -lscale, lscale*0.707106781186548) + _ptcl_cm.pos;
86 _ptcl_tt[i].vel = _ptcl_cm.vel;
88 _ptcl_tt[i].mass = 0.0;
95 for (
int k=0; k<
getParticleN(); k++) _ptcl_tt[k].acc -= _ptcl_cm.acc;
127 template<
class Tptcl>
128 void fit(Tptcl* _ptcl_tt, Tptcl& _ptcl_cm,
const PS::F64 _size) {
137 for (
PS::S32 i=0; i<n_point; i++) fi[i] = _ptcl_tt[i].acc;
141 T1[0] = T1[1] = T1[2] = 0.0;
143 #ifdef TIDAL_TENSOR_3RD
146 T2[0] = 0.250000000000000*fi[0][0] + -0.250000000000000*fi[2][0] + 0.250000000000000*fi[4][0] + -0.250000000000000*fi[6][0];
147 T2[1] = 0.125000000000000*fi[0][1] + 0.125000000000000*fi[1][0] + -0.125000000000000*fi[2][1] + -0.125000000000000*fi[3][0]
148 + 0.125000000000000*fi[4][1] + 0.125000000000000*fi[5][0] + -0.125000000000000*fi[6][1] + -0.125000000000000*fi[7][0];
149 T2[2] = -0.083333333333333*fi[0][0] + 0.083333333333333*fi[0][2] + -0.083333333333333*fi[1][0] + -0.083333333333333*fi[2][0]
150 + -0.083333333333333*fi[2][2] + -0.083333333333333*fi[3][0] + 0.083333333333333*fi[4][0] + 0.083333333333333*fi[4][2]
151 + 0.083333333333333*fi[5][0] + 0.083333333333333*fi[6][0] + -0.083333333333333*fi[6][2] + 0.083333333333333*fi[7][0];
155 T2[4] = 0.250000000000000*fi[1][1] + -0.250000000000000*fi[3][1] + 0.250000000000000*fi[5][1] + -0.250000000000000*fi[7][1];
156 T2[5] = -0.083333333333333*fi[0][1] + -0.083333333333333*fi[1][1] + 0.083333333333333*fi[1][2] + -0.083333333333333*fi[2][1]
157 + -0.083333333333333*fi[3][1] + -0.083333333333333*fi[3][2] + 0.083333333333333*fi[4][1] + 0.083333333333333*fi[5][1]
158 + 0.083333333333333*fi[5][2] + 0.083333333333333*fi[6][1] + 0.083333333333333*fi[7][1] + -0.083333333333333*fi[7][2];
163 T2[8] = -0.125000000000000*fi[0][2] + -0.125000000000000*fi[1][2] + -0.125000000000000*fi[2][2] + -0.125000000000000*fi[3][2]
164 + 0.125000000000000*fi[4][2] + 0.125000000000000*fi[5][2] + 0.125000000000000*fi[6][2] + 0.125000000000000*fi[7][2];
167 T3[0] = 0.250000000000000*fi[0][0] + 0.125000000000000*fi[0][2] + 0.250000000000000*fi[2][0] + -0.125000000000000*fi[2][2]
168 + 0.250000000000000*fi[4][0] + -0.125000000000000*fi[4][2] + 0.250000000000000*fi[6][0] + 0.125000000000000*fi[6][2];
169 T3[1] = 0.250000000000000*fi[0][1] + 0.125000000000000*fi[1][2] + 0.250000000000000*fi[2][1] + -0.125000000000000*fi[3][2]
170 + 0.250000000000000*fi[4][1] + -0.125000000000000*fi[5][2] + 0.250000000000000*fi[6][1] + 0.125000000000000*fi[7][2];
171 T3[2] = -0.112500000000000*fi[0][0] + 0.025000000000000*fi[0][2] + -0.012500000000000*fi[1][1] + -0.025000000000000*fi[1][2]
172 + 0.112500000000000*fi[2][0] + 0.025000000000000*fi[2][2] + 0.012500000000000*fi[3][1] + -0.025000000000000*fi[3][2]
173 + 0.112500000000000*fi[4][0] + 0.025000000000000*fi[4][2] + 0.012500000000000*fi[5][1] + -0.025000000000000*fi[5][2]
174 + -0.112500000000000*fi[6][0] + 0.025000000000000*fi[6][2] + -0.012500000000000*fi[7][1] + -0.025000000000000*fi[7][2];
175 T3[3] = 0.125000000000000*fi[0][2] + 0.250000000000000*fi[1][0] + -0.125000000000000*fi[2][2] + 0.250000000000000*fi[3][0]
176 + -0.125000000000000*fi[4][2] + 0.250000000000000*fi[5][0] + 0.125000000000000*fi[6][2] + 0.250000000000000*fi[7][0];
177 T3[4] = -0.062500000000000*fi[0][1] + -0.062500000000000*fi[1][0] + 0.062500000000000*fi[2][1] + 0.062500000000000*fi[3][0]
178 + 0.062500000000000*fi[4][1] + 0.062500000000000*fi[5][0] + -0.062500000000000*fi[6][1] + -0.062500000000000*fi[7][0];
179 T3[5] = -0.125000000000000*fi[0][2] + 0.125000000000000*fi[2][2] + 0.125000000000000*fi[4][2] + -0.125000000000000*fi[6][2];
180 T3[6] = 0.250000000000000*fi[1][1] + 0.125000000000000*fi[1][2] + 0.250000000000000*fi[3][1] + -0.125000000000000*fi[3][2]
181 + 0.250000000000000*fi[5][1] + -0.125000000000000*fi[5][2] + 0.250000000000000*fi[7][1] + 0.125000000000000*fi[7][2];
182 T3[7] = -0.012500000000000*fi[0][0] + -0.025000000000000*fi[0][2] + -0.112500000000000*fi[1][1] + 0.025000000000000*fi[1][2]
183 + 0.012500000000000*fi[2][0] + -0.025000000000000*fi[2][2] + 0.112500000000000*fi[3][1] + 0.025000000000000*fi[3][2]
184 + 0.012500000000000*fi[4][0] + -0.025000000000000*fi[4][2] + 0.112500000000000*fi[5][1] + 0.025000000000000*fi[5][2]
185 + -0.012500000000000*fi[6][0] + -0.025000000000000*fi[6][2] + -0.112500000000000*fi[7][1] + 0.025000000000000*fi[7][2];
186 T3[8] = -0.125000000000000*fi[1][2] + 0.125000000000000*fi[3][2] + 0.125000000000000*fi[5][2] + -0.125000000000000*fi[7][2];
187 T3[9] = 0.062500000000000*fi[0][0] + 0.125000000000000*fi[0][2] + 0.062500000000000*fi[1][1] + 0.125000000000000*fi[1][2]
188 + -0.062500000000000*fi[2][0] + 0.125000000000000*fi[2][2] + -0.062500000000000*fi[3][1] + 0.125000000000000*fi[3][2]
189 + -0.062500000000000*fi[4][0] + 0.125000000000000*fi[4][2] + -0.062500000000000*fi[5][1] + 0.125000000000000*fi[5][2]
190 + 0.062500000000000*fi[6][0] + 0.125000000000000*fi[6][2] + 0.062500000000000*fi[7][1] + 0.125000000000000*fi[7][2];
194 PS::F64 T2S = 1.0/(_size*0.16);
196 for (
PS::S32 i=0; i<9; i++) T2[i] *= T2S;
197 for (
PS::S32 i=0; i<10; i++) T3[i] *= T3S;
203 T2[0] = 0.500000000000000*fi[0][0]+ -0.500000000000000*fi[1][0];
204 T2[1] = 0.250000000000000*fi[0][1]+ -0.250000000000000*fi[1][1]+ 0.250000000000000*fi[2][0]+ -0.250000000000000*fi[3][0];
205 T2[2] = -0.176776695296637*fi[0][0]+ 0.250000000000000*fi[0][2]+ -0.176776695296637*fi[1][0]+ -0.250000000000000*fi[1][2]+ 0.176776695296637*fi[2][0]+ 0.176776695296637*fi[3][0];
209 T2[4] = 0.500000000000000*fi[2][1]+ -0.500000000000000*fi[3][1];
210 T2[5] = -0.176776695296637*fi[0][1]+ -0.176776695296637*fi[1][1]+ 0.176776695296637*fi[2][1]+ 0.250000000000000*fi[2][2]+ 0.176776695296637*fi[3][1]+ -0.250000000000000*fi[3][2];
215 T2[8] = -0.353553390593274*fi[0][2]+ -0.353553390593274*fi[1][2]+ 0.353553390593274*fi[2][2]+ 0.353553390593274*fi[3][2];
219 for (
PS::S32 i=0; i<9; i++) T2[i] *= T2S;
248 #ifdef TIDAL_TENSOR_3RD
272 T2[0] += 2.0*(T3[0]*x + T3[1]*y + T3[2]*z);
273 T2[1] += 2.0*(T3[1]*x + T3[3]*y + T3[4]*z);
274 T2[2] += 2.0*(T3[2]*x + T3[4]*y + T3[5]*z);
276 T2[3] += 2.0*(T3[1]*x + T3[3]*y + T3[4]*z);
277 T2[4] += 2.0*(T3[3]*x + T3[6]*y + T3[7]*z);
278 T2[5] += 2.0*(T3[4]*x + T3[7]*y + T3[8]*z);
280 T2[6] += 2.0*(T3[2]*x + T3[4]*y + T3[5]*z);
281 T2[7] += 2.0*(T3[4]*x + T3[7]*y + T3[8]*z);
282 T2[8] += 2.0*(T3[5]*x + T3[8]*y + T3[9]*z);
298 #ifdef TIDAL_TENSOR_3RD
326 acc0 += T1[0] + T2[0]*x + T2[1]*y + T2[2]*z
327 + T3[0]*x2 + 2*T3[1]*xy + 2*T3[2]*xz + T3[3]*y2 + 2*T3[4]*yz + T3[5]*z2;
328 acc1 += T1[1] + T2[3]*x + T2[4]*y + T2[5]*z
329 + T3[1]*x2 + 2*T3[3]*xy + 2*T3[4]*xz + T3[6]*y2 + 2*T3[7]*yz + T3[8]*z2;
330 acc2 += T1[2] + T2[6]*x + T2[7]*y + T2[8]*z
331 + T3[2]*x2 + 2*T3[4]*xy + 2*T3[5]*xz + T3[7]*y2 + 2*T3[8]*yz + T3[9]*z2;
334 acc0 += T1[0] + T2[0]*x + T2[1]*y + T2[2]*z;
335 acc1 += T1[1] + T2[3]*x + T2[4]*y + T2[5]*z;
336 acc2 += T1[2] + T2[6]*x + T2[7]*y + T2[8]*z;
348 #ifdef TIDAL_TENSOR_3RD
356 PS::F64 acc0 = T1[0] + 0.5*(T2[0]*x + T2[1]*y + T2[2]*z)
357 + (T3[0]*x2 + 2*T3[1]*xy + 2*T3[2]*xz + T3[3]*y2 + 2*T3[4]*yz + T3[5]*z2)/3.0;
358 PS::F64 acc1 = T1[1] + 0.5*(T2[3]*x + T2[4]*y + T2[5]*z)
359 + (T3[1]*x2 + 2*T3[3]*xy + 2*T3[4]*xz + T3[6]*y2 + 2*T3[7]*yz + T3[8]*z2)/3.0;
360 PS::F64 acc2 = T1[2] + 0.5*(T2[6]*x + T2[7]*y + T2[8]*z)
361 + (T3[2]*x2 + 2*T3[4]*xy + 2*T3[5]*xz + T3[7]*y2 + 2*T3[8]*yz + T3[9]*z2)/3.0;
363 PS::F64 acc0 = T1[0] + 0.5*(T2[0]*x + T2[1]*y + T2[2]*z);
364 PS::F64 acc1 = T1[1] + 0.5*(T2[3]*x + T2[4]*y + T2[5]*z);
365 PS::F64 acc2 = T1[2] + 0.5*(T2[6]*x + T2[7]*y + T2[8]*z);
368 return - x*acc0 - y*acc1 - z*acc2;
371 void print(std::ostream & _fout,
const int _width)
const{
373 <<std::setw(_width)<<T1[0]
374 <<std::setw(_width)<<T1[1]
375 <<std::setw(_width)<<T1[2]
378 <<std::setw(_width)<<T2[0]
379 <<std::setw(_width)<<T2[1]
380 <<std::setw(_width)<<T2[2]
382 <<std::setw(_width)<<T2[3]
383 <<std::setw(_width)<<T2[4]
384 <<std::setw(_width)<<T2[5]
386 <<std::setw(_width)<<T2[6]
387 <<std::setw(_width)<<T2[7]
388 <<std::setw(_width)<<T2[8]
389 #ifdef TIDAL_TENSOR_3RD
393 <<std::setw(_width)<<T3[0]
394 <<std::setw(_width)<<T3[1]
395 <<std::setw(_width)<<T3[2]
397 <<std::setw(_width)<<T3[1]
398 <<std::setw(_width)<<T3[3]
399 <<std::setw(_width)<<T3[4]
401 <<std::setw(_width)<<T3[2]
402 <<std::setw(_width)<<T3[4]
403 <<std::setw(_width)<<T3[5]
406 <<std::setw(_width)<<T3[1]
407 <<std::setw(_width)<<T3[3]
408 <<std::setw(_width)<<T3[4]
410 <<std::setw(_width)<<T3[3]
411 <<std::setw(_width)<<T3[6]
412 <<std::setw(_width)<<T3[7]
414 <<std::setw(_width)<<T3[4]
415 <<std::setw(_width)<<T3[7]
416 <<std::setw(_width)<<T3[8]
419 <<std::setw(_width)<<T3[2]
420 <<std::setw(_width)<<T3[4]
421 <<std::setw(_width)<<T3[5]
423 <<std::setw(_width)<<T3[4]
424 <<std::setw(_width)<<T3[7]
425 <<std::setw(_width)<<T3[8]
427 <<std::setw(_width)<<T3[5]
428 <<std::setw(_width)<<T3[8]
429 <<std::setw(_width)<<T3[9]
436 #ifdef TIDAL_TENSOR_3RD