From d47c3dc7324ddc372d3eeeab1366a38064fd70db Mon Sep 17 00:00:00 2001 From: knarf Date: Fri, 26 Mar 2010 10:25:54 +0000 Subject: use trunc structure git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Multipole/trunk@53 4f5cb9a8-4dd8-4c2d-9bbd-173fa4467843 --- src/utils.cc | 251 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 251 insertions(+) create mode 100644 src/utils.cc (limited to 'src/utils.cc') 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 +#include +#include + +#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); + } +} -- cgit v1.2.3