multipole.h 93.8 KB
Newer Older
Pedro Gonnet's avatar
Pedro Gonnet committed
1
2
/*******************************************************************************
 * This file is part of SWIFT.
3
 * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
4
 *               2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
5
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
6
7
8
9
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
10
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
11
12
13
14
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
15
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
16
17
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
18
 *
Pedro Gonnet's avatar
Pedro Gonnet committed
19
 ******************************************************************************/
20
21
#ifndef SWIFT_MULTIPOLE_H
#define SWIFT_MULTIPOLE_H
Pedro Gonnet's avatar
Pedro Gonnet committed
22

23
24
25
/* Config parameters. */
#include "../config.h"

Pedro Gonnet's avatar
Pedro Gonnet committed
26
27
/* Some standard headers. */
#include <math.h>
28
#include <string.h>
Pedro Gonnet's avatar
Pedro Gonnet committed
29

30
/* Includes. */
31
#include "align.h"
32
#include "const.h"
33
34
#include "error.h"
#include "gravity_derivatives.h"
35
#include "gravity_properties.h"
36
#include "gravity_softened_derivatives.h"
Pedro Gonnet's avatar
Pedro Gonnet committed
37
#include "inline.h"
38
#include "kernel_gravity.h"
39
#include "part.h"
40
#include "periodic.h"
41
#include "vector_power.h"
42

43
44
#define multipole_align 128

45
struct grav_tensor {
46

47
  /* 0th order terms */
48
  float F_000;
49

50
51
52
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0

  /* 1st order terms */
53
  float F_100, F_010, F_001;
54
55
56
57
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1

  /* 2nd order terms */
58
59
  float F_200, F_020, F_002;
  float F_110, F_101, F_011;
60
61
62
63
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2

  /* 3rd order terms */
64
65
66
67
68
  float F_300, F_030, F_003;
  float F_210, F_201;
  float F_120, F_021;
  float F_102, F_012;
  float F_111;
69
70
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
71
72

  /* 4th order terms */
73
74
75
76
77
78
  float F_400, F_040, F_004;
  float F_310, F_301;
  float F_130, F_031;
  float F_103, F_013;
  float F_220, F_202, F_022;
  float F_211, F_121, F_112;
79
80
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
81
82

  /* 5th order terms */
83
84
85
86
87
88
89
  float F_005, F_014, F_023;
  float F_032, F_041, F_050;
  float F_104, F_113, F_122;
  float F_131, F_140, F_203;
  float F_212, F_221, F_230;
  float F_302, F_311, F_320;
  float F_401, F_410, F_500;
90
91
92
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
93
94
95
#endif
#ifdef SWIFT_DEBUG_CHECKS

96
97
  /* Total number of gpart this field tensor interacted with */
  long long num_interacted;
98

99
100
101
  /* Last time this tensor was zeroed */
  integertime_t ti_init;

102
#endif
103
104
105

  /* Has this tensor received any contribution? */
  char interacted;
106
107
};

108
struct multipole {
109

110
  /* Bulk velocity */
111
  float vel[3];
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136

  /* 0th order terms */
  float M_000;

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0

  /* 1st order terms */
  float M_100, M_010, M_001;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1

  /* 2nd order terms */
  float M_200, M_020, M_002;
  float M_110, M_101, M_011;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2

  /* 3rd order terms */
  float M_300, M_030, M_003;
  float M_210, M_201;
  float M_120, M_021;
  float M_102, M_012;
  float M_111;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
137
138
139
140
141
142
143
144
145
146

  /* 4th order terms */
  float M_400, M_040, M_004;
  float M_310, M_301;
  float M_130, M_031;
  float M_103, M_013;
  float M_220, M_202, M_022;
  float M_211, M_121, M_112;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
147
148

  /* 5th order terms */
149
150
151
152
153
154
155
  float M_005, M_014, M_023;
  float M_032, M_041, M_050;
  float M_104, M_113, M_122;
  float M_131, M_140, M_203;
  float M_212, M_221, M_230;
  float M_302, M_311, M_320;
  float M_401, M_410, M_500;
156
157
158
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
159
#endif
160
161
162
163
164
165
166

#ifdef SWIFT_DEBUG_CHECKS

  /* Total number of gpart in this multipole */
  long long num_gpart;

#endif
167
};
168

169
170
171
172
/**
 * @brief The multipole expansion of a mass distribution.
 */
struct gravity_tensors {
173

174
  union {
175

176
177
    /*! Linking pointer for "memory management". */
    struct gravity_tensors *next;
178

179
180
    /*! The actual content */
    struct {
181

182
183
184
185
186
187
      /*! Multipole mass */
      struct multipole m_pole;

      /*! Field tensor for the potential */
      struct grav_tensor pot;

188
189
      /*! Centre of mass of the matter dsitribution */
      double CoM[3];
190

191
192
193
      /*! Centre of mass of the matter dsitribution at the last time-step */
      double CoM_old[3];

194
195
196
      /*! Centre of mass of the matter dsitribution at the last rebuild */
      double CoM_rebuild[3];

197
198
199
      /*! Upper limit of the CoM<->gpart distance */
      double r_max;

200
201
202
      /*! Upper limit of the CoM<->gpart distance at the last time-step */
      double r_max_old;

203
204
      /*! Upper limit of the CoM<->gpart distance at the last rebuild */
      double r_max_rebuild;
205
    };
206
  };
207
208
209
210
211
212
213
} SWIFT_STRUCT_ALIGN;

/**
 * @brief Reset the data of a #multipole.
 *
 * @param m The #multipole.
 */
214
INLINE static void gravity_reset(struct gravity_tensors *m) {
215
216

  /* Just bzero the struct. */
217
218
219
  bzero(m, sizeof(struct gravity_tensors));
}

220
221
222
223
224
/**
 * @brief Drifts a #multipole forward in time.
 *
 * @param m The #multipole.
 * @param dt The drift time-step.
225
226
 * @param x_diff The maximal distance moved by any particle since the last
 * rebuild.
227
 */
228
229
INLINE static void gravity_drift(struct gravity_tensors *m, double dt,
                                 float x_diff) {
230

231
232
233
  const double dx = m->m_pole.vel[0] * dt;
  const double dy = m->m_pole.vel[1] * dt;
  const double dz = m->m_pole.vel[2] * dt;
234

235
  /* Move the whole thing according to bulk motion */
236
237
238
  m->CoM[0] += dx;
  m->CoM[1] += dy;
  m->CoM[2] += dz;
239
240

  /* Conservative change in maximal radius containing all gpart */
241
  m->r_max = m->r_max_rebuild + 0. * x_diff;
242
243
}

244
245
246
247
/**
 * @brief Zeroes all the fields of a field tensor
 *
 * @param l The field tensor.
248
 * @param ti_current The current (integer) time (for debugging only).
249
 */
250
251
INLINE static void gravity_field_tensors_init(struct grav_tensor *l,
                                              integertime_t ti_current) {
252

253
  bzero(l, sizeof(struct grav_tensor));
254
255
256
257

#ifdef SWIFT_DEBUG_CHECKS
  l->ti_init = ti_current;
#endif
258
259
}

260
/**
261
 * @brief Adds a field tensor to another one (i.e. does la += lb).
262
263
264
265
 *
 * @param la The gravity tensors to add to.
 * @param lb The gravity tensors to add.
 */
266
267
INLINE static void gravity_field_tensors_add(struct grav_tensor *la,
                                             const struct grav_tensor *lb) {
268
#ifdef SWIFT_DEBUG_CHECKS
269
270
  if (lb->num_interacted == 0) error("Adding tensors that did not interact");
  la->num_interacted += lb->num_interacted;
271
#endif
272

273
274
  la->interacted = 1;

275
  /* Add 0th order terms */
276
  la->F_000 += lb->F_000;
277
278
279

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  /* Add 1st order terms */
280
281
282
  la->F_100 += lb->F_100;
  la->F_010 += lb->F_010;
  la->F_001 += lb->F_001;
283
284
285
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
  /* Add 2nd order terms */
286
287
288
289
290
291
  la->F_200 += lb->F_200;
  la->F_020 += lb->F_020;
  la->F_002 += lb->F_002;
  la->F_110 += lb->F_110;
  la->F_101 += lb->F_101;
  la->F_011 += lb->F_011;
292
293
294
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
  /* Add 3rd order terms */
295
296
297
298
299
300
301
302
303
304
  la->F_300 += lb->F_300;
  la->F_030 += lb->F_030;
  la->F_003 += lb->F_003;
  la->F_210 += lb->F_210;
  la->F_201 += lb->F_201;
  la->F_120 += lb->F_120;
  la->F_021 += lb->F_021;
  la->F_102 += lb->F_102;
  la->F_012 += lb->F_012;
  la->F_111 += lb->F_111;
305
#endif
306
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
  /* Add 4th order terms */
  la->F_400 += lb->F_400;
  la->F_040 += lb->F_040;
  la->F_004 += lb->F_004;
  la->F_310 += lb->F_310;
  la->F_301 += lb->F_301;
  la->F_130 += lb->F_130;
  la->F_031 += lb->F_031;
  la->F_103 += lb->F_103;
  la->F_013 += lb->F_013;
  la->F_220 += lb->F_220;
  la->F_202 += lb->F_202;
  la->F_022 += lb->F_022;
  la->F_211 += lb->F_211;
  la->F_121 += lb->F_121;
  la->F_112 += lb->F_112;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
  /* 5th order terms */
  la->F_005 += lb->F_005;
  la->F_014 += lb->F_014;
  la->F_023 += lb->F_023;
  la->F_032 += lb->F_032;
  la->F_041 += lb->F_041;
  la->F_050 += lb->F_050;
  la->F_104 += lb->F_104;
  la->F_113 += lb->F_113;
  la->F_122 += lb->F_122;
  la->F_131 += lb->F_131;
  la->F_140 += lb->F_140;
  la->F_203 += lb->F_203;
  la->F_212 += lb->F_212;
  la->F_221 += lb->F_221;
  la->F_230 += lb->F_230;
  la->F_302 += lb->F_302;
  la->F_311 += lb->F_311;
  la->F_320 += lb->F_320;
  la->F_401 += lb->F_401;
  la->F_410 += lb->F_410;
  la->F_500 += lb->F_500;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
350
#endif
351
352
353
354
355
356
357
358
359
360
361
362
}

/**
 * @brief Prints the content of a #grav_tensor to stdout.
 *
 * Note: Uses directly printf(), not a call to message().
 *
 * @param l The #grav_tensor to print.
 */
INLINE static void gravity_field_tensors_print(const struct grav_tensor *l) {

  printf("-------------------------\n");
363
  printf("Interacted: %d\n", l->interacted);
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
  printf("F_000= %12.5e\n", l->F_000);
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  printf("-------------------------\n");
  printf("F_100= %12.5e F_010= %12.5e F_001= %12.5e\n", l->F_100, l->F_010,
         l->F_001);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
  printf("-------------------------\n");
  printf("F_200= %12.5e F_020= %12.5e F_002= %12.5e\n", l->F_200, l->F_020,
         l->F_002);
  printf("F_110= %12.5e F_101= %12.5e F_011= %12.5e\n", l->F_110, l->F_101,
         l->F_011);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
  printf("-------------------------\n");
  printf("F_300= %12.5e F_030= %12.5e F_003= %12.5e\n", l->F_300, l->F_030,
         l->F_003);
  printf("F_210= %12.5e F_201= %12.5e F_120= %12.5e\n", l->F_210, l->F_201,
         l->F_120);
  printf("F_021= %12.5e F_102= %12.5e F_012= %12.5e\n", l->F_021, l->F_102,
         l->F_012);
  printf("F_111= %12.5e\n", l->F_111);
#endif
387
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
388
389
390
391
392
393
394
395
396
397
398
399
400
  printf("-------------------------\n");
  printf("F_400= %12.5e F_040= %12.5e F_004= %12.5e\n", l->F_400, l->F_040,
         l->F_004);
  printf("F_310= %12.5e F_301= %12.5e F_130= %12.5e\n", l->F_310, l->F_301,
         l->F_130);
  printf("F_031= %12.5e F_103= %12.5e F_013= %12.5e\n", l->F_031, l->F_103,
         l->F_013);
  printf("F_220= %12.5e F_202= %12.5e F_022= %12.5e\n", l->F_220, l->F_202,
         l->F_022);
  printf("F_211= %12.5e F_121= %12.5e F_112= %12.5e\n", l->F_211, l->F_121,
         l->F_112);
#endif
  printf("-------------------------\n");
401
402
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
403
#endif
404
405
}

406
407
408
409
410
/**
 * @brief Zeroes all the fields of a multipole.
 *
 * @param m The multipole
 */
411
412
413
414
415
INLINE static void gravity_multipole_init(struct multipole *m) {

  bzero(m, sizeof(struct multipole));
}

416
417
418
419
420
421
422
/**
 * @brief Prints the content of a #multipole to stdout.
 *
 * Note: Uses directly printf(), not a call to message().
 *
 * @param m The #multipole to print.
 */
423
INLINE static void gravity_multipole_print(const struct multipole *m) {
424

425
426
  printf("Vel= [%12.5e %12.5e %12.5e]\n", m->vel[0], m->vel[1], m->vel[2]);
  printf("-------------------------\n");
427
428
  printf("M_000= %12.5e\n", m->M_000);
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
429
  printf("-------------------------\n");
430
  printf("M_100= %12.5e M_010= %12.5e M_001= %12.5e\n", m->M_100, m->M_010,
431
432
433
         m->M_001);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
434
  printf("-------------------------\n");
435
  printf("M_200= %12.5e M_020= %12.5e M_002= %12.5e\n", m->M_200, m->M_020,
436
         m->M_002);
437
  printf("M_110= %12.5e M_101= %12.5e M_011= %12.5e\n", m->M_110, m->M_101,
438
439
440
         m->M_011);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
441
  printf("-------------------------\n");
442
  printf("M_300= %12.5e M_030= %12.5e M_003= %12.5e\n", m->M_300, m->M_030,
443
         m->M_003);
444
  printf("M_210= %12.5e M_201= %12.5e M_120= %12.5e\n", m->M_210, m->M_201,
445
         m->M_120);
446
  printf("M_021= %12.5e M_102= %12.5e M_012= %12.5e\n", m->M_021, m->M_102,
447
         m->M_012);
448
  printf("M_111= %12.5e\n", m->M_111);
449
#endif
450
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
451
452
453
454
455
456
457
458
459
460
461
462
463
  printf("-------------------------\n");
  printf("M_400= %12.5e M_040= %12.5e M_004= %12.5e\n", m->M_400, m->M_040,
         m->M_004);
  printf("M_310= %12.5e M_301= %12.5e M_130= %12.5e\n", m->M_310, m->M_301,
         m->M_130);
  printf("M_031= %12.5e M_103= %12.5e M_013= %12.5e\n", m->M_031, m->M_103,
         m->M_013);
  printf("M_220= %12.5e M_202= %12.5e M_022= %12.5e\n", m->M_220, m->M_202,
         m->M_022);
  printf("M_211= %12.5e M_121= %12.5e M_112= %12.5e\n", m->M_211, m->M_121,
         m->M_112);
#endif
  printf("-------------------------\n");
464
465
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
466
#endif
467
468
469
470
471
472
473
474
}

/**
 * @brief Adds a #multipole to another one (i.e. does ma += mb).
 *
 * @param ma The multipole to add to.
 * @param mb The multipole to add.
 */
475
476
INLINE static void gravity_multipole_add(struct multipole *ma,
                                         const struct multipole *mb) {
477

478
479
  const float M_000 = ma->M_000 + mb->M_000;
  const float inv_M_000 = 1.f / M_000;
480
481

  /* Add the bulk velocities */
482
483
484
  ma->vel[0] = (ma->vel[0] * ma->M_000 + mb->vel[0] * mb->M_000) * inv_M_000;
  ma->vel[1] = (ma->vel[1] * ma->M_000 + mb->vel[1] * mb->M_000) * inv_M_000;
  ma->vel[2] = (ma->vel[2] * ma->M_000 + mb->vel[2] * mb->M_000) * inv_M_000;
485

486
487
  /* Add 0th order terms */
  ma->M_000 = M_000;
488
489
490
491
492
493
494

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  /* Add 1st order terms */
  ma->M_100 += mb->M_100;
  ma->M_010 += mb->M_010;
  ma->M_001 += mb->M_001;
#endif
495
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
496
  /* Add 2nd order terms */
497
498
499
500
501
502
503
504
  ma->M_200 += mb->M_200;
  ma->M_020 += mb->M_020;
  ma->M_002 += mb->M_002;
  ma->M_110 += mb->M_110;
  ma->M_101 += mb->M_101;
  ma->M_011 += mb->M_011;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
505
  /* Add 3rd order terms */
506
507
508
509
510
511
512
513
514
515
516
  ma->M_300 += mb->M_300;
  ma->M_030 += mb->M_030;
  ma->M_003 += mb->M_003;
  ma->M_210 += mb->M_210;
  ma->M_201 += mb->M_201;
  ma->M_120 += mb->M_120;
  ma->M_021 += mb->M_021;
  ma->M_102 += mb->M_102;
  ma->M_012 += mb->M_012;
  ma->M_111 += mb->M_111;
#endif
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
  /* Add 4th order terms */
  ma->M_400 += mb->M_400;
  ma->M_040 += mb->M_040;
  ma->M_004 += mb->M_004;
  ma->M_310 += mb->M_310;
  ma->M_301 += mb->M_301;
  ma->M_130 += mb->M_130;
  ma->M_031 += mb->M_031;
  ma->M_103 += mb->M_103;
  ma->M_013 += mb->M_013;
  ma->M_220 += mb->M_220;
  ma->M_202 += mb->M_202;
  ma->M_022 += mb->M_022;
  ma->M_211 += mb->M_211;
  ma->M_121 += mb->M_121;
  ma->M_112 += mb->M_112;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
  /* 5th order terms */
  ma->M_005 += mb->M_005;
  ma->M_014 += mb->M_014;
  ma->M_023 += mb->M_023;
  ma->M_032 += mb->M_032;
  ma->M_041 += mb->M_041;
  ma->M_050 += mb->M_050;
  ma->M_104 += mb->M_104;
  ma->M_113 += mb->M_113;
  ma->M_122 += mb->M_122;
  ma->M_131 += mb->M_131;
  ma->M_140 += mb->M_140;
  ma->M_203 += mb->M_203;
  ma->M_212 += mb->M_212;
  ma->M_221 += mb->M_221;
  ma->M_230 += mb->M_230;
  ma->M_302 += mb->M_302;
  ma->M_311 += mb->M_311;
  ma->M_320 += mb->M_320;
  ma->M_401 += mb->M_401;
  ma->M_410 += mb->M_410;
  ma->M_500 += mb->M_500;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
561
#endif
562

563
564
565
566
567
  // MATTHIEU
  ma->M_100 = 0.f;
  ma->M_010 = 0.f;
  ma->M_001 = 0.f;

568
569
570
#ifdef SWIFT_DEBUG_CHECKS
  ma->num_gpart += mb->num_gpart;
#endif
571
572
573
574
575
}

/**
 * @brief Verifies whether two #multipole's are equal or not.
 *
576
577
 * @param ga The first #multipole.
 * @param gb The second #multipole.
578
 * @param tolerance The maximal allowed relative difference for the fields.
579
 * @return 1 if the multipoles are equal, 0 otherwise
580
 */
581
582
INLINE static int gravity_multipole_equal(const struct gravity_tensors *ga,
                                          const struct gravity_tensors *gb,
583
                                          double tolerance) {
584
585

  /* Check CoM */
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
  if (fabs(ga->CoM[0] - gb->CoM[0]) / fabs(ga->CoM[0] + gb->CoM[0]) >
      tolerance) {
    message("CoM[0] different");
    return 0;
  }
  if (fabs(ga->CoM[1] - gb->CoM[1]) / fabs(ga->CoM[1] + gb->CoM[1]) >
      tolerance) {
    message("CoM[1] different");
    return 0;
  }
  if (fabs(ga->CoM[2] - gb->CoM[2]) / fabs(ga->CoM[2] + gb->CoM[2]) >
      tolerance) {
    message("CoM[2] different");
    return 0;
  }
601
602
603
604

  /* Helper pointers */
  const struct multipole *ma = &ga->m_pole;
  const struct multipole *mb = &gb->m_pole;
605

606
607
608
609
  const double v2 = ma->vel[0] * ma->vel[0] + ma->vel[1] * ma->vel[1] +
                    ma->vel[2] * ma->vel[2];

  /* Check bulk velocity (if non-zero and component > 1% of norm)*/
610
  if (fabsf(ma->vel[0] + mb->vel[0]) > 1e-10 &&
611
      (ma->vel[0] * ma->vel[0]) > 0.0001 * v2 &&
612
      fabsf(ma->vel[0] - mb->vel[0]) / fabsf(ma->vel[0] + mb->vel[0]) >
613
614
615
616
          tolerance) {
    message("v[0] different");
    return 0;
  }
617
  if (fabsf(ma->vel[1] + mb->vel[1]) > 1e-10 &&
618
      (ma->vel[1] * ma->vel[1]) > 0.0001 * v2 &&
619
      fabsf(ma->vel[1] - mb->vel[1]) / fabsf(ma->vel[1] + mb->vel[1]) >
620
621
622
623
          tolerance) {
    message("v[1] different");
    return 0;
  }
624
  if (fabsf(ma->vel[2] + mb->vel[2]) > 1e-10 &&
625
      (ma->vel[2] * ma->vel[2]) > 0.0001 * v2 &&
626
      fabsf(ma->vel[2] - mb->vel[2]) / fabsf(ma->vel[2] + mb->vel[2]) >
627
628
629
630
          tolerance) {
    message("v[2] different");
    return 0;
  }
631

632
633
634
635
636
637
  /* Manhattan Norm of 0th order terms */
  const float order0_norm = fabsf(ma->M_000) + fabsf(mb->M_000);

  /* Compare 0th order terms above 1% of norm */
  if (fabsf(ma->M_000 + mb->M_000) > 0.01f * order0_norm &&
      fabsf(ma->M_000 - mb->M_000) / fabsf(ma->M_000 + mb->M_000) > tolerance) {
638
639
640
    message("M_000 term different");
    return 0;
  }
641
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
642
643
644
645
646
647
648
649
650
  /* Manhattan Norm of 1st order terms */
  const float order1_norm = fabsf(ma->M_001) + fabsf(mb->M_001) +
                            fabsf(ma->M_010) + fabsf(mb->M_010) +
                            fabsf(ma->M_100) + fabsf(mb->M_100);

  /* Compare 1st order terms above 1% of norm */
  if (fabsf(ma->M_001 + mb->M_001) > 0.01f * order1_norm &&
      fabsf(ma->M_001 - mb->M_001) / fabsf(ma->M_001 + mb->M_001) > tolerance) {
    message("M_001 term different");
651
652
    return 0;
  }
653
  if (fabsf(ma->M_010 + mb->M_010) > 0.01f * order1_norm &&
654
655
656
657
      fabsf(ma->M_010 - mb->M_010) / fabsf(ma->M_010 + mb->M_010) > tolerance) {
    message("M_010 term different");
    return 0;
  }
658
659
660
  if (fabsf(ma->M_100 + mb->M_100) > 0.01f * order1_norm &&
      fabsf(ma->M_100 - mb->M_100) / fabsf(ma->M_100 + mb->M_100) > tolerance) {
    message("M_100 term different");
661
662
    return 0;
  }
663
#endif
664
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
  /* Manhattan Norm of 2nd order terms */
  const float order2_norm =
      fabsf(ma->M_002) + fabsf(mb->M_002) + fabsf(ma->M_011) +
      fabsf(mb->M_011) + fabsf(ma->M_020) + fabsf(mb->M_020) +
      fabsf(ma->M_101) + fabsf(mb->M_101) + fabsf(ma->M_110) +
      fabsf(mb->M_110) + fabsf(ma->M_200) + fabsf(mb->M_200);

  /* Compare 2nd order terms above 1% of norm */
  if (fabsf(ma->M_002 + mb->M_002) > 0.01f * order2_norm &&
      fabsf(ma->M_002 - mb->M_002) / fabsf(ma->M_002 + mb->M_002) > tolerance) {
    message("M_002 term different");
    return 0;
  }
  if (fabsf(ma->M_011 + mb->M_011) > 0.01f * order2_norm &&
      fabsf(ma->M_011 - mb->M_011) / fabsf(ma->M_011 + mb->M_011) > tolerance) {
    message("M_011 term different");
681
682
    return 0;
  }
683
  if (fabsf(ma->M_020 + mb->M_020) > 0.01f * order2_norm &&
684
685
686
687
      fabsf(ma->M_020 - mb->M_020) / fabsf(ma->M_020 + mb->M_020) > tolerance) {
    message("M_020 term different");
    return 0;
  }
688
689
690
  if (fabsf(ma->M_101 + mb->M_101) > 0.01f * order2_norm &&
      fabsf(ma->M_101 - mb->M_101) / fabsf(ma->M_101 + mb->M_101) > tolerance) {
    message("M_101 term different");
691
692
    return 0;
  }
693
  if (fabsf(ma->M_110 + mb->M_110) > 0.01f * order2_norm &&
694
695
696
697
      fabsf(ma->M_110 - mb->M_110) / fabsf(ma->M_110 + mb->M_110) > tolerance) {
    message("M_110 term different");
    return 0;
  }
698
699
700
  if (fabsf(ma->M_200 + mb->M_200) > 0.01f * order2_norm &&
      fabsf(ma->M_200 - mb->M_200) / fabsf(ma->M_200 + mb->M_200) > tolerance) {
    message("M_200 term different");
701
702
    return 0;
  }
703
#endif
704
  tolerance *= 10.;
705
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
706
707
708
709
710
711
712
713
714
715
716
717
  /* Manhattan Norm of 3rd order terms */
  const float order3_norm =
      fabsf(ma->M_003) + fabsf(mb->M_003) + fabsf(ma->M_012) +
      fabsf(mb->M_012) + fabsf(ma->M_021) + fabsf(mb->M_021) +
      fabsf(ma->M_030) + fabsf(mb->M_030) + fabsf(ma->M_102) +
      fabsf(mb->M_102) + fabsf(ma->M_111) + fabsf(mb->M_111) +
      fabsf(ma->M_120) + fabsf(mb->M_120) + fabsf(ma->M_201) +
      fabsf(mb->M_201) + fabsf(ma->M_210) + fabsf(mb->M_210) +
      fabsf(ma->M_300) + fabsf(mb->M_300);

  /* Compare 3rd order terms above 1% of norm */
  if (fabsf(ma->M_003 + mb->M_003) > 0.01f * order3_norm &&
718
719
720
721
      fabsf(ma->M_003 - mb->M_003) / fabsf(ma->M_003 + mb->M_003) > tolerance) {
    message("M_003 term different");
    return 0;
  }
722
723
724
  if (fabsf(ma->M_012 + mb->M_012) > 0.01f * order3_norm &&
      fabsf(ma->M_012 - mb->M_012) / fabsf(ma->M_012 + mb->M_012) > tolerance) {
    message("M_012 term different");
725
726
    return 0;
  }
727
  if (fabsf(ma->M_021 + mb->M_021) > 0.01f * order3_norm &&
728
729
730
731
      fabsf(ma->M_021 - mb->M_021) / fabsf(ma->M_021 + mb->M_021) > tolerance) {
    message("M_021 term different");
    return 0;
  }
732
733
734
  if (fabsf(ma->M_030 + mb->M_030) > 0.01f * order3_norm &&
      fabsf(ma->M_030 - mb->M_030) / fabsf(ma->M_030 + mb->M_030) > tolerance) {
    message("M_030 term different");
735
736
    return 0;
  }
737
738
739
  if (fabsf(ma->M_102 + mb->M_102) > 0.01f * order3_norm &&
      fabsf(ma->M_102 - mb->M_102) / fabsf(ma->M_102 + mb->M_102) > tolerance) {
    message("M_102 term different");
740
741
    return 0;
  }
742
  if (fabsf(ma->M_111 + mb->M_111) > 0.01f * order3_norm &&
743
744
745
746
      fabsf(ma->M_111 - mb->M_111) / fabsf(ma->M_111 + mb->M_111) > tolerance) {
    message("M_111 term different");
    return 0;
  }
747
748
749
  if (fabsf(ma->M_120 + mb->M_120) > 0.01f * order3_norm &&
      fabsf(ma->M_120 - mb->M_120) / fabsf(ma->M_120 + mb->M_120) > tolerance) {
    message("M_120 term different");
750
751
    return 0;
  }
752
753
754
  if (fabsf(ma->M_201 + mb->M_201) > 0.01f * order3_norm &&
      fabsf(ma->M_201 - mb->M_201) / fabsf(ma->M_201 + mb->M_201) > tolerance) {
    message("M_201 term different");
755
756
    return 0;
  }
757
758
759
  if (fabsf(ma->M_210 + mb->M_210) > 0.01f * order3_norm &&
      fabsf(ma->M_210 - mb->M_210) / fabsf(ma->M_210 + mb->M_210) > tolerance) {
    message("M_210 term different");
760
761
    return 0;
  }
762
763
764
  if (fabsf(ma->M_300 + mb->M_300) > 0.01f * order3_norm &&
      fabsf(ma->M_300 - mb->M_300) / fabsf(ma->M_300 + mb->M_300) > tolerance) {
    message("M_300 term different");
765
766
767
    return 0;
  }
#endif
768
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
  /* Manhattan Norm of 4th order terms */
  const float order4_norm =
      fabsf(ma->M_004) + fabsf(mb->M_004) + fabsf(ma->M_013) +
      fabsf(mb->M_013) + fabsf(ma->M_022) + fabsf(mb->M_022) +
      fabsf(ma->M_031) + fabsf(mb->M_031) + fabsf(ma->M_040) +
      fabsf(mb->M_040) + fabsf(ma->M_103) + fabsf(mb->M_103) +
      fabsf(ma->M_112) + fabsf(mb->M_112) + fabsf(ma->M_121) +
      fabsf(mb->M_121) + fabsf(ma->M_130) + fabsf(mb->M_130) +
      fabsf(ma->M_202) + fabsf(mb->M_202) + fabsf(ma->M_211) +
      fabsf(mb->M_211) + fabsf(ma->M_220) + fabsf(mb->M_220) +
      fabsf(ma->M_301) + fabsf(mb->M_301) + fabsf(ma->M_310) +
      fabsf(mb->M_310) + fabsf(ma->M_400) + fabsf(mb->M_400);

  /* Compare 4th order terms above 1% of norm */
  if (fabsf(ma->M_004 + mb->M_004) > 0.01f * order4_norm &&
      fabsf(ma->M_004 - mb->M_004) / fabsf(ma->M_004 + mb->M_004) > tolerance) {
    message("M_004 term different");
    return 0;
  }
  if (fabsf(ma->M_013 + mb->M_013) > 0.01f * order4_norm &&
      fabsf(ma->M_013 - mb->M_013) / fabsf(ma->M_013 + mb->M_013) > tolerance) {
    message("M_013 term different");
    return 0;
  }
  if (fabsf(ma->M_022 + mb->M_022) > 0.01f * order4_norm &&
      fabsf(ma->M_022 - mb->M_022) / fabsf(ma->M_022 + mb->M_022) > tolerance) {
    message("M_022 term different");
    return 0;
  }
  if (fabsf(ma->M_031 + mb->M_031) > 0.01f * order4_norm &&
      fabsf(ma->M_031 - mb->M_031) / fabsf(ma->M_031 + mb->M_031) > tolerance) {
    message("M_031 term different");
    return 0;
  }
  if (fabsf(ma->M_040 + mb->M_040) > 0.01f * order4_norm &&
      fabsf(ma->M_040 - mb->M_040) / fabsf(ma->M_040 + mb->M_040) > tolerance) {
    message("M_040 term different");
    return 0;
  }
  if (fabsf(ma->M_103 + mb->M_103) > 0.01f * order4_norm &&
      fabsf(ma->M_103 - mb->M_103) / fabsf(ma->M_103 + mb->M_103) > tolerance) {
    message("M_103 term different");
    return 0;
  }
  if (fabsf(ma->M_112 + mb->M_112) > 0.01f * order4_norm &&
      fabsf(ma->M_112 - mb->M_112) / fabsf(ma->M_112 + mb->M_112) > tolerance) {
    message("M_112 term different");
    return 0;
  }
  if (fabsf(ma->M_121 + mb->M_121) > 0.01f * order4_norm &&
      fabsf(ma->M_121 - mb->M_121) / fabsf(ma->M_121 + mb->M_121) > tolerance) {
    message("M_121 term different");
    return 0;
  }
  if (fabsf(ma->M_130 + mb->M_130) > 0.01f * order4_norm &&
      fabsf(ma->M_130 - mb->M_130) / fabsf(ma->M_130 + mb->M_130) > tolerance) {
    message("M_130 term different");
    return 0;
  }
  if (fabsf(ma->M_202 + mb->M_202) > 0.01f * order4_norm &&
      fabsf(ma->M_202 - mb->M_202) / fabsf(ma->M_202 + mb->M_202) > tolerance) {
    message("M_202 term different");
    return 0;
  }
  if (fabsf(ma->M_211 + mb->M_211) > 0.01f * order4_norm &&
      fabsf(ma->M_211 - mb->M_211) / fabsf(ma->M_211 + mb->M_211) > tolerance) {
    message("M_211 term different");
    return 0;
  }
  if (fabsf(ma->M_220 + mb->M_220) > 0.01f * order4_norm &&
      fabsf(ma->M_220 - mb->M_220) / fabsf(ma->M_220 + mb->M_220) > tolerance) {
    message("M_220 term different");
    return 0;
  }
  if (fabsf(ma->M_301 + mb->M_301) > 0.01f * order4_norm &&
      fabsf(ma->M_301 - mb->M_301) / fabsf(ma->M_301 + mb->M_301) > tolerance) {
    message("M_301 term different");
    return 0;
  }
  if (fabsf(ma->M_310 + mb->M_310) > 0.01f * order4_norm &&
      fabsf(ma->M_310 - mb->M_310) / fabsf(ma->M_310 + mb->M_310) > tolerance) {
    message("M_310 term different");
    return 0;
  }
  if (fabsf(ma->M_400 + mb->M_400) > 0.01f * order4_norm &&
      fabsf(ma->M_400 - mb->M_400) / fabsf(ma->M_400 + mb->M_400) > tolerance) {
    message("M_400 term different");
    return 0;
  }
858
#endif
859
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
  /* Manhattan Norm of 5th order terms */
  const float order5_norm =
      fabsf(ma->M_005) + fabsf(mb->M_005) + fabsf(ma->M_014) +
      fabsf(mb->M_014) + fabsf(ma->M_023) + fabsf(mb->M_023) +
      fabsf(ma->M_032) + fabsf(mb->M_032) + fabsf(ma->M_041) +
      fabsf(mb->M_041) + fabsf(ma->M_050) + fabsf(mb->M_050) +
      fabsf(ma->M_104) + fabsf(mb->M_104) + fabsf(ma->M_113) +
      fabsf(mb->M_113) + fabsf(ma->M_122) + fabsf(mb->M_122) +
      fabsf(ma->M_131) + fabsf(mb->M_131) + fabsf(ma->M_140) +
      fabsf(mb->M_140) + fabsf(ma->M_203) + fabsf(mb->M_203) +
      fabsf(ma->M_212) + fabsf(mb->M_212) + fabsf(ma->M_221) +
      fabsf(mb->M_221) + fabsf(ma->M_230) + fabsf(mb->M_230) +
      fabsf(ma->M_302) + fabsf(mb->M_302) + fabsf(ma->M_311) +
      fabsf(mb->M_311) + fabsf(ma->M_320) + fabsf(mb->M_320) +
      fabsf(ma->M_401) + fabsf(mb->M_401) + fabsf(ma->M_410) +
      fabsf(mb->M_410) + fabsf(ma->M_500) + fabsf(mb->M_500);

  /* Compare 5th order terms above 1% of norm */
  if (fabsf(ma->M_005 + mb->M_005) > 0.01f * order5_norm &&
      fabsf(ma->M_005 - mb->M_005) / fabsf(ma->M_005 + mb->M_005) > tolerance) {
    message("M_005 term different");
    return 0;
  }
  if (fabsf(ma->M_014 + mb->M_014) > 0.01f * order5_norm &&
      fabsf(ma->M_014 - mb->M_014) / fabsf(ma->M_014 + mb->M_014) > tolerance) {
    message("M_014 term different");
    return 0;
  }
  if (fabsf(ma->M_023 + mb->M_023) > 0.01f * order5_norm &&
      fabsf(ma->M_023 - mb->M_023) / fabsf(ma->M_023 + mb->M_023) > tolerance) {
    message("M_023 term different");
    return 0;
  }
  if (fabsf(ma->M_032 + mb->M_032) > 0.01f * order5_norm &&
      fabsf(ma->M_032 - mb->M_032) / fabsf(ma->M_032 + mb->M_032) > tolerance) {
    message("M_032 term different");
    return 0;
  }
  if (fabsf(ma->M_041 + mb->M_041) > 0.01f * order5_norm &&
      fabsf(ma->M_041 - mb->M_041) / fabsf(ma->M_041 + mb->M_041) > tolerance) {
    message("M_041 term different");
    return 0;
  }
  if (fabsf(ma->M_050 + mb->M_050) > 0.01f * order5_norm &&
      fabsf(ma->M_050 - mb->M_050) / fabsf(ma->M_050 + mb->M_050) > tolerance) {
    message("M_050 term different");
    return 0;
  }
  if (fabsf(ma->M_104 + mb->M_104) > 0.01f * order5_norm &&
      fabsf(ma->M_104 - mb->M_104) / fabsf(ma->M_104 + mb->M_104) > tolerance) {
    message("M_104 term different");
    return 0;
  }
  if (fabsf(ma->M_113 + mb->M_113) > 0.01f * order5_norm &&
      fabsf(ma->M_113 - mb->M_113) / fabsf(ma->M_113 + mb->M_113) > tolerance) {
    message("M_113 term different");
    return 0;
  }
  if (fabsf(ma->M_122 + mb->M_122) > 0.01f * order5_norm &&
      fabsf(ma->M_122 - mb->M_122) / fabsf(ma->M_122 + mb->M_122) > tolerance) {
    message("M_122 term different");
    return 0;
  }
  if (fabsf(ma->M_131 + mb->M_131) > 0.01f * order5_norm &&
      fabsf(ma->M_131 - mb->M_131) / fabsf(ma->M_131 + mb->M_131) > tolerance) {
    message("M_131 term different");
    return 0;
  }
  if (fabsf(ma->M_140 + mb->M_140) > 0.01f * order5_norm &&
      fabsf(ma->M_140 - mb->M_140) / fabsf(ma->M_140 + mb->M_140) > tolerance) {
    message("M_140 term different");
    return 0;
  }
  if (fabsf(ma->M_203 + mb->M_203) > 0.01f * order5_norm &&
      fabsf(ma->M_203 - mb->M_203) / fabsf(ma->M_203 + mb->M_203) > tolerance) {
    message("M_203 term different");
    return 0;
  }
  if (fabsf(ma->M_212 + mb->M_212) > 0.01f * order5_norm &&
      fabsf(ma->M_212 - mb->M_212) / fabsf(ma->M_212 + mb->M_212) > tolerance) {
    message("M_212 term different");
    return 0;
  }
  if (fabsf(ma->M_221 + mb->M_221) > 0.01f * order5_norm &&
      fabsf(ma->M_221 - mb->M_221) / fabsf(ma->M_221 + mb->M_221) > tolerance) {
    message("M_221 term different");
    return 0;
  }
  if (fabsf(ma->M_230 + mb->M_230) > 0.01f * order5_norm &&
      fabsf(ma->M_230 - mb->M_230) / fabsf(ma->M_230 + mb->M_230) > tolerance) {
    message("M_230 term different");
    return 0;
  }
  if (fabsf(ma->M_302 + mb->M_302) > 0.01f * order5_norm &&
      fabsf(ma->M_302 - mb->M_302) / fabsf(ma->M_302 + mb->M_302) > tolerance) {
    message("M_302 term different");
    return 0;
  }
  if (fabsf(ma->M_311 + mb->M_311) > 0.01f * order5_norm &&
      fabsf(ma->M_311 - mb->M_311) / fabsf(ma->M_311 + mb->M_311) > tolerance) {
    message("M_311 term different");
    return 0;
  }
  if (fabsf(ma->M_320 + mb->M_320) > 0.01f * order5_norm &&
      fabsf(ma->M_320 - mb->M_320) / fabsf(ma->M_320 + mb->M_320) > tolerance) {
    message("M_320 term different");
    return 0;
  }
  if (fabsf(ma->M_401 + mb->M_401) > 0.01f * order5_norm &&
      fabsf(ma->M_401 - mb->M_401) / fabsf(ma->M_401 + mb->M_401) > tolerance) {
    message("M_401 term different");
    return 0;
  }
  if (fabsf(ma->M_410 + mb->M_410) > 0.01f * order5_norm &&
      fabsf(ma->M_410 - mb->M_410) / fabsf(ma->M_410 + mb->M_410) > tolerance) {
    message("M_410 term different");
    return 0;
  }
  if (fabsf(ma->M_500 + mb->M_500) > 0.01f * order5_norm &&
      fabsf(ma->M_500 - mb->M_500) / fabsf(ma->M_500 + mb->M_500) > tolerance) {
    message("M_500 term different");
    return 0;
  }
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
986
987
#endif

988
989
990
991
992
993
994
995
996
997
  /* All is good */
  return 1;
}

/**
 * @brief Constructs the #multipole of a bunch of particles around their
 * centre of mass.
 *
 * Corresponds to equation (28c).
 *
998
 * @param multi The #multipole (content will  be overwritten).
999
1000
1001
 * @param gparts The #gpart.
 * @param gcount The number of particles.
 */
1002
INLINE static void gravity_P2M(struct gravity_tensors *multi,
1003
                               const struct gpart *gparts, int gcount) {
Pedro Gonnet's avatar
Pedro Gonnet committed
1004

1005
1006
1007
  /* Temporary variables */
  double mass = 0.0;
  double com[3] = {0.0, 0.0, 0.0};
1008
  double vel[3] = {0.f, 0.f, 0.f};
1009

1010
  /* Collect the particle data for CoM. */
1011
  for (int k = 0; k < gcount; k++) {
1012
    const double m = gparts[k].mass;
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022

    mass += m;
    com[0] += gparts[k].x[0] * m;
    com[1] += gparts[k].x[1] * m;
    com[2] += gparts[k].x[2] * m;
    vel[0] += gparts[k].v_full[0] * m;
    vel[1] += gparts[k].v_full[1] * m;
    vel[2] += gparts[k].v_full[2] * m;
  }

1023
  /* Final operation on CoM */
1024
  const double imass = 1.0 / mass;
1025
1026
1027
1028
1029
1030
1031
  com[0] *= imass;
  com[1] *= imass;
  com[2] *= imass;
  vel[0] *= imass;
  vel[1] *= imass;
  vel[2] *= imass;

1032
1033
  /* Prepare some local counters */
  double r_max2 = 0.;
1034
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
1035
  double M_100 = 0., M_010 = 0., M_001 = 0.;
1036
1037
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
1038
1039
  double M_200 = 0., M_020 = 0., M_002 = 0.;
  double M_110 = 0., M_101 = 0., M_011 = 0.;
1040
1041
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
1042
1043
1044
1045
  double M_300 = 0., M_030 = 0., M_003 = 0.;
  double M_210 = 0., M_201 = 0., M_120 = 0.;
  double M_021 = 0., M_102 = 0., M_012 = 0.;
  double M_111 = 0.;
1046
1047
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
1048
1049
1050
1051
1052
  double M_400 = 0., M_040 = 0., M_004 = 0.;
  double M_310 = 0., M_301 = 0., M_130 = 0.;
  double M_031 = 0., M_103 = 0., M_013 = 0.;
  double M_220 = 0., M_202 = 0., M_022 = 0.;
  double M_211 = 0., M_121 = 0., M_112 = 0.;
1053
1054
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
1055
1056
1057
1058
1059
1060
1061
  double M_005 = 0., M_014 = 0., M_023 = 0.;
  double M_032 = 0., M_041 = 0., M_050 = 0.;
  double M_104 = 0., M_113 = 0., M_122 = 0.;
  double M_131 = 0., M_140 = 0., M_203 = 0.;
  double M_212 = 0., M_221 = 0., M_230 = 0.;
  double M_302 = 0., M_311 = 0., M_320 = 0.;
  double M_401 = 0., M_410 = 0., M_500 = 0.;
1062
1063
1064
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
1065
1066
1067
1068
#endif

  /* Construce the higher order terms */
  for (int k = 0; k < gcount; k++) {
1069

1070
1071
1072
    const double dx[3] = {gparts[k].x[0] - com[0], gparts[k].x[1] - com[1],
                          gparts[k].x[2] - com[2]};

1073
1074
1075
1076
1077
1078
    /* Maximal distance CoM<->gpart */
    r_max2 = max(r_max2, dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);

#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
    const double m = gparts[k].mass;

1079
    /* 1st order terms */
1080
1081
1082
    M_100 += -m * X_100(dx);
    M_010 += -m * X_010(dx);
    M_001 += -m * X_001(dx);
1083
1084
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
1085
1086

    /* 2nd order terms */
1087
1088
1089
1090
1091
1092
    M_200 += m * X_200(dx);
    M_020 += m * X_020(dx);
    M_002 += m * X_002(dx);
    M_110 += m * X_110(dx);
    M_101 += m * X_101(dx);
    M_011 += m * X_011(dx);
1093
1094
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
1095
1096

    /* 3rd order terms */
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
    M_300 += -m * X_300(dx);
    M_030 += -m * X_030(dx);
    M_003 += -m * X_003(dx);
    M_210 += -m * X_210(dx);
    M_201 += -m * X_201(dx);
    M_120 += -m * X_120(dx);
    M_021 += -m * X_021(dx);
    M_102 += -m * X_102(dx);
    M_012 += -m * X_012(dx);
    M_111 += -m * X_111(dx);
1107
1108
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
1109
1110

    /* 4th order terms */
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
    M_400 += m * X_400(dx);
    M_040 += m * X_040(dx);
    M_004 += m * X_004(dx);
    M_310 += m * X_310(dx);
    M_301 += m * X_301(dx);
    M_130 += m * X_130(dx);
    M_031 += m * X_031(dx);
    M_103 += m * X_103(dx);
    M_013 += m * X_013(dx);
    M_220 += m * X_220(dx);
    M_202 += m * X_202(dx);
    M_022 += m * X_022(dx);
    M_211 += m * X_211(dx);
    M_121 += m * X_121(dx);
    M_112 += m * X_112(dx);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
1128

1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
    /* 5th order terms */
    M_005 += -m * X_005(dx);
    M_014 += -m * X_014(dx);
    M_023 += -m * X_023(dx);
    M_032 += -m * X_032(dx);
    M_041 += -m * X_041(dx);
    M_050 += -m * X_050(dx);
    M_104 += -m * X_104(dx);
    M_113 += -m * X_113(dx);
    M_122 += -m * X_122(dx);
    M_131 += -m * X_131(dx);
    M_140 += -m * X_140(dx);
    M_203 += -m * X_203(dx);
    M_212 += -m * X_212(dx);
    M_221 += -m * X_221(dx);
    M_230 += -m * X_230(dx);
    M_302 += -m * X_302(dx);
    M_311 += -m * X_311(dx);
    M_320 += -m * X_320(dx);
    M_401 += -m * X_401(dx);
    M_410 += -m * X_410(dx);
    M_500 += -m * X_500(dx);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
1154
1155
#endif
  }
1156

1157
1158
1159
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
  M_100 = M_010 = M_001 = 0.f; /* Matthieu */
#endif
1160

1161
  /* Store the data on the multipole. */
1162
1163
1164
1165
1166
1167
1168
1169
  multi->m_pole.M_000 = mass;
  multi->r_max = sqrt(r_max2);
  multi->CoM[0] = com[0];
  multi->CoM[1] = com[1];
  multi->CoM[2] = com[2];
  multi->m_pole.vel[0] = vel[0];
  multi->m_pole.vel[1] = vel[1];
  multi->m_pole.vel[2] = vel[2];
1170

1171
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
1172
1173

  /* 1st order terms */
1174
1175
1176
  multi->m_pole.M_100 = M_100;
  multi->m_pole.M_010 = M_010;
  multi->m_pole.M_001 = M_001;
1177
1178
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
1179
1180

  /* 2nd order terms */
1181
1182
1183
1184
1185
1186
  multi->m_pole.M_200 = M_200;
  multi->m_pole.M_020 = M_020;
  multi->m_pole.M_002 = M_002;
  multi->m_pole.M_110 = M_110;
  multi->m_pole.M_101 = M_101;
  multi->m_pole.M_011 = M_011;
1187
1188
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
1189
1190

  /* 3rd order terms */
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
  multi->m_pole.M_300 = M_300;
  multi->m_pole.M_030 = M_030;
  multi->m_pole.M_003 = M_003;
  multi->m_pole.M_210 = M_210;
  multi->m_pole.M_201 = M_201;
  multi->m_pole.M_120 = M_120;
  multi->m_pole.M_021 = M_021;
  multi->m_pole.M_102 = M_102;
  multi->m_pole.M_012 = M_012;
  multi->m_pole.M_111 = M_111;
1201
#endif
1202
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
1203
1204

  /* 4th order terms */
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
  multi->m_pole.M_400 = M_400;
  multi->m_pole.M_040 = M_040;
  multi->m_pole.M_004 = M_004;
  multi->m_pole.M_310 = M_310;
  multi->m_pole.M_301 = M_301;
  multi->m_pole.M_130 = M_130;
  multi->m_pole.M_031 = M_031;
  multi->m_pole.M_103 = M_103;
  multi->m_pole.M_013 = M_013;
  multi->m_pole.M_220 = M_220;
  multi->m_pole.M_202 = M_202;
  multi->m_pole.M_022 = M_022;
  multi->m_pole.M_211 = M_211;
  multi->m_pole.M_121 = M_121;
  multi->m_pole.M_112 = M_112;
1220