From 63c9c726455f85651610fdb50982c5737dc65645 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Tue, 21 Apr 2020 23:07:02 +0200
Subject: [PATCH] Added generation of the multipole power calculation to the
 python script

---
 src/multipole.h                               | 98 +++++++++++++++++++
 src/multipole_struct.h                        |  5 +
 .../generate_multipoles/multipoles.py         | 36 ++++++-
 3 files changed, 138 insertions(+), 1 deletion(-)

diff --git a/src/multipole.h b/src/multipole.h
index 1ac413111f..aa8c7f84a1 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -872,6 +872,104 @@ __attribute__((nonnull)) INLINE static int gravity_multipole_equal(
   return 1;
 }
 
+/**
+ * @brief Compute the multipole power of a #multipole.
+ *
+ * @param m The #multipole.
+ */
+__attribute__((nonnull)) INLINE static void gravity_multipole_compute_power(
+    struct multipole *m) {
+
+  double power[SELF_GRAVITY_MULTIPOLE_ORDER + 1] = {0.};
+
+  /* 0th order terms */
+  m->power[0] = m->M_000;
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+  /* 1st order terms */
+  power[1] += m->M_001 * m->M_001;
+  power[1] += m->M_010 * m->M_010;
+  power[1] += m->M_100 * m->M_100;
+
+  m->power[1] = sqrt(power[1]);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+  /* 2nd order terms */
+  power[2] += m->M_002 * m->M_002;
+  power[2] += 5.000000000000000e-01 * m->M_011 * m->M_011;
+  power[2] += m->M_020 * m->M_020;
+  power[2] += 5.000000000000000e-01 * m->M_101 * m->M_101;
+  power[2] += 5.000000000000000e-01 * m->M_110 * m->M_110;
+  power[2] += m->M_200 * m->M_200;
+
+  m->power[2] = sqrt(power[2]);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+  /* 3rd order terms */
+  power[3] += m->M_003 * m->M_003;
+  power[3] += 3.333333333333333e-01 * m->M_012 * m->M_012;
+  power[3] += 3.333333333333333e-01 * m->M_021 * m->M_021;
+  power[3] += m->M_030 * m->M_030;
+  power[3] += 3.333333333333333e-01 * m->M_102 * m->M_102;
+  power[3] += 1.666666666666667e-01 * m->M_111 * m->M_111;
+  power[3] += 3.333333333333333e-01 * m->M_120 * m->M_120;
+  power[3] += 3.333333333333333e-01 * m->M_201 * m->M_201;
+  power[3] += 3.333333333333333e-01 * m->M_210 * m->M_210;
+  power[3] += m->M_300 * m->M_300;
+
+  m->power[3] = sqrt(power[3]);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+  /* 4th order terms */
+  power[4] += m->M_004 * m->M_004;
+  power[4] += 2.500000000000000e-01 * m->M_013 * m->M_013;
+  power[4] += 1.666666666666667e-01 * m->M_022 * m->M_022;
+  power[4] += 2.500000000000000e-01 * m->M_031 * m->M_031;
+  power[4] += m->M_040 * m->M_040;
+  power[4] += 2.500000000000000e-01 * m->M_103 * m->M_103;
+  power[4] += 8.333333333333333e-02 * m->M_112 * m->M_112;
+  power[4] += 8.333333333333333e-02 * m->M_121 * m->M_121;
+  power[4] += 2.500000000000000e-01 * m->M_130 * m->M_130;
+  power[4] += 1.666666666666667e-01 * m->M_202 * m->M_202;
+  power[4] += 8.333333333333333e-02 * m->M_211 * m->M_211;
+  power[4] += 1.666666666666667e-01 * m->M_220 * m->M_220;
+  power[4] += 2.500000000000000e-01 * m->M_301 * m->M_301;
+  power[4] += 2.500000000000000e-01 * m->M_310 * m->M_310;
+  power[4] += m->M_400 * m->M_400;
+
+  m->power[4] = sqrt(power[4]);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+  /* 5th order terms */
+  power[5] += m->M_005 * m->M_005;
+  power[5] += 2.000000000000000e-01 * m->M_014 * m->M_014;
+  power[5] += 1.000000000000000e-01 * m->M_023 * m->M_023;
+  power[5] += 1.000000000000000e-01 * m->M_032 * m->M_032;
+  power[5] += 2.000000000000000e-01 * m->M_041 * m->M_041;
+  power[5] += m->M_050 * m->M_050;
+  power[5] += 2.000000000000000e-01 * m->M_104 * m->M_104;
+  power[5] += 5.000000000000000e-02 * m->M_113 * m->M_113;
+  power[5] += 3.333333333333333e-02 * m->M_122 * m->M_122;
+  power[5] += 5.000000000000000e-02 * m->M_131 * m->M_131;
+  power[5] += 2.000000000000000e-01 * m->M_140 * m->M_140;
+  power[5] += 1.000000000000000e-01 * m->M_203 * m->M_203;
+  power[5] += 3.333333333333333e-02 * m->M_212 * m->M_212;
+  power[5] += 3.333333333333333e-02 * m->M_221 * m->M_221;
+  power[5] += 1.000000000000000e-01 * m->M_230 * m->M_230;
+  power[5] += 1.000000000000000e-01 * m->M_302 * m->M_302;
+  power[5] += 5.000000000000000e-02 * m->M_311 * m->M_311;
+  power[5] += 1.000000000000000e-01 * m->M_320 * m->M_320;
+  power[5] += 2.000000000000000e-01 * m->M_401 * m->M_401;
+  power[5] += 2.000000000000000e-01 * m->M_410 * m->M_410;
+  power[5] += m->M_500 * m->M_500;
+
+  m->power[5] = sqrt(power[5]);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
+#error "Missing implementation for order >5"
+#endif
+}
+
 /**
  * @brief Constructs the #multipole of a bunch of particles around their
  * centre of mass.
diff --git a/src/multipole_struct.h b/src/multipole_struct.h
index ee5e525e28..5d1c1ce218 100644
--- a/src/multipole_struct.h
+++ b/src/multipole_struct.h
@@ -121,6 +121,9 @@ struct multipole {
   /*! Maximal co-moving softening of all the #gpart in the mulipole */
   float max_softening;
 
+  /*! Mulipole power for the different orders */
+  float power[SELF_GRAVITY_MULTIPOLE_ORDER + 1];
+
   /* 0th order term */
   float M_000;
 
@@ -128,12 +131,14 @@ struct multipole {
 
   /* 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
 
diff --git a/theory/Multipoles/generate_multipoles/multipoles.py b/theory/Multipoles/generate_multipoles/multipoles.py
index ef263d09f2..0eca7e25c0 100644
--- a/theory/Multipoles/generate_multipoles/multipoles.py
+++ b/theory/Multipoles/generate_multipoles/multipoles.py
@@ -95,6 +95,7 @@ for i in range(order+1):
         for k in range(order+1):
             if i + j + k == order:
                 print "ma->M_%d%d%d += mb->M_%d%d%d;"%(i,j,k,i,j,k)
+                
 if order > 0:
     print "#endif"
 
@@ -140,6 +141,36 @@ if order > 0:
 print ""
 print "-------------------------------------------------"
 
+print "gravity_multipole_compute_power():"
+print "-------------------------------------------------\n"
+
+if order > 0:
+    print "#if SELF_GRAVITY_MULTIPOLE_ORDER > %d"%(order-1)
+
+print "/* %s order terms */"%ordinal(order)
+
+# Add the terms to the multipole power
+for i in range(order+1):
+    for j in range(order+1):
+        for k in range(order+1):
+            if i + j + k == order:
+                fact1 = factorial(i) * factorial(j) * factorial(k)
+                fact2 = float(factorial(i + j + k))
+                frac = fact1 / fact2
+                if frac == 1.0:
+                    print "power[%d] += m->M_%d%d%d * m->M_%d%d%d;"%(order, i, j, k, i, j, k)
+                else:
+                    print "power[%d] += %12.15e * m->M_%d%d%d * m->M_%d%d%d;"%(order, frac, i, j, k, i, j, k)
+
+print ""
+print "m->power[%d] = sqrt(power[%d]);"%(order, order)
+                    
+if order > 0:
+    print "#endif"
+                    
+print ""
+print "-------------------------------------------------"
+
 print "gravity_P2M(): (loop)"
 print "-------------------------------------------------\n"
 
@@ -157,12 +188,13 @@ for i in range(order+1):
                     print "M_%d%d%d += m * X_%d%d%d(dx);"%(i,j,k,i,j,k)
                 else:
                     print "M_%d%d%d += -m * X_%d%d%d(dx);"%(i,j,k,i,j,k)
+                    
 if order > 0:
     print "#endif"
 
 print ""
 print "-------------------------------------------------"
-    
+
 print "gravity_P2M(): (storing)"
 print "-------------------------------------------------\n"
 
@@ -177,6 +209,7 @@ for i in range(order+1):
         for k in range(order+1):
             if i + j + k == order:
                 print "m->m_pole.M_%d%d%d = M_%d%d%d;"%(i,j,k,i,j,k)
+                
 if order > 0:
     print "#endif"
 
@@ -212,6 +245,7 @@ for i in range(order+1):
 
                                         
                 print ";"
+                
 if order > 0:
     print "#endif"
 
-- 
GitLab