Correct theta in MSIS00IceCubeCentered

MCEq/geometry/density_profiles.py

@@ -1049,7 +1049,7 @@ def latitude(self, det_zenith_deg):
           float: latitude of the impact point in degrees
         """
         r = self.geom.r_E
-        d = 1948  # m
+        d = 1948 * 1e2 # cm (same as r_E)
 
         theta_rad = det_zenith_deg / 180.0 * np.pi
 
@@ -1058,31 +1058,25 @@ def latitude(self, det_zenith_deg):
         ) * np.cos(theta_rad)
 
         return (
-            -90.0
-            + np.arctan2(x * np.sin(theta_rad), r - d + x * np.cos(theta_rad))
+            np.arctan2(x * np.sin(theta_rad), r - d + x * np.cos(theta_rad))
             / np.pi
             * 180.0
         )
 
     def set_theta(self, theta_deg):
 
-        self._msis.set_location_coord(longitude=0.0, latitude=self.latitude(theta_deg))
+        alpha_deg = self.latitude(theta_deg)
+        theta_deg = theta_deg-alpha_deg
+
+        self._msis.set_location_coord(longitude=0.0, latitude=alpha_deg-90.0)
         info(
             1,
             "latitude = {0:5.2f} for zenith angle = {1:5.2f}".format(
-                self.latitude(theta_deg), theta_deg
+                alpha_deg-90.0, theta_deg
             ),
         )
-        downgoing_theta_deg = theta_deg
-        if theta_deg > 90.0:
-            downgoing_theta_deg = 180.0 - theta_deg
-            info(
-                1,
-                "theta = {0:5.2f} below horizon. using theta = {1:5.2f}".format(
-                    theta_deg, downgoing_theta_deg
-                ),
-            )
-        MSIS00Atmosphere.set_theta(self, downgoing_theta_deg)
+
+        MSIS00Atmosphere.set_theta(self, theta_deg)
 
         self.theta_deg = theta_deg