From 829cb257effb57c1c36ae6956896f321196794fc Mon Sep 17 00:00:00 2001
From: Jooheon Yoo <jy884@cornell.edu>
Date: Wed, 28 Aug 2024 13:25:39 -0400
Subject: [PATCH] add inverse map

---
 .../AnalyticData/GrMhd/SphericalTorus.cpp     | 26 ++++++++++++++++++-
 1 file changed, 25 insertions(+), 1 deletion(-)

diff --git a/src/PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.cpp b/src/PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.cpp
index 49ddfee7dbfb..fd52568645f4 100644
--- a/src/PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.cpp
+++ b/src/PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.cpp
@@ -7,9 +7,11 @@
 
 #include "Options/ParseError.hpp"
 #include "Utilities/ErrorHandling/Assert.hpp"
+#include "Utilities/Functional.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/MakeWithValue.hpp"
+#include "Utilities/Math.hpp"
 
 namespace grmhd::AnalyticData {
 
@@ -83,6 +85,8 @@ tnsr::I<T, 3> SphericalTorus::operator()(
 
 tnsr::I<double, 3> SphericalTorus::inverse(
     const tnsr::I<double, 3>& target_coords) const {
+  using std::abs;
+
   tnsr::I<double, 3> result;
   const double r = std::hypot(get<0>(target_coords), get<1>(target_coords),
                               get<2>(target_coords));
@@ -91,7 +95,27 @@ tnsr::I<double, 3> SphericalTorus::inverse(
                  get<2>(target_coords));
   const double phi = std::atan2(get<1>(target_coords), get<0>(target_coords));
   get<0>(result) = radius_inverse(r);
-  get<1>(result) = (M_PI_2 - theta) / pi_over_2_minus_theta_min_;
+
+  // if we are using compression (compression_level_ !=0), there is an extra
+  // step note this is just compressed eta; need to invert the cubic map back to
+  // eta
+  if (compression_level_ == 0.0) {
+    get<1>(result) = (M_PI_2 - theta) / pi_over_2_minus_theta_min_;
+  } else {
+    // using p.228 of Numerical Recipe
+    get<1>(result) = (M_PI_2 - theta) / pi_over_2_minus_theta_min_;
+    const double Q = (compression_level_ - 1.0) / (3.0 * compression_level_);
+    const double R = -0.5 * get<1>(result) / (compression_level_);
+    const double A =
+        -sgn(R) * pow(abs(R) + sqrt(square(R) - pow<3>(Q)), 1.0 / 3.0);
+    double B;
+    if (A != 0.0) {
+      B = Q / A;
+    } else {
+      B = 0.0;
+    }
+    get<1>(result) = A + B;
+  }
   get<2>(result) = phi / (M_PI * fraction_of_torus_);
   return result;
 }