done upto inv_jacob

src/Domain/CoordinateMaps/SphericalTorus.cpp

@@ -87,8 +87,8 @@ SphericalTorus::jacobian(const std::array<T, 3>& source_coords) const {
   // In order to reduce number of memory allocations we use some slots of
   // jacobian for storing temp variables as below.
 
-  auto& sin_theta = get<2, 0>(jacobian);
-  auto& cos_theta = get<2, 1>(jacobian);
+  auto& sin_theta = get<2, 1>(jacobian);
+  auto& cos_theta = get<2, 0>(jacobian);
   get<2, 2>(jacobian) =
       M_PI_2 - (pi_over_2_minus_theta_min_ * source_coords[1]);
   sin_theta = sin(get<2, 2>(jacobian));
@@ -128,15 +128,15 @@ SphericalTorus::inv_jacobian(const std::array<T, 3>& source_coords) const {
   // In order to reduce number of memory allocations we use some slots of
   // jacobian for storing temp variables as below.
 
-  auto& sin_theta = get<1, 1>(inv_jacobian);
-  auto& cos_theta = get<2, 1>(inv_jacobian);
+  auto& sin_theta = get<1, 2>(inv_jacobian);
+  auto& cos_theta = get<0, 2>(inv_jacobian);
   get<2, 2>(inv_jacobian) =
       M_PI_2 - (pi_over_2_minus_theta_min_ * source_coords[1]);
   cos_theta = cos(get<2, 2>(inv_jacobian));
   sin_theta = sin(get<2, 2>(inv_jacobian));
 
-  auto& sin_phi = get<0, 2>(inv_jacobian);
-  auto& cos_phi = get<1, 2>(inv_jacobian);
+  auto& sin_phi = get<1, 1>(inv_jacobian);
+  auto& cos_phi = get<2, 1>(inv_jacobian);
   get<2, 2>(inv_jacobian) = M_PI * fraction_of_torus_ * source_coords[2];
   cos_phi = cos(get<2, 2>(inv_jacobian));
   sin_phi = sin(get<2, 2>(inv_jacobian));
@@ -146,18 +146,18 @@ SphericalTorus::inv_jacobian(const std::array<T, 3>& source_coords) const {
 
   // Note : execution order matters here since we are overwriting each temp
   // variables with jacobian values corresponding to the slot.
-  get<0, 0>(inv_jacobian) = 2.0 / (r_max_ - r_min_) * cos_theta * cos_phi;
-  get<0, 1>(inv_jacobian) = 2.0 / (r_max_ - r_min_) * sin_theta * cos_phi;
+  get<0, 0>(inv_jacobian) = 2.0 / (r_max_ - r_min_) * sin_theta * cos_phi;
+  get<0, 1>(inv_jacobian) = 2.0 / (r_max_ - r_min_) * sin_theta * sin_phi;
   get<1, 0>(inv_jacobian) =
-      -(1.0 / pi_over_2_minus_theta_min_) * cos_theta * sin_phi / r;
+      -(1.0 / pi_over_2_minus_theta_min_) * cos_theta * cos_phi / r;
   get<2, 0>(inv_jacobian) =
-      -(1.0 / (M_PI * fraction_of_torus_)) * sin_theta / (r * cos_phi);
+      -(1.0 / (M_PI * fraction_of_torus_)) * sin_phi / (r * sin_theta);
   get<1, 1>(inv_jacobian) =
       -(1.0 / pi_over_2_minus_theta_min_) * sin_theta * sin_phi / r;
   get<2, 1>(inv_jacobian) =
-      (1.0 / (M_PI * fraction_of_torus_)) * cos_theta / (r * cos_phi);
-  get<0, 2>(inv_jacobian) = 2.0 / (r_max_ - r_min_) * sin_phi;
-  get<1, 2>(inv_jacobian) = (1.0 / pi_over_2_minus_theta_min_) * cos_phi / r;
+      (1.0 / (M_PI * fraction_of_torus_)) * cos_phi / (r * sin_theta);
+  get<0, 2>(inv_jacobian) = 2.0 / (r_max_ - r_min_) * cos_theta;
+  get<1, 2>(inv_jacobian) = (1.0 / pi_over_2_minus_theta_min_) * sin_theta / r;
   get<2, 2>(inv_jacobian) = 0.0;
 
   return inv_jacobian;