Multipole: add spherical harmonics to class Sphere

Multipole/src/multipole.cxx

@@ -95,8 +95,8 @@ static void get_spin_weights(vector<VariableParse> vars,
 }
 
 static void
-setup_harmonics(const vector<int> spin_weights, int lmax, CCTK_REAL th[],
-                CCTK_REAL ph[], int array_size,
+setup_harmonics(const vector<int> spin_weights, int lmax, vector<CCTK_REAL> th,
+                vector<CCTK_REAL> ph, int array_size,
                 vector<vector<vector<vector<CCTK_REAL> > > > &realY,
                 vector<vector<vector<vector<CCTK_REAL> > > > &imagY) {
   for (size_t si = 0; si < spin_weights.size(); si++) {
@@ -197,8 +197,7 @@ extern "C" void Multipole_Setup(CCTK_ARGUMENTS) {
   parse_variables_string(string(variables), vars);
   get_spin_weights(vars, spin_weights);
   CoordSetup(xhat.data(), yhat.data(), zhat.data(), th.data(), ph.data());
-  setup_harmonics(spin_weights, l_max, th.data(), ph.data(), array_size, realY,
-                  imagY);
+  setup_harmonics(spin_weights, l_max, th, ph, array_size, realY, imagY);
   CCTK_VINFO("initialized arrays");
 }
 

Multipole/src/sphericalharmonic.cxx

@@ -48,8 +48,9 @@ void SphericalHarmonic(int s, int l, int m, CCTK_REAL th, CCTK_REAL ph,
   *imY = all_coeff * sum * sin(m * ph);
 }
 
-void HarmonicSetup(int s, int l, int m, int array_size, CCTK_REAL const th[],
-                   CCTK_REAL const ph[], std::vector<CCTK_REAL> &reY,
+void HarmonicSetup(int s, int l, int m, int array_size,
+                   const std::vector<CCTK_REAL> th,
+                   const std::vector<CCTK_REAL> ph, std::vector<CCTK_REAL> &reY,
                    std::vector<CCTK_REAL> &imY) {
   for (int i = 0; i < array_size; i++) {
     SphericalHarmonic(s, l, m, th[i], ph[i], &reY[i], &imY[i]);

Multipole/src/sphericalharmonic.hxx

@@ -12,8 +12,9 @@ namespace Multipole {
 void SphericalHarmonic(int s, int l, int m, CCTK_REAL th, CCTK_REAL ph,
                        CCTK_REAL *reY, CCTK_REAL *imY);
 
-void HarmonicSetup(int s, int l, int m, int array_size, CCTK_REAL const th[],
-                   CCTK_REAL const ph[], std::vector<CCTK_REAL> &reY,
+void HarmonicSetup(int s, int l, int m, int array_size,
+                   const std::vector<CCTK_REAL> th,
+                   const std::vector<CCTK_REAL> ph, std::vector<CCTK_REAL> &reY,
                    std::vector<CCTK_REAL> &imY);
 
 } // namespace Multipole

Multipole/src/surface.hxx

@@ -4,6 +4,8 @@
 #ifndef MULTIPOLE_SURFACE_HXX
 #define MULTIPOLE_SURFACE_HXX
 
+#include "sphericalharmonic.hxx"
+
 #include <cctk.h>
 
 #include <vector>
@@ -46,15 +48,18 @@ protected:
   inline int index2D(int it, int ip) { return it + (nTheta_ + 1) * ip; }
 
   const int nTheta_, nPhi_;
+
   std::vector<CCTK_REAL> theta_, phi_;
   std::vector<CCTK_REAL> x_, y_, z_; // embedding map
 };
 
 // 2D Sphere
 class Sphere : public Surface {
 public:
-  Sphere(int nTheta, int nPhi, bool isMidPoint)
-      : Surface(nTheta, nPhi, isMidPoint) {
+  Sphere(int nTheta, int nPhi, bool isMidPoint, std::vector<int> spinWeights,
+         int lmax)
+      : Surface(nTheta, nPhi, isMidPoint), spin_weights_(spinWeights),
+        lmax_(lmax) {
     xhat_.resize(theta_.size());
     yhat_.resize(theta_.size());
     zhat_.resize(theta_.size());
@@ -63,6 +68,29 @@ public:
       yhat_[i] = std::sin(phi_[i]) * std::sin(theta_[i]);
       zhat_[i] = std::cos(theta_[i]);
     }
+
+    realY_.resize(spinWeights.size());
+    imagY_.resize(spinWeights.size());
+
+    for (size_t si = 0; si < spinWeights.size(); ++si) {
+      int sw = spinWeights[si];
+      realY_[si].resize(lmax + 1);
+      imagY_[si].resize(lmax + 1);
+
+      for (int l = 0; l <= lmax; ++l) {
+        realY_[si][l].resize(2 * l + 1);
+        imagY_[si][l].resize(2 * l + 1);
+
+        for (int m = -l; m <= l; ++m) {
+          std::vector<CCTK_REAL> realY_s_l_m(theta_.size());
+          std::vector<CCTK_REAL> imagY_s_l_m(theta_.size());
+          HarmonicSetup(sw, l, m, theta_.size(), theta_, phi_, realY_s_l_m,
+                        imagY_s_l_m);
+          realY_[si][l][m + l] = realY_s_l_m;
+          imagY_[si][l][m + l] = imagY_s_l_m;
+        }
+      }
+    }
   }
 
   ~Sphere() override = default;
@@ -77,6 +105,11 @@ public:
 
 private:
   std::vector<CCTK_REAL> xhat_, yhat_, zhat_; // unit sphere
+  std::vector<int> spin_weights_;             // spin weights
+  // spin weighted spherical harmonics
+  int lmax_;
+  std::vector<std::vector<std::vector<std::vector<CCTK_REAL> > > > realY_;
+  std::vector<std::vector<std::vector<std::vector<CCTK_REAL> > > > imagY_;
 };
 
 } // namespace Multipole

Multipole/src/tests.cxx

@@ -134,13 +134,14 @@ extern "C" void Multipole_TestOrthonormality(CCTK_ARGUMENTS) {
   /* Campute Cartesian coordinates of points on the sphere */
   int array_size = (ntheta + 1) * (nphi + 1);
 
-  CCTK_REAL *th = new CCTK_REAL[array_size];
-  CCTK_REAL *ph = new CCTK_REAL[array_size];
+  vector<CCTK_REAL> th, ph;
+  th.resize(array_size);
+  ph.resize(array_size);
   CCTK_REAL *xhat = new CCTK_REAL[array_size];
   CCTK_REAL *yhat = new CCTK_REAL[array_size];
   CCTK_REAL *zhat = new CCTK_REAL[array_size];
 
-  CoordSetup(xhat, yhat, zhat, th, ph);
+  CoordSetup(xhat, yhat, zhat, th.data(), ph.data());
 
   /* Populate spherical-harmonic array */
   vector<vector<vector<vector<CCTK_REAL> > > > reY, imY;
@@ -184,7 +185,8 @@ extern "C" void Multipole_TestOrthonormality(CCTK_ARGUMENTS) {
             CCTK_REAL real_lm = 0.0, imag_lm = 0.0;
             Integrate(array_size, ntheta, reY[sw][li][mi + li],
                       imY[sw][li][mi + li], reY[sw][l][m + l],
-                      imY[sw][l][m + l], th, ph, &real_lm, &imag_lm);
+                      imY[sw][l][m + l], th.data(), ph.data(), &real_lm,
+                      &imag_lm);
 
             assert(idx < 1 * N * (N + 1) / 2);
             test_orthonormality[idx++] =
@@ -199,8 +201,6 @@ extern "C" void Multipole_TestOrthonormality(CCTK_ARGUMENTS) {
   delete[] zhat;
   delete[] yhat;
   delete[] xhat;
-  delete[] ph;
-  delete[] th;
 
   return;
 }