add inverse map

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;
 }