Skip to content

Correct theta in MSIS00IceCubeCentered #51

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
pochoarus wants to merge 6 commits into
base: main
Choose a base branch
from
Open

Correct theta in MSIS00IceCubeCentered #51

Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
8d5b58c
Correct not only latitude but also theta_deg in MSIS00IceCubeCentered
pochoarus Jan 27, 2023
26d0cc6
no need to computed the downgoing since theta_deg<90 with new correction
pochoarus Jan 27, 2023
e7a9aa1
use theta_deg instead of downgoing_theta_deg
pochoarus Jan 27, 2023
3d4c41f
change print
pochoarus Jan 28, 2023
c76b734
depth must be in cm as r_E
pochoarus Jan 28, 2023
a9b7582
change name of the function
pochoarus Jan 28, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 12 additions & 17 deletions MCEq/geometry/density_profiles.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1034,8 +1034,9 @@ def __init__(self, location, season):
# Allow for upgoing zenith angles
self.max_theta = 180.0

def latitude(self, det_zenith_deg):
"""Returns the geographic latitude of the shower impact point.
def depth_correction(self, det_zenith_deg):
"""Returns the geographic latitude (plus 90 deg.) of the
shower impact point.

Assumes a spherical earth. The detector is 1948m under the
surface.
Expand All @@ -1049,7 +1050,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

Expand All @@ -1058,31 +1059,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.depth_correction(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

Expand Down