1 files changed, 251 insertions, 0 deletions
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);
+  }
+}