From bc7d557b77d270c20d25922208632db7884a20c6 Mon Sep 17 00:00:00 2001
From: Liwei Ji <jiliwei.phys@gmail.com>
Date: Wed, 26 Jun 2024 11:28:30 +0000
Subject: [PATCH] Multipole: beautify class Surface and Sphere

---
 Multipole/src/surface.cxx | 44 ---------------------------
 Multipole/src/surface.hxx | 62 +++++++++++++++++++++++++++++++--------
 2 files changed, 50 insertions(+), 56 deletions(-)

diff --git a/Multipole/src/surface.cxx b/Multipole/src/surface.cxx
index 07e429b5..27b44e12 100644
--- a/Multipole/src/surface.cxx
+++ b/Multipole/src/surface.cxx
@@ -5,48 +5,4 @@
 
 namespace Multipole {
 
-Surface::Surface(int nTheta, int nPhi, bool isMidPoint)
-    : nTheta_(nTheta), nPhi_(nPhi) {
-  constexpr CCTK_REAL PI = acos(-1.0);
-  const int arraySize = (nTheta + 1) * (nPhi + 1);
-  theta_.resize(arraySize);
-  phi_.resize(arraySize);
-  x_.resize(arraySize);
-  y_.resize(arraySize);
-  z_.resize(arraySize);
-
-  // Add an offset for midpoint integration.
-  const CCTK_REAL is_midpoint = static_cast<CCTK_REAL>(isMidPoint);
-  const CCTK_REAL dTheta = PI / (nTheta + is_midpoint);
-  const CCTK_REAL dPhi = 2 * PI / (nPhi + is_midpoint);
-
-  for (int it = 0; it <= nTheta; ++it) {
-    for (int ip = 0; ip <= nPhi; ++ip) {
-      const int i = index2D(it, ip);
-      theta_[i] = it * dTheta + 0.5 * dTheta * is_midpoint;
-      phi_[i] = ip * dPhi + 0.5 * dPhi * is_midpoint;
-    }
-  }
-}
-
-Sphere::Sphere(int nTheta, int nPhi, bool isMidPoint)
-    : Surface(nTheta, nPhi, isMidPoint) {
-  xhat_.resize(theta_.size());
-  yhat_.resize(theta_.size());
-  zhat_.resize(theta_.size());
-  for (size_t i = 0; i < theta_.size(); ++i) {
-    xhat_[i] = cos(phi_[i]) * sin(theta_[i]);
-    yhat_[i] = sin(phi_[i]) * sin(theta_[i]);
-    zhat_[i] = cos(theta_[i]);
-  }
-}
-
-void Sphere::setRadius(CCTK_REAL radius) {
-  for (size_t i = 0; i < theta_.size(); ++i) {
-    x_[i] = radius * xhat_[i];
-    y_[i] = radius * yhat_[i];
-    z_[i] = radius * zhat_[i];
-  }
-}
-
 } // namespace Multipole
diff --git a/Multipole/src/surface.hxx b/Multipole/src/surface.hxx
index ab85e704..8f5a6237 100644
--- a/Multipole/src/surface.hxx
+++ b/Multipole/src/surface.hxx
@@ -13,29 +13,67 @@ namespace Multipole {
 // 2D Surface Embedding in 3D Space
 class Surface {
 public:
-  Surface(int nTheta, int nPhi, bool isMidPoint);
-  virtual ~Surface() {} // Ready for potential inheritance
+  Surface(int nTheta, int nPhi, bool isMidPoint)
+      : nTheta_(nTheta), nPhi_(nPhi) {
+    constexpr CCTK_REAL PI = std::acos(-1.0);
+    const int arraySize = (nTheta + 1) * (nPhi + 1);
+    theta_.resize(arraySize);
+    phi_.resize(arraySize);
+    x_.resize(arraySize);
+    y_.resize(arraySize);
+    z_.resize(arraySize);
 
-  std::vector<CCTK_REAL> getTheta() const { return theta_; }
-  std::vector<CCTK_REAL> getPhi() const { return phi_; }
+    // Add an offset for midpoint integration.
+    const CCTK_REAL is_midpoint = static_cast<CCTK_REAL>(isMidPoint);
+    const CCTK_REAL dTheta = PI / (nTheta + is_midpoint);
+    const CCTK_REAL dPhi = 2 * PI / (nPhi + is_midpoint);
 
-private:
-  friend class Sphere; // allow Sphere to access private members
+    for (int it = 0; it <= nTheta; ++it) {
+      for (int ip = 0; ip <= nPhi; ++ip) {
+        const int i = index2D(it, ip);
+        theta_[i] = it * dTheta + 0.5 * dTheta * is_midpoint;
+        phi_[i] = ip * dPhi + 0.5 * dPhi * is_midpoint;
+      }
+    }
+  }
 
-  const int nTheta_, nPhi_;
-  std::vector<CCTK_REAL> theta_, phi_;
-  std::vector<CCTK_REAL> x_, y_, z_; // embedding information
+  virtual ~Surface() = default; // Use default for trivial destructors
+
+  const std::vector<CCTK_REAL> &getTheta() const { return theta_; }
+  const std::vector<CCTK_REAL> &getPhi() const { return phi_; }
 
+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);
-  ~Sphere() {} // Ready for potential inheritance
+  Sphere(int nTheta, int nPhi, bool isMidPoint)
+      : Surface(nTheta, nPhi, isMidPoint) {
+    xhat_.resize(theta_.size());
+    yhat_.resize(theta_.size());
+    zhat_.resize(theta_.size());
+    for (size_t i = 0; i < theta_.size(); ++i) {
+      xhat_[i] = std::cos(phi_[i]) * std::sin(theta_[i]);
+      yhat_[i] = std::sin(phi_[i]) * std::sin(theta_[i]);
+      zhat_[i] = std::cos(theta_[i]);
+    }
+  }
+
+  ~Sphere() override = default;
 
-  void setRadius(CCTK_REAL radius);
+  void setRadius(CCTK_REAL radius) {
+    for (size_t i = 0; i < theta_.size(); ++i) {
+      x_[i] = radius * xhat_[i];
+      y_[i] = radius * yhat_[i];
+      z_[i] = radius * zhat_[i];
+    }
+  }
 
 private:
   std::vector<CCTK_REAL> xhat_, yhat_, zhat_; // unit sphere