From f724e5d083c3371b001ce9e9e2e568be72abd04d Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Sun, 2 Apr 2017 14:50:48 +0800
Subject: [PATCH] Added script to generate the multipole_equal() function of
any order
---
theory/Multipoles/multipoles.py | 39 +++++++++++++++++++++++++++++++++
1 file changed, 39 insertions(+)
diff --git a/theory/Multipoles/multipoles.py b/theory/Multipoles/multipoles.py
index 06a886e70a..4f665949e5 100644
--- a/theory/Multipoles/multipoles.py
+++ b/theory/Multipoles/multipoles.py
@@ -101,6 +101,45 @@ if order > 0:
print ""
print "-------------------------------------------------"
+print "gravity_multipole_equal():"
+print "-------------------------------------------------\n"
+
+if order > 0:
+ print "#if SELF_GRAVITY_MULTIPOLE_ORDER > %d"%(order-1)
+
+# Create all the terms relevent for this order
+print "/* Manhattan Norm of %s order terms */"%ordinal(order)
+print "const float order%d_norm = "%order,
+first = True
+for i in range(order+1):
+ for j in range(order+1):
+ for k in range(order+1):
+ if i + j + k == order:
+ if first:
+ first = False
+ else:
+ print "+",
+ print "fabsf(ma->M_%d%d%d)"%(i,j,k),
+ print "+ fabsf(mb->M_%d%d%d)"%(i,j,k),
+print ";\n"
+print "/* Compare %s order terms above 1%% of norm */"%ordinal(order)
+for i in range(order+1):
+ for j in range(order+1):
+ for k in range(order+1):
+ if i + j + k == order:
+ print "if (fabsf(ma->M_%d%d%d + mb->M_%d%d%d) > order%d_norm &&"%(i,j,k,i,j,k,order)
+ print " fabsf(ma->M_%d%d%d - mb->M_%d%d%d) / fabsf(ma->M_%d%d%d + mb->M_%d%d%d) > tolerance) {"%(i,j,k,i,j,k,i,j,k,i,j,k)
+ print " message(\"M_%d%d%d term different\");"%(i,j,k)
+ print " return 0;"
+ print "}"
+
+if order > 0:
+ print "#endif"
+
+
+print ""
+print "-------------------------------------------------"
+
print "gravity_P2M(): (loop)"
print "-------------------------------------------------\n"
--
GitLab