use trunc structure

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Multipole/trunk@53 4f5cb9a8-4dd8-4c2d-9bbd-173fa4467843

14 files changed, 1202 insertions, 0 deletions

diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..463a8ce
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,10 @@
+#interface.ccl for thorn Multipole
+implements: multipole
+inherits:   Grid 
+
+public:
+
+CCTK_REAL harmonics type=GF timelevels=1
+{
+  harmonic_re, harmonic_im
+} "Spherical harmonics"
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..e752c8d
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,111 @@
+#param.ccl for thorn Multipole
+
+restricted:
+#Interpolator params:
+CCTK_STRING interpolator_name "Which interpolator should I use"
+{
+    ".+" :: "Any nonempty string"
+} "Hermite polynomial interpolation"
+
+CCTK_STRING interpolator_pars "Parameters for the interpolator"
+{
+    ".*" :: "Any string that Util_TableSetFromString() will take"
+} "order=3"
+
+CCTK_STRING coord_system "What is the coord system?"
+{
+    ".*" :: "Any smart string will do"
+} "cart3d"
+
+KEYWORD integration_method "How to do surface integrals" STEERABLE=always
+{
+  "midpoint" :: "Midpoint rule (1st order)"
+  "Simpson"  :: "Simpson's rule (4th order)"
+} "midpoint"
+
+CCTK_INT out_every "How often to output" \
+STEERABLE=recover
+{
+  * :: "Any integer"
+} 1
+
+CCTK_INT out_1d_every "How often to output 1d data" \
+STEERABLE=recover
+{
+  * :: "Any integer"
+} 0
+
+#physical params:
+CCTK_INT nradii "How many extraction radii?" \
+STEERABLE=recover
+{
+  0:101 :: "A positive integer less than 101"
+} 1
+
+
+CCTK_REAL radius[101] "The radii for extraction" \
+STEERABLE=recover
+{
+  0.0:* :: "Please keep it in the grid"
+} 0.0
+
+CCTK_INT ntheta "How many points in the theta direction?" \
+STEERABLE=recover
+{
+  0:* :: "Positive please"
+} 50
+
+CCTK_INT nphi "How many points in the phi direction?" \
+STEERABLE=recover
+{
+  0:* :: "Positive please"
+} 100
+
+CCTK_STRING variables "What variables to decompose"
+{
+    ".*" :: "A list of variables"
+} ""
+
+CCTK_BOOLEAN verbose "Output detailed information about what is happening"
+{
+} "no"
+
+KEYWORD mode_type "Which type of mode extraction do we have" STEERABLE=always
+{
+  "all modes"     :: "Extract all modes up to (l_mode, m_mode)."
+  "specific mode" :: "Select one specific (l_mode, m_mode) mode"
+} "all modes"
+
+CCTK_INT l_min "all modes: above which l mode to calculate/ specific mode: which l mode to extract" STEERABLE=always
+{
+  0:* :: "l >= 0"
+} 2
+
+CCTK_INT l_mode "all modes: Up to which l mode to calculate/ specific mode: which l mode to extract" STEERABLE=always
+{
+  0:* :: "l >= 0"
+} 2
+                                                                                                         
+CCTK_INT m_mode "all modes: Up to which m mode to calculate/ specific mode: which m mode to extract " STEERABLE=always
+{
+  0:* :: "Positive Please"
+} 2
+
+CCTK_BOOLEAN enable_test "whether to set a spherical harmonic in the 'harmonic' grid functions"
+{
+} "no"
+
+CCTK_INT test_l "which mode to put into the test variables"
+{
+  * :: "Any integer"
+} 2
+
+CCTK_INT test_m "which mode to put into the test variables"
+{
+  * :: "Any integer"
+} 2
+
+CCTK_INT test_sw "which spin weight to put into the test variables"
+{
+  * :: "Any integer"
+} -2
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..a7123ee
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,20 @@
+#schedule.ccl for thorn Multipole
+
+if (enable_test)
+{
+  STORAGE: harmonics[1]
+}
+
+schedule Multipole_Calc at CCTK_ANALYSIS after (calc_np,PsiKadelia)
+{
+  LANG: C
+  OPTIONS: GLOBAL
+} "Calculate Multipoles"
+
+if (enable_test)
+{
+    schedule Multipole_SetHarmonic at CCTK_INITIAL
+    {
+      LANG: C
+    } "Populate grid functions with spherical harmonics"
+}
diff --git a/src/integrate.cc b/src/integrate.cc
new file mode 100644
index 0000000..0f1c392
--- /dev/null
+++ b/src/integrate.cc
@@ -0,0 +1,152 @@
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <assert.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+/*
+
+We will want to integrate functions F(th,ph) from th = 0 to pi, ph = 0
+to 2 pi with a weighting function sin(th).  Alternatively, we might
+want to use u = cos(th) as the variable, in which case we will go from
+u = -1 to 1 and ph = 0 to 2 pi.  For simplicity, we implement an
+integration routine with a weight function of 1, and require the user
+to multiply the integrand by their own weight function.  We divide the
+interval [a,b] into nx subintervals of spacing h = (b-a)/nx.  These
+have coordinates [x_i-1, xi] where x_i = x_0 + i h. So i runs from 0
+to nx.  We require the function to integrate at the points F[x_i,
+y_i].  We have x_0 = a and x_n = b.  Check: x_n = x_0 + n (b-a)/n = a
++ b - a = b.  Good.  
+
+If we are given these points in an array, we also need the width and
+height of the array.  To get an actual integral, we also need the grid
+spacing hx and hy, but these are just multiplied by the result to give
+the integral.
+
+*/
+
+
+
+#define idx(xx,yy) (assert((xx) <= nx), assert((xx) >= 0), assert((yy) <= ny), assert((yy) >= 0), ((xx) + (yy) * (nx+1)))
+
+// Hard coded 2D integrals
+
+static CCTK_REAL Trapezoidal2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, CCTK_REAL hy)
+{
+  CCTK_REAL integrand_sum = 0.0;
+  int ix = 0, iy = 0;
+
+  assert(nx > 0); assert(ny > 0); assert (f);
+
+  // Corners
+  integrand_sum += f[idx(0,0)] + f[idx(nx,0)] + f[idx(0,ny)] + f[idx(nx,ny)];
+
+  // Edges
+  for (ix = 1; ix <= nx-1; ix++)
+    integrand_sum += 2 * f[idx(ix,0)] + 2 * f[idx(ix,ny)];
+
+  for (iy = 1; iy <= ny-1; iy++)
+    integrand_sum += 2 * f[idx(0,iy)] + 2 * f[idx(nx,iy)];
+
+  // Interior
+  for (iy = 1; iy <= ny-1; iy++)
+    for (ix = 1; ix <= nx-1; ix++)
+      integrand_sum += 4 * f[idx(ix,iy)];
+
+  return (double) 1 / (double) 4 * hx * hy * integrand_sum;
+}
+
+CCTK_REAL Simpson2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, CCTK_REAL hy)
+{
+  CCTK_REAL integrand_sum = 0;
+  int ix = 0, iy = 0;
+
+  assert(nx > 0); assert(ny > 0); assert (f);
+  assert(nx % 2 == 0);
+  assert(ny % 2 == 0);
+
+  int px = nx / 2;
+  int py = ny / 2;
+
+  // Corners
+  integrand_sum += f[idx(0,0)] + f[idx(nx,0)] + f[idx(0,ny)] + f[idx(nx,ny)];
+
+  // Edges
+  for (iy = 1; iy <= py; iy++)
+    integrand_sum += 4 * f[idx(0,2*iy-1)] + 4 * f[idx(nx,2*iy-1)];
+
+  for (iy = 1; iy <= py-1; iy++)
+    integrand_sum += 2 * f[idx(0,2*iy)] + 2 * f[idx(nx,2*iy)];
+
+  for (ix = 1; ix <= px; ix++)
+    integrand_sum += 4 * f[idx(2*ix-1,0)] + 4 * f[idx(2*ix-1,ny)];
+
+  for (ix = 1; ix <= px-1; ix++)
+    integrand_sum += 2 * f[idx(2*ix,0)] + 2 * f[idx(2*ix,ny)];
+
+  // Interior
+  for (iy = 1; iy <= py; iy++)
+    for (ix = 1; ix <= px; ix++)
+      integrand_sum += 16 * f[idx(2*ix-1,2*iy-1)];
+
+  for (iy = 1; iy <= py-1; iy++)
+    for (ix = 1; ix <= px; ix++)
+      integrand_sum += 8 * f[idx(2*ix-1,2*iy)];
+
+  for (iy = 1; iy <= py; iy++)
+    for (ix = 1; ix <= px-1; ix++)
+      integrand_sum += 8 * f[idx(2*ix,2*iy-1)];
+
+  for (iy = 1; iy <= py-1; iy++)
+    for (ix = 1; ix <= px-1; ix++)
+      integrand_sum += 4 * f[idx(2*ix,2*iy)];
+
+  return ((double) 1 / (double) 9) * hx * hy * integrand_sum;
+}
+
+// 1D integrals
+
+static CCTK_REAL Simpson1DIntegral(CCTK_REAL const *f, int n, CCTK_REAL h)
+{
+  CCTK_REAL integrand_sum = 0;
+  int i = 0;
+
+  assert(n > 0); assert(f);
+  assert(n % 2 == 0);
+
+  int p = n / 2;
+
+  integrand_sum += f[0] + f[n];
+
+  for (i = 1; i <= p-1; i++)
+    integrand_sum += 4 * f[2*i-1] + 2 * f[2*i];
+
+  integrand_sum += 4 * f[2*p-1];
+
+  return 1.0/3.0 * h * integrand_sum;
+}
+
+// 2D integral built up from 1D
+
+static CCTK_REAL Composite2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, CCTK_REAL hy)
+{
+  CCTK_REAL integrand_sum = 0;
+
+  assert(nx > 0); assert(ny > 0); assert (f);
+  assert(nx % 2 == 0);
+  assert(ny % 2 == 0);
+
+  CCTK_REAL *g = new CCTK_REAL[ny+1];
+
+  for (int i = 0; i <= ny; i++)
+  {
+    g[i] = Simpson1DIntegral(&f[idx(0,i)], nx, hx);
+  }
+
+  integrand_sum = Simpson1DIntegral(g, ny, hy);
+  delete [] g;
+  return integrand_sum;
+}
diff --git a/src/integrate.hh b/src/integrate.hh
new file mode 100644
index 0000000..c45917e
--- /dev/null
+++ b/src/integrate.hh
@@ -0,0 +1,9 @@
+
+#ifndef __integrate_h
+#define __integrate_h
+#include "cctk.h"
+
+CCTK_REAL Simpson2DIntegral(CCTK_REAL const *f, int nx, int ny, 
+                            CCTK_REAL hx, CCTK_REAL hy);
+
+#endif
diff --git a/src/interpolate.cc b/src/interpolate.cc
new file mode 100644
index 0000000..5c9946c
--- /dev/null
+++ b/src/interpolate.cc
@@ -0,0 +1,146 @@
+#include <stdio.h>
+
+#include "interpolate.hh"
+
+static
+void report_interp_error(int ierr)
+{
+  if (ierr < 0)
+  {
+    CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+        "CCTK_InterpGridArrays returned error code %d",ierr);
+  }
+}
+  
+void Multipole_Interp(CCTK_ARGUMENTS,
+                      CCTK_REAL xs[], CCTK_REAL ys[], CCTK_REAL zs[],
+                      int real_idx, int imag_idx,
+                      CCTK_REAL sphere_real[], CCTK_REAL sphere_imag[])
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  // Need parameters for the following:
+  // ntheta (dtheta = pi/(ntheta) 
+  // nphi (dphi = 2pi/(nphi)
+  // r (radius of sphere)
+  // NOTE: depending on the interval of integration, denominator above may
+  // need to be modified to avoid double counting
+
+
+  CCTK_INT num_input_arrays  = imag_idx == -1 ? 1 : 2;
+  CCTK_INT num_output_arrays = imag_idx == -1 ? 1 : 2;
+  const CCTK_INT num_dims = 3;
+  int ierr = -1;
+
+  const void* interp_coords[num_dims] 
+    = { (const void *) xs,
+        (const void *) ys,
+        (const void *) zs };
+
+  const CCTK_INT input_array_indices[2]
+    = { real_idx,
+        imag_idx };
+
+  const CCTK_INT output_array_types[2]
+    = { CCTK_VARIABLE_REAL,
+        CCTK_VARIABLE_REAL };
+
+  void * output_arrays[2]
+    = { (void *) sphere_real,
+        (void *) sphere_imag };
+
+  const int operator_handle = CCTK_InterpHandle(interpolator_name);
+
+  int param_table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT);
+  ierr = Util_TableSetFromString(param_table_handle, interpolator_pars);
+
+  const int coord_system_handle = CCTK_CoordSystemHandle(coord_system);
+
+  ierr = CCTK_InterpGridArrays(
+      cctkGH,
+      num_dims,
+      operator_handle,
+      param_table_handle,
+      coord_system_handle,
+      CCTK_MyProc(cctkGH) == 0 ? (ntheta+1)*(nphi+1) : 0, // Only the 0 processor needs the points
+      CCTK_VARIABLE_REAL,
+      interp_coords,
+      num_input_arrays,
+      input_array_indices,
+      num_output_arrays,
+      output_array_types,
+      output_arrays);
+
+  report_interp_error(ierr);
+
+  if (imag_idx == -1)
+  {
+    for (int i = 0; i < (ntheta+1)*(nphi+1); i++)
+    {
+      sphere_imag[i] = 0;
+    }
+  }
+
+  Util_TableDestroy(param_table_handle);
+}
+
+// // Debugging routine
+// void Multipole_InterpVar(CCTK_ARGUMENTS,
+//                          CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[], const char *var_name,
+//                          CCTK_REAL sphere_var[])
+// {
+//   DECLARE_CCTK_ARGUMENTS;
+//   DECLARE_CCTK_PARAMETERS;
+
+//   // Need parameters for the following:
+//   // ntheta (dtheta = pi/(ntheta) 
+//   // nphi (dphi = 2pi/(nphi)
+//   // r (radius of sphere)
+//   // NOTE: depending on the interval of integration, denominator above may
+//   // need to be modified to avoid double counting
+
+// //  printf("Multipole: Interpolating %s\n", var_name);
+
+//   const CCTK_INT num_input_arrays = 1;
+//   const CCTK_INT num_output_arrays = 1;
+//   const CCTK_INT num_dims = 3;
+//   int ierr = -1;
+
+//   const void* interp_coords[num_dims] 
+//     = { (const void *) x,
+//         (const void *) y,
+//         (const void *) z };
+
+//   const CCTK_INT input_array_indices[num_input_arrays]
+//     = { CCTK_VarIndex(var_name) };
+
+//   const CCTK_INT output_array_types[num_output_arrays]
+//     = { CCTK_VARIABLE_REAL };
+
+//   void * output_arrays[num_output_arrays]
+//     = { (void *) sphere_var };
+
+//   const int operator_handle = CCTK_InterpHandle(interpolator_name);
+
+//   int param_table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT);
+//   ierr = Util_TableSetFromString(param_table_handle, interpolator_pars);
+
+//   const int coord_system_handle = CCTK_CoordSystemHandle(coord_system);
+
+//   ierr = CCTK_InterpGridArrays(
+//       cctkGH,
+//       num_dims,
+//       operator_handle,
+//       param_table_handle,
+//       coord_system_handle,
+//       (ntheta+1)*(nphi+1),
+//       CCTK_VARIABLE_REAL,
+//       interp_coords,
+//       num_input_arrays,
+//       input_array_indices,
+//       num_output_arrays,
+//       output_array_types,
+//       output_arrays);
+//   report_interp_error(ierr);
+// }
diff --git a/src/interpolate.hh b/src/interpolate.hh
new file mode 100644
index 0000000..6668abc
--- /dev/null
+++ b/src/interpolate.hh
@@ -0,0 +1,19 @@
+#include <math.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "cctk_Functions.h"
+#include "util_Table.h"
+
+// Multipole_Interp:
+//      This function interpolates psi4 onto the sphere in cartesian 
+//      coordinates as created by Multipole_CoordSetup.
+void Multipole_Interp(CCTK_ARGUMENTS,
+                      CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[], 
+                      int real_idx, int imag_idx,
+                      CCTK_REAL psi4r[], CCTK_REAL psi4i[]);
+
+void Multipole_InterpVar(CCTK_ARGUMENTS,
+                         CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[], const char *var_name,
+                         CCTK_REAL sphere_var[]);
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..cccefbd
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1 @@
+SRCS=interpolate.cc multipole.cc utils.cc sphericalharmonic.cc integrate.cc
diff --git a/src/multipole.cc b/src/multipole.cc
new file mode 100644
index 0000000..30a6624
--- /dev/null
+++ b/src/multipole.cc
@@ -0,0 +1,320 @@
+#include <stdio.h>
+#include <string.h>
+#include <assert.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "cctk_Functions.h"
+
+#include "multipole.hh"
+#include "interpolate.hh"
+#include "utils.hh"
+#include "sphericalharmonic.hh"
+
+static const int var_name_length = 30;
+static const int max_vars = 10;
+
+static const int max_spin_weights = 3;
+static const int max_l_modes = 10;
+static const int max_m_modes = 2 * max_l_modes + 1;
+
+typedef struct
+{
+  int index;
+  int imag_index;
+  int spin_weight;
+  char name[var_name_length];
+}
+variable_desc;
+
+typedef struct
+{
+  int n_vars;
+  variable_desc *vars;
+}
+variables_desc;
+
+
+static void fill_variable(int idx, const char *optstring, void *callback_arg)
+{
+  assert(idx >= 0);
+  assert(callback_arg != 0);
+
+  variables_desc *vs = (variables_desc * ) callback_arg;
+
+  assert(vs->n_vars < max_vars); // Too many variables in the variables list
+  variable_desc *v = &vs->vars[vs->n_vars];
+
+  v->index = idx;
+
+  // Default values if there is no option string or if the options are
+  // not present
+  v->imag_index = -1;
+  v->spin_weight = 0;
+  strcpy(v->name, CCTK_VarName(v->index));
+
+  if (optstring != 0)
+  {
+    int table = Util_TableCreateFromString(optstring);
+
+    if (table >= 0)
+    {
+      const int buffer_length = 256;
+      char imag_name[buffer_length];
+      
+      Util_TableGetInt(table, &v->spin_weight , "sw");
+/////////////////////////////////////////////////////////////
+CCTK_VInfo(CCTK_THORNSTRING,"spinweight %d", v->spin_weight);
+/////////////////////////////////////////////////////////////
+      if (Util_TableGetString(table, buffer_length, imag_name , "cmplx") >= 0)
+      {
+        v->imag_index = CCTK_VarIndex(imag_name);
+      }
+      Util_TableGetString(table, var_name_length, v->name , "name");
+    }
+  }
+  vs->n_vars++;
+}
+
+static void parse_variables_string(const char *var_string, variable_desc v[max_vars], int *n_variables)
+{
+  variables_desc  vars;
+
+  vars.n_vars = 0;
+  vars.vars = v;
+
+  const char *actual_var_string = strlen(var_string) == 0 ? 
+    "WeylScal4::Psi4r{sw=-2 cmplx=\'WeylScal4::Psi4i\' name=\'psi4\'}" : 
+    var_string;
+
+  if (strlen(var_string) == 0)
+  {
+    static bool have_warned = false;
+    if (!have_warned)
+    {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "WARNING: No Multipole::variables string has been specified.  For compatibility, the current default is to use \"%s\" but in future, this default will change to \"\".  Please update your parameter files to include this new parameter.", actual_var_string);
+      have_warned = true;
+    }
+  }
+    
+  int ierr = CCTK_TraverseString(actual_var_string, fill_variable, &vars, CCTK_GROUP_OR_VAR);
+  assert(ierr > 0);
+
+  *n_variables = vars.n_vars;
+}
+
+static void output_mode(CCTK_ARGUMENTS, const variable_desc *v, CCTK_REAL rad, 
+                     CCTK_REAL real_lm, CCTK_REAL imag_lm, int l, int m)
+{
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  char name_tmp[1000];
+
+  if (CCTK_MyProc(cctkGH) == 0)
+  {
+    sprintf(name_tmp, "mp_%s_l%d_m%d_r%1.2f.asc", v->name, l, m, rad);
+    Multipole_OutputComplexToFile(CCTK_PASS_CTOC, name_tmp, real_lm, imag_lm);
+  }
+}
+
+static void output_1D(CCTK_ARGUMENTS, const variable_desc *v, CCTK_REAL rad, 
+                      CCTK_REAL *th, CCTK_REAL *ph,
+                      CCTK_REAL *real, CCTK_REAL *imag, int array_size)
+{
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  char name_tmp[1000];
+
+  if (CCTK_MyProc(cctkGH) == 0)
+  {
+    if (out_1d_every != 0 && (cctk_iteration) % out_1d_every == 0)
+    {
+      const char *real_name = CCTK_VarName(v->index);
+      const char *imag_name = CCTK_VarName(v->imag_index);
+
+      sprintf(name_tmp, "mp_%s_r%1.2f.th.asc", real_name, rad);
+      Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, array_size, th, ph, mp_theta, real);
+      sprintf(name_tmp, "mp_%s_r%1.2f.th.asc", imag_name, rad);
+      Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, array_size, th, ph, mp_theta, imag);
+
+      sprintf(name_tmp, "mp_%s_r%1.2f.ph.asc", real_name, rad);
+      Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, array_size, th, ph, mp_phi, real);
+      sprintf(name_tmp, "mp_%s_r%1.2f.ph.asc", imag_name, rad);
+      Multipole_Output1D(CCTK_PASS_CTOC, name_tmp, array_size, th, ph, mp_phi, imag);
+    }
+  }
+}
+
+bool int_in_array(int a, const int array[], int len)
+{
+  for (int i = 0; i < len; i++)
+  {
+    if (array[i] == a)
+      return true;
+  }
+  return false;
+}
+
+int find_int_in_array(int a, const int array[], int len)
+{
+  for (int i = 0; i < len; i++)
+  {
+    if (array[i] == a)
+      return i;
+  }
+  return -1;
+}
+
+
+static
+void get_spin_weights(variable_desc vars[], int n_vars, int spin_weights[max_spin_weights], int *n_weights)
+{
+  int n_spin_weights = 0;
+
+  for (int i = 0; i < n_vars; i++)
+  {
+    if (!int_in_array(vars[i].spin_weight, spin_weights, n_spin_weights))
+    {
+      assert(n_spin_weights < max_spin_weights);
+      spin_weights[n_spin_weights] = vars[i].spin_weight;
+      n_spin_weights++;
+    }
+  }
+  *n_weights = n_spin_weights;
+}
+
+
+static
+void setup_harmonics(const int spin_weights[max_spin_weights], int n_spin_weights, int lmin, int lmax, 
+                     int m_mode, bool allmodes, CCTK_REAL th[], CCTK_REAL ph[], int array_size, 
+                     CCTK_REAL *reY[max_spin_weights][max_l_modes][max_m_modes],
+                     CCTK_REAL *imY[max_spin_weights][max_l_modes][max_m_modes])
+{
+  for (int si = 0; si < max_spin_weights; si++)
+  {
+    for (int l = 0; l < max_l_modes; l++)
+    {
+      for (int mi = 0; mi < max_m_modes; mi++)
+      {
+        reY[si][l][mi] = 0;
+        imY[si][l][mi] = 0;
+      }
+    }
+  }
+  for (int si = 0; si < n_spin_weights; si++)
+  {
+    int sw = spin_weights[si];
+
+    for (int l=lmin; l <= lmax; l++)
+    {
+      int mmin = (allmodes ? (m_mode > l ? -l : -m_mode) : m_mode);
+      int mmax = (allmodes ? (m_mode > l ? l : m_mode) : m_mode);
+
+      for (int m=mmin; m <= mmax; m++)
+      {
+        reY[si][l-lmin][m-mmin] = new CCTK_REAL[array_size];
+        imY[si][l-lmin][m-mmin] = new CCTK_REAL[array_size];
+
+        Multipole_HarmonicSetup(sw, l, m, array_size, th, ph, reY[si][l-lmin][m-mmin], imY[si][l-lmin][m-mmin]);
+      }
+    }
+  } 
+}
+
+extern "C" void Multipole_Calc(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  static CCTK_REAL *xs, *ys, *zs;
+  static CCTK_REAL *xhat, *yhat, *zhat;
+  static CCTK_REAL *th, *ph;
+  static CCTK_REAL *real = 0, *imag = 0;
+  static CCTK_REAL *reY[max_spin_weights][max_l_modes][max_m_modes];
+  static CCTK_REAL *imY[max_spin_weights][max_l_modes][max_m_modes];
+  static variable_desc   vars[max_vars];
+  static int n_variables = 0;
+  static int spin_weights[max_spin_weights];
+  static int n_spin_weights = 0;
+
+  static bool initialized = false;
+
+  int i, array_size=(ntheta+1)*(nphi+1);
+  CCTK_REAL real_lm = 0.0, imag_lm = 0.0;
+  int lmin, lmax, mmin, mmax, allmodes=0;
+
+  if (cctk_iteration % out_every != 0)
+    return;
+
+  // Determine modes to compute
+  if (CCTK_EQUALS(mode_type, "all modes"))
+  {
+    allmodes=1;
+  }
+  //if mode_type is set to "all modes", all modes with l <= l_mode, 
+  //abs(m) <= min{l_mode, m_mode} will be output.  Otherwise, only
+  //the l_mode, m_mode mode will be output.
+  lmin = (allmodes ? l_min : l_mode);
+  lmax = l_mode;
+
+  assert(lmax + 1 <= max_l_modes);
+
+  if (!initialized)
+  {
+    printf("Multipole: Allocating memory in Multipole_Calc\n");
+    real = new CCTK_REAL[array_size];
+    imag = new CCTK_REAL[array_size];
+    th = new CCTK_REAL[array_size];  
+    ph = new CCTK_REAL[array_size];
+    xs = new CCTK_REAL[array_size];
+    ys = new CCTK_REAL[array_size];
+    zs = new CCTK_REAL[array_size];
+    xhat = new CCTK_REAL[array_size];
+    yhat = new CCTK_REAL[array_size];
+    zhat = new CCTK_REAL[array_size];
+  
+    printf("Multipole: parsing variables string.\n");
+    parse_variables_string(variables, vars, &n_variables);
+    get_spin_weights(vars, n_variables, spin_weights, &n_spin_weights);
+    Multipole_CoordSetup(ntheta, nphi, xhat, yhat, zhat, th, ph);
+    setup_harmonics(spin_weights, n_spin_weights, lmin, lmax, m_mode, allmodes, th, ph, array_size, reY, imY);
+    initialized = true;
+  }
+
+  for (int v = 0; v < n_variables; v++)
+  {
+    //assert(vars[v].spin_weight == -2);
+
+    int si = find_int_in_array(vars[v].spin_weight, spin_weights, n_spin_weights);
+    assert(si != -1);
+
+    for(i=0; i<nradii; i++)
+    {
+      // Compute x^i = r * \hat x^i
+      Multipole_ScaleCartesian(ntheta, nphi, radius[i], xhat, yhat, zhat, xs, ys, zs);
+      
+      // Interpolate Psi4r and Psi4i
+      Multipole_Interp(CCTK_PASS_CTOC, xs, ys, zs, vars[v].index, vars[v].imag_index, 
+                       real, imag);
+      for (int l=lmin; l <= lmax; l++)
+      {
+        mmin = (allmodes ? (m_mode > l ? -l : -m_mode) : m_mode) ;
+        mmax = (allmodes ? (m_mode > l ? l : m_mode) : m_mode);
+        
+        for (int m=mmin; m <= mmax; m++)
+        {
+          // Integrate sYlm (real + i imag) over the sphere at radius r
+          Multipole_Integrate(array_size, ntheta, reY[si][l-lmin][m-mmin], imY[si][l-lmin][m-mmin], real, imag, th, ph,
+                              &real_lm, &imag_lm);
+    
+          output_mode(CCTK_PASS_CTOC, &vars[v], radius[i], real_lm, imag_lm, l, m);
+        }//loop over m
+      }//loop over l
+      output_1D(CCTK_PASS_CTOC, &vars[v], radius[i], th, ph, real, imag, array_size);
+    }//loop over radii
+  }//loop over variables
+}
diff --git a/src/multipole.hh b/src/multipole.hh
new file mode 100644
index 0000000..004d351
--- /dev/null
+++ b/src/multipole.hh
@@ -0,0 +1,15 @@
+#include "cctk.h"
+#include "cctk_Arguments.h"
+
+// Multipole_Calc
+//      This is the main scheduling file.  Because we are completely local here
+//      and do not use cactus arrays etc, we schedule only one function and then
+//      like program like one would in C, C++ with this function taking the 
+//      place of int main(void).
+//
+//      This function calls functions to accomplish 3 things:
+//        1) Interpolate psi4 onto a sphere
+//        2) Integrate psi4 with the ylm's over that sphere
+//        2) Output the mode decomposed psi4
+extern "C" void Multipole_Calc(CCTK_ARGUMENTS);
+
diff --git a/src/sphericalharmonic.cc b/src/sphericalharmonic.cc
new file mode 100644
index 0000000..116613c
--- /dev/null
+++ b/src/sphericalharmonic.cc
@@ -0,0 +1,97 @@
+#include <math.h>
+#include <assert.h>
+#include <iostream>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+static const CCTK_REAL PI = acos(-1.0);
+
+static double factorial(int n)
+{
+  double returnval = 1;
+  for (int i = n; i >= 1; i--)
+  {
+    returnval *= i;
+  }
+  return returnval;
+}
+
+static inline double combination(int n, int m)
+{
+  // Binomial coefficient is undefined if these conditions do not hold
+  assert(n >= 0);
+  assert(m >= 0);
+  assert(m <= n);
+
+  return factorial(n) / (factorial(m) * factorial(n-m));
+}
+
+static inline int imin(int a, int b)
+{
+  return a < b ? a : b;
+}
+
+static inline int imax(int a, int b)
+{
+  return a > b ? a : b;
+}
+
+static
+void Multipole_SphericalHarmonic(int s, int l, int m, 
+                                 CCTK_REAL th, CCTK_REAL ph,
+                                 CCTK_REAL *reY, CCTK_REAL *imY)
+{
+//  assert(s == -2 && l == 2 && m == 2);
+//  *reY = 1.0/2.0 * sqrt(5/PI) * pow(cos(th/2), 4) * cos(2*ph);
+//  *imY = 1.0/2.0 * sqrt(5/PI) * pow(cos(th/2), 4) * sin(2*ph);
+  double all_coeff = 0, sum = 0;
+  all_coeff = pow(-1.0, m);
+  all_coeff *= sqrt(factorial(l+m)*factorial(l-m)*(2*l+1) / (4.*PI*factorial(l+s)*factorial(l-s)));
+  sum = 0.;
+  for (int i = imax(m - s, 0); i <= imin(l + m, l - s); i++)
+  {
+    double sum_coeff = combination(l-s, i) * combination(l+s, i+s-m);
+    sum += sum_coeff * pow(-1.0, l-i-s) * pow(cos(th/2.), 2 * i + s - m) * 
+      pow(sin(th/2.), 2*(l-i)+m-s);
+  }
+  *reY = all_coeff*sum*cos(m*ph);
+  *imY = all_coeff*sum*sin(m*ph);
+}
+
+void Multipole_HarmonicSetup(int s, int l, int m,
+                             int array_size, CCTK_REAL const th[], CCTK_REAL const ph[],
+                             CCTK_REAL reY[], CCTK_REAL imY[])
+{
+  for (int i = 0; i < array_size; i++)
+  {
+    Multipole_SphericalHarmonic(s,l,m,th[i],ph[i],&reY[i], &imY[i]);
+  }
+}
+
+
+// Fill a grid function with a given spherical harmonic
+extern "C" void Multipole_SetHarmonic(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  for (int k = 0; k < cctk_lsh[2]; k++)
+  {
+    for (int j = 0; j < cctk_lsh[1]; j++)
+    {
+      for (int i = 0; i < cctk_lsh[0]; i++)
+      {
+        int index = CCTK_GFINDEX3D(cctkGH,i,j,k) ;
+
+        CCTK_REAL theta = acos(z[index]/r[index]);
+        CCTK_REAL phi = atan2(y[index],x[index]);
+
+        Multipole_SphericalHarmonic(test_sw,test_l,test_m,theta,phi,
+                                    &harmonic_re[index], &harmonic_im[index]);
+      }
+    }
+  }
+  return;
+}
diff --git a/src/sphericalharmonic.hh b/src/sphericalharmonic.hh
new file mode 100644
index 0000000..32b1bba
--- /dev/null
+++ b/src/sphericalharmonic.hh
@@ -0,0 +1,11 @@
+
+#ifndef __sphericalharmonic_h
+#define __sphericalharmonic_h
+#include "cctk.h"
+
+void Multipole_HarmonicSetup(int s, int l, int m,
+                             int array_size, CCTK_REAL const th[], CCTK_REAL const ph[],
+                             CCTK_REAL reY[], CCTK_REAL imY[]);
+
+
+#endif
diff --git a/src/utils.cc b/src/utils.cc
new file mode 100644
index 0000000..00c4752
--- /dev/null
+++ b/src/utils.cc
@@ -0,0 +1,251 @@
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "utils.hh"
+#include "integrate.hh"
+
+////////////////////////////////////////////////////////////////////////
+// File manipulation
+////////////////////////////////////////////////////////////////////////
+
+FILE *Multipole_OpenOutputFile(CCTK_ARGUMENTS, char const *name)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  const int buf_size = 1024;
+  bool first_time = cctk_iteration == 0;
+  const char *mode = first_time ? "w" : "a";
+  char output_name[buf_size];
+  CCTK_STRING *out_dir  = (CCTK_STRING *) CCTK_ParameterGet("out_dir", "IOUtil", NULL);
+
+  if (*out_dir == 0)
+  {
+    CCTK_WARN(1,"Parameter IOUtil::out_dir not found");
+    return 0;
+  }
+
+  if ((int)strlen(*out_dir)+1 > buf_size)
+  {
+    CCTK_WARN(1,"Parameter IOUtil::out_dir is too long");
+    return 0;
+  }
+
+  sprintf(output_name, "%s/%s", *out_dir, name);
+
+  FILE *fp = fopen(output_name, mode);
+
+  if (fp == 0)
+  {
+    CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "Could not open output file %s", output_name);
+  }
+
+  return fp;
+}
+
+////////////////////////////////////////////////////////////////////////
+// Unused
+////////////////////////////////////////////////////////////////////////
+
+void Multipole_OutputArray(CCTK_ARGUMENTS, FILE *f, int array_size, 
+                           CCTK_REAL const th[], CCTK_REAL const ph[], 
+                           CCTK_REAL const xs[], CCTK_REAL const ys[], CCTK_REAL const zs[],
+                           CCTK_REAL const data[])
+{
+  DECLARE_CCTK_PARAMETERS;
+  DECLARE_CCTK_ARGUMENTS;
+
+  CCTK_REAL last_ph = ph[0];
+
+  for (int i = 0; i < array_size; i++)
+  {
+    if (ph[i] != last_ph) // Separate blocks for gnuplot
+      fprintf(f, "\n");
+    fprintf(f, "%f %f %f %f %f %f %.19g\n", cctk_time, th[i], ph[i], xs[i], ys[i], zs[i], data[i]);
+    last_ph = ph[i];
+  }
+}
+
+
+void Multipole_OutputArrayToFile(CCTK_ARGUMENTS, char *name, int array_size,
+                                 CCTK_REAL const th[], CCTK_REAL const ph[],
+                                 CCTK_REAL const xs[], CCTK_REAL const ys[], CCTK_REAL const zs[],
+                                 CCTK_REAL const data[])
+{
+  DECLARE_CCTK_ARGUMENTS;
+
+  if (FILE *fp = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name))
+  {
+    Multipole_OutputArray(CCTK_PASS_CTOC, fp, array_size, th, ph, xs, ys, zs, data);
+    fclose(fp);
+  }
+}
+
+////////////////////////////////////////////////////////////////////////
+// Misc
+////////////////////////////////////////////////////////////////////////
+
+void Multipole_Output1D(CCTK_ARGUMENTS, char const *name, int array_size,
+                        CCTK_REAL const th[], CCTK_REAL const ph[], mp_coord coord,
+                        CCTK_REAL const data[])
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  if (FILE *f = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name))
+  {
+    fprintf(f, "\"Time = %.19g\n", cctk_time);
+
+    if (coord == mp_theta)
+    {
+      for (int i = 0; i <= ntheta; i++)
+      {
+        int idx = Multipole_Index(i, 0, ntheta);
+        fprintf(f, "%f %.19g\n", th[idx], data[idx]);
+      }
+    }
+    else if (coord == mp_phi)
+    {
+      for (int i = 0; i <= nphi; i++)
+      {
+        int idx = Multipole_Index(ntheta / 4, i, ntheta);
+        fprintf(f, "%f %.19g\n", ph[idx], data[idx]);
+      }
+    }
+    fprintf(f, "\n\n");
+    fclose(f);
+  }
+}
+
+////////////////////////////////////////////////////////////////////////
+// Complex IO
+////////////////////////////////////////////////////////////////////////
+
+void Multipole_OutputComplex(CCTK_ARGUMENTS, FILE *fp, CCTK_REAL redata, CCTK_REAL imdata)
+{
+  DECLARE_CCTK_PARAMETERS;
+  DECLARE_CCTK_ARGUMENTS;
+  fprintf(fp, "%f %.19g %.19g\n", cctk_time, redata, imdata);
+}
+
+void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, char const *name, CCTK_REAL redata, CCTK_REAL imdata)
+{
+  DECLARE_CCTK_ARGUMENTS;
+
+  if (FILE *fp = Multipole_OpenOutputFile(CCTK_PASS_CTOC, name))
+  {
+    Multipole_OutputComplex(CCTK_PASS_CTOC, fp, redata, imdata);
+    fclose(fp);
+  }
+}
+
+
+////////////////////////////////////////////////////////////////////////
+// Coordinates
+////////////////////////////////////////////////////////////////////////
+
+void Multipole_CoordSetup(int ntheta, int nphi, 
+                          CCTK_REAL xhat[], CCTK_REAL yhat[], 
+                          CCTK_REAL zhat[], CCTK_REAL th[], 
+                          CCTK_REAL ph[])
+{
+  const CCTK_REAL PI = acos(-1.0);
+
+  for (int it = 0; it <= ntheta; it++)
+  {
+    for (int ip = 0; ip <= nphi; ip++)
+    {
+      const int i = Multipole_Index(it, ip, ntheta);
+
+      th[i] = it * PI / (ntheta);
+      ph[i] = ip * 2 * PI / nphi;
+      xhat[i] = cos(ph[i])*sin(th[i]);
+      yhat[i] = sin(ph[i])*sin(th[i]);
+      zhat[i] = cos(th[i]);
+    }
+  }
+}
+
+void Multipole_ScaleCartesian(int ntheta, int nphi, CCTK_REAL r,
+                              CCTK_REAL const xhat[], CCTK_REAL const yhat[], CCTK_REAL const zhat[],
+                              CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[])
+{
+  for (int it = 0; it <= ntheta; it++)
+  {
+    for (int ip = 0; ip <= nphi; ip++)
+    {
+      const int i = Multipole_Index(it, ip, ntheta);
+
+      x[i] = r * xhat[i];
+      y[i] = r * yhat[i];
+      z[i] = r * zhat[i];
+    }
+  }
+}
+
+////////////////////////////////////////////////////////////////////////
+// Integration
+////////////////////////////////////////////////////////////////////////
+
+void 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)
+{
+  DECLARE_CCTK_PARAMETERS
+
+  int il = Multipole_Index(0,0,ntheta);
+  int iu = Multipole_Index(1,0,ntheta);
+  CCTK_REAL dth = th[iu] - th[il];
+  iu = Multipole_Index(0,1,ntheta);
+  CCTK_REAL dph = ph[iu] - ph[il];
+
+  if (CCTK_Equals(integration_method, "midpoint"))
+  {
+    CCTK_REAL tempr = 0.0;
+    CCTK_REAL tempi = 0.0;
+      
+    for (int i=0; i < array_size; i++) {
+      // the below calculations take the integral of conj(array1)*array2
+      tempr += ( array1r[i]*array2r[i] + array1i[i]*array2i[i] )
+        *sin(th[i])*dth*dph;
+      tempi += ( array1r[i]*array2i[i] - array1i[i]*array2r[i] )
+        *sin(th[i])*dth*dph;
+      
+      *outre = tempr;
+      *outim = tempi;
+    } 
+  }
+  else if (CCTK_Equals(integration_method, "simpson"))
+  {
+    static CCTK_REAL *fr = 0;
+    static CCTK_REAL *fi = 0;
+    static bool allocated_memory = false;
+
+    // Construct an array for the real integrand
+    if (!allocated_memory)
+    {
+      printf("Multipole: Allocating memory for integrand\n");
+      fr = new CCTK_REAL[array_size];
+      fi = new CCTK_REAL[array_size];
+      allocated_memory = true;
+    }
+    
+    for (int i = 0; i < array_size; i++)
+    {
+      fr[i] = (array1r[i] * array2r[i] + 
+               array1i[i] * array2i[i] ) * sin(th[i]);
+      fi[i] = (array1r[i] * array2i[i] - 
+               array1i[i] * array2r[i] ) * sin(th[i]);
+    }
+    
+    *outre = Simpson2DIntegral(fr, ntheta, nphi, dth, dph);
+    *outim = Simpson2DIntegral(fi, ntheta, nphi, dth, dph);
+  }
+} 
diff --git a/src/utils.hh b/src/utils.hh
new file mode 100644
index 0000000..2f7d2be
--- /dev/null
+++ b/src/utils.hh
@@ -0,0 +1,40 @@
+
+#ifndef __utils_h
+#define __utils_h
+
+#include "cctk.h"
+
+enum mp_coord {mp_theta, mp_phi};
+
+void Multipole_OutputArrayToFile(CCTK_ARGUMENTS, char const *name, int array_size,
+                                 CCTK_REAL const th[], CCTK_REAL const ph[],
+                                 CCTK_REAL const x[], CCTK_REAL const y[], CCTK_REAL const z[],
+                                 CCTK_REAL const data[]);
+
+void Multipole_Output1D(CCTK_ARGUMENTS, char const *name, int array_size,
+                        CCTK_REAL const th[], CCTK_REAL const ph[], mp_coord coord,
+                        CCTK_REAL const data[]);
+
+void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, char const *name, CCTK_REAL redata, CCTK_REAL imdata);
+
+void Multipole_CoordSetup(int ntheta, int nphi, 
+                          CCTK_REAL xhat[], CCTK_REAL yhat[], 
+                          CCTK_REAL zhat[], CCTK_REAL th[], 
+                          CCTK_REAL ph[]);
+
+void Multipole_ScaleCartesian(int ntheta, int nphi, CCTK_REAL r,
+                              CCTK_REAL const xhat[], CCTK_REAL const yhat[], CCTK_REAL const zhat[],
+                              CCTK_REAL x[], CCTK_REAL y[], CCTK_REAL z[]);
+
+static inline int Multipole_Index(int it, int ip, int ntheta)
+{
+  return it + (ntheta+1)*ip;
+}
+
+void Multipole_Integrate(int array_size, int ntheta,
+    CCTK_REAL const array1r[], CCTK_REAL const array1i[],
+    CCTK_REAL const array2r[], CCTK_REAL const array2i[],
+    CCTK_REAL const th[], CCTK_REAL const pph[], 
+    CCTK_REAL out_arrayr[], CCTK_REAL out_arrayi[]);
+
+#endif