1 files changed, 46 insertions, 0 deletions
diff --git a/src/tests.cc b/src/tests.cc
index 42e79b3..de711eb 100644
--- a/src/tests.cc
+++ b/src/tests.cc
@@ -40,6 +40,38 @@ static CCTK_REAL test_integral(int n, CCTK_REAL (*integration_fn) (const CCTK_RE
   return result;
 }
 
+
+static CCTK_REAL test_pi_symmetric_sphere_integral(CCTK_REAL (*integration_fn) (const CCTK_REAL *, int, int, CCTK_REAL, CCTK_REAL))
+{
+  const int n = 100;
+  const int nth = n;
+  const int nph = n;
+  const int array_size=(nth+1)*(nph+1);
+  const CCTK_REAL PI = acos(-1.0);
+
+  CCTK_REAL *f = new CCTK_REAL[array_size];
+
+  const CCTK_REAL dth = PI/nth;
+  const CCTK_REAL dph = 2*PI/nph;
+
+  for (int ith = 0; ith <= nth; ith++)
+  {
+    for (int iph = 0; iph <= nph; iph++)
+    {
+      const int i = Multipole_Index(ith, iph, nth);
+
+      const CCTK_REAL th = ith*dth;
+      const CCTK_REAL ph = iph*dph;
+
+      f[i] = -(cos(ph)*sqrt(5/PI)*pow(cos(th/2.),3)*sin(th/2.));
+    }
+  }
+
+  const CCTK_REAL result = integration_fn(f, nth, nph, dth, dph);
+  delete [] f;
+  return result;
+}
+
 CCTK_REAL integration_convergence_order(CCTK_REAL (*integration_fn) (const CCTK_REAL *, int, int, CCTK_REAL, CCTK_REAL))
 {
   const int n1 = 100;
@@ -61,6 +93,20 @@ void Multipole_TestIntegrationConvergence(CCTK_ARGUMENTS)
   *test_midpoint_convergence_order = integration_convergence_order(&Midpoint2DIntegral);
 }
 
+void Multipole_TestIntegrationSymmetry(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+
+  *test_simpson_pi_symmetry = test_pi_symmetric_sphere_integral(&Simpson2DIntegral);
+  *test_midpoint_pi_symmetry = test_pi_symmetric_sphere_integral(&Midpoint2DIntegral);
+  *test_trapezoidal_pi_symmetry = test_pi_symmetric_sphere_integral(&Trapezoidal2DIntegral);
+  *test_driscollhealy_pi_symmetry = test_pi_symmetric_sphere_integral(&DriscollHealy2DIntegral);
+  printf("Pi symmetry Simpson integral: %.19g\n", *test_simpson_pi_symmetry);
+  printf("Pi symmetry midpoint integral: %.19g\n", *test_midpoint_pi_symmetry);
+  printf("Pi symmetry trapezoidal integral: %.19g\n", *test_trapezoidal_pi_symmetry);
+  printf("Pi symmetry Driscoll and Healy integral: %.19g\n", *test_driscollhealy_pi_symmetry);
+}
+
 // void Multipole_TestIntegrate(CCTK_ARGUMENTS)
 // {
 //   const int n = 100;