From 86fa90915444a7c4b9571a27978f5e11724ba2a0 Mon Sep 17 00:00:00 2001
From: hinder <hinder@4f5cb9a8-4dd8-4c2d-9bbd-173fa4467843>
Date: Wed, 9 Oct 2013 17:56:00 +0000
Subject: Add convergence test for Simpson integration method

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Multipole/trunk@94 4f5cb9a8-4dd8-4c2d-9bbd-173fa4467843
---
 src/make.code.defn |  2 +-
 src/tests.cc       | 97 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 98 insertions(+), 1 deletion(-)
 create mode 100644 src/tests.cc

(limited to 'src')

diff --git a/src/make.code.defn b/src/make.code.defn
index cccefbd..4d49656 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -1 +1 @@
-SRCS=interpolate.cc multipole.cc utils.cc sphericalharmonic.cc integrate.cc
+SRCS=interpolate.cc multipole.cc utils.cc sphericalharmonic.cc integrate.cc tests.cc
diff --git a/src/tests.cc b/src/tests.cc
new file mode 100644
index 0000000..b63e18b
--- /dev/null
+++ b/src/tests.cc
@@ -0,0 +1,97 @@
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <assert.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "utils.hh"
+#include "integrate.hh"
+
+static CCTK_REAL test_simpson_integral(int n)
+{
+  const int nx = n;
+  const int ny = n;
+  const int array_size=(nx+1)*(ny+1);
+
+  CCTK_REAL *f = new CCTK_REAL[array_size];
+
+  const CCTK_REAL dx = 1./nx;
+  const CCTK_REAL dy = 1./ny;
+
+  for (int ix = 0; ix <= nx; ix++)
+  {
+    for (int iy = 0; iy <= ny; iy++)
+    {
+      const int i = Multipole_Index(ix, iy, nx);
+
+      const CCTK_REAL x = ix*dx;
+      const CCTK_REAL y = iy*dy;
+      const CCTK_REAL PI = acos(-1.0);
+
+      f[i] = x*pow(y,2)*pow(cos(2*PI*y),2)*pow(sin(2*PI*x),2);
+    }
+  }
+
+  const CCTK_REAL result = Simpson2DIntegral(f, nx, ny, dx, dy);
+  delete [] f;
+  return result;
+}
+
+
+void Multipole_TestSimpson(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+
+  const int n1 = 100;
+  const int n2 = 200;
+  const CCTK_REAL PI = acos(-1.0);
+  const CCTK_REAL result1 = test_simpson_integral(100);
+  const CCTK_REAL result2 = test_simpson_integral(200);
+  const CCTK_REAL exact = 1./24 + 1./(64 * pow(PI,2));
+  const CCTK_REAL error1 = fabs(result1 - exact);
+  const CCTK_REAL error2 = fabs(result2 - exact);
+  *test_simpson_convergence_order = log10(error1/error2) / log10((CCTK_REAL) n2/n1);
+}
+
+// void Multipole_TestIntegrate(CCTK_ARGUMENTS)
+// {
+//   const int n = 100;
+
+//   const int nth = n;
+//   const int nph = n;
+//   const int array_size=(nth+1)*(nph+1);
+
+//   CCTK_REAL *f = new CCTK_REAL[array_size];
+
+//   const CCTK_REAL dth = 1/nx;
+//   const CCTK_REAL dph = 1/ny;
+
+//   for (int ix = 0; ix <= nx; ix++)
+//   {
+//     for (int iy = 0; iy <= ny; iy++)
+//     {
+//       const int i = Multipole_Index(ix, iy, nx);
+
+//       const CCTK_REAL x = ix*dx;
+//       const CCTK_REAL y = iy*dy;
+
+//       f[i] = sin(2*PI*x)*cos(2*PI*y);
+//     }
+//   }
+
+
+
+
+//   CCTK_REAL result = Multipole_Integrate(int array_size, int nthetap,
+//                                          CCTK_REAL const array1r[], CCTK_REAL const array1i[],
+//                                          CCTK_REAL const array2r[], CCTK_REAL const array2i[],
+//                                          CCTK_REAL const th[], CCTK_REAL const ph[], 
+//                                          CCTK_REAL *outre, CCTK_REAL *outim)
+
+
+
+//   printf("Integration result: %.19g\n", result);
+// }
-- 
cgit v1.2.3