# -*- coding: utf-8 -*-
# SPDX-License-Identifier: GPL-3.0-only
# Copyright (C) 2026 Mauricio Bustamante
r"""earth.py
Contains helper functions related to the Earth: its internal matter
density and the geometry of neutrino trajectories through it.
Routine listings
----------------
* density_matter_func_prem - Returns the density inside the Earth
using the Preliminary Reference Earth Model (PREM), with an
optional override of the outermost shell's density
* prem_layer_edges_along_chord - Returns the positions at which a
chord through the Earth crosses the PREM layer boundaries
* distance_traveled_inside_earth - Returns the chord length for a
given neutrino direction, with either endpoint optionally
underground
* earth_radial_distance_from_depth - Converts position along a
chord to radial distance from the center of the Earth
* dms_to_decimal - Converts (degree, minute, second) coordinates to
decimal degrees
* chord_length_inside_earth - Returns the chord length between two
locations on the surface of the Earth
* costhz_between_points_on_surface - Returns the zenith angle of
the chord between two locations on the surface of the Earth
* coordinates_of_named_location - Returns the coordinates of a
predefined location (e.g., a neutrino detector site)
* electron_fraction_func_prem - Returns Y_e at one or more radii,
resolved per PREM layer (iron core, rock mantle, crust, ocean)
* neutron_to_proton_ratio_from_electron_fraction - Returns the
neutron-to-proton ratio implied by an electron fraction,
r = (1 - Y_e)/Y_e
"""
__author__ = "Mauricio Bustamante"
__email__ = "mbustamante@gmail.com"
import numpy as np
from typing import Optional, Union
import magnus.globaldefs as gd
# Predefined locations in ISO 6709:
# North latitudes are positive, South latitudes are negative
# East longitudes are positive, West longitudes are negative
[docs]
loc_coords_dms = {
'baikal': {'lat': (51, 45, 54), 'lon': (104, 24, 54)},
'cern': {'lat': (46, 14, 1.80), 'lon': (6, 3, 11.40)},
'desy': {'lat': (53, 34, 19.79), 'lon': (9, 52, 27.59)},
'ess': {'lat': (55, 44, 6), 'lon': (13, 15, 5.04)},
'fermilab': {'lat': (41, 49, 55), 'lon': (-88, 15, 26)},
'gran_sasso': {'lat': (42, 25, 15.8), 'lon': (13, 30, 58.43)},
'homestake': {'lat': (44, 21, 5.76), 'lon': (-103, 45, 4.68)},
'kamioka': {'lat': (36, 25, 50.05), 'lon': (137, 18, 41.15)}, # Mozumi mine
'km3net_arca': {'lat': (36, 16, 0), 'lon': (16, 6, 0)},
'km3net_orca': {'lat': (42, 48, 0), 'lon': (6, 2, 0)},
'north_pole': {'lat': (90, 0, 0), 'lon': (0, 0, 0)},
'pyhaasalmi': {'lat': (63, 39, 31), 'lon': (26, 2, 28)},
'snolab': {'lat': (46, 28, 18), 'lon': (-81, 11, 12)},
'south_pole': {'lat': (-90, 0, 0), 'lon': (0, 0, 0)},
'tokai': {'lat': (36, 27, 59), 'lon': (140, 36, 24)},
}
# PREM layers: inner radial boundary of each shell [km] (the last shell ends
# at the surface, gd.EARTH_RADIUS), and the coefficients (c0, c1, c2, c3) of
# the density polynomial rho(x) = c0 + c1*x + c2*x^2 + c3*x^3, with
# x = r/EARTH_RADIUS, inside each shell (Dziewonski & Anderson 1981).
[docs]
PREM_BOUNDARIES = np.array([1221.5, 3480.0, 5701.0, 5771.0, 5971.0, 6151.0,
6346.6, 6356.0, 6368.0])
_PREM_COEFFS = np.array([
[13.0885, 0.0, -8.8381, 0.0],
[12.5815, -1.2638, -3.6426, -5.5281],
[ 7.9565, -6.4761, 5.5283, -3.0807],
[ 5.3197, -1.4836, 0.0, 0.0],
[11.2494, -8.0298, 0.0, 0.0],
[ 7.1089, -3.8045, 0.0, 0.0],
[ 2.6910, 0.6924, 0.0, 0.0],
[ 2.900, 0.0, 0.0, 0.0],
[ 2.600, 0.0, 0.0, 0.0],
[ 1.020, 0.0, 0.0, 0.0],
])
[docs]
def density_matter_func_prem(r: Union[float, np.ndarray],
tol: Optional[float]=1.e-8,
density_matter_ocean: Optional[float]=None) -> Union[float, np.ndarray]:
r"""Returns the matter density inside the Earth according to the
Preliminary Reference Earth Model (PREM) [1]_.
Returns the matter density inside the Earth according to the PREM,
for a given radial distance measured from the center of the Earth.
Accepts a single radial distance or an array of radial distances;
array input is evaluated in a single vectorized pass.
.. versionadded:: 1.0.0
.. versionchanged:: 1.1.1
Added ``density_matter_ocean``, which replaces the density of
PREM's outermost shell. Left as None, the profile is unchanged.
Parameters
----------
r : float or np.ndarray
Radial distance(s) measured from the center of the Earth [km].
tol : float, optional
Relative tolerance by which a radial distance may exceed
``globaldefs.EARTH_RADIUS`` before a ValueError is raised;
radii within the tolerance are clamped onto the surface.
Default: 1e-8.
density_matter_ocean : float, optional
Density of the outermost PREM shell, :math:`r > 6368` km
[:math:`\text{g cm}^{-3}`]. PREM puts a global-average ocean
there, at 1.020; a detector under continental rock sits under
about 2.6 instead, and one under Antarctic ice under about 0.92.
Pass ``electron_fraction_ocean`` alongside it to set the
composition of the same shell, which
:func:`electron_fraction_func_prem` handles. Default: None, i.e.
PREM's own ocean.
Returns
-------
float or np.ndarray
Matter density [:math:`\text{g cm}^{-3}`].
Raises
------
ValueError
If any radial distance exceeds globaldefs.EARTH_RADIUS by more
than the relative tolerance tol.
References
----------
.. [1] Adam M. Dziewonski & Don L. Anderson, "Preliminary Reference
Earth Model", Physics of the Earth and Planetary Interiors, 25,
297 (1981).
Examples
--------
.. jupyter-execute::
from magnus import earth
for r in (0.0, 3000.0, 5000.0, 6371.0):
print('r = %6.0f km -> %6.2f g/cm^3'
% (r, earth.density_matter_func_prem(r)))
"""
scalar_input = (np.ndim(r) == 0)
r = np.asarray(r, dtype=float)
x = r/gd.EARTH_RADIUS
if np.any(x - 1.0 > tol):
raise ValueError('Error in magnus: earth.density_matter_func_prem: value of r cannot exceed ' + \
'globaldefs.EARTH_RADIUS = ' + str(gd.EARTH_RADIUS) + ' km by more than the ' + \
'desired tolerance of tol = ' + str(tol))
# Clamp radii within tolerance of the surface onto the surface
r = np.minimum(r, gd.EARTH_RADIUS)
x = np.minimum(x, 1.0)
# Look up the PREM shell of each radius (side='left' reproduces the
# right-closed bins of the piecewise definition, e.g., r <= 1221.5) and
# evaluate the density polynomial via Horner's rule. This is ~10x
# faster than an np.select over the ten shells.
shell = np.searchsorted(PREM_BOUNDARIES, r, side='left')
c = _PREM_COEFFS[shell]
density = c[..., 0] + x*(c[..., 1] + x*(c[..., 2] + x*c[..., 3]))
# PREM's outermost shell is 3 km of global-average ocean. A detector under rock or
# ice has none, and for a trajectory close to horizontal that shell can be the whole
# path, so a caller who knows what their outermost 3 km is made of can say so. The
# substitution is by shell index rather than by radius, so it tracks the boundary at
# 6368 km wherever the lookup puts it. Left as None, nothing is substituted and the
# returned array is the one every earlier version returned.
if density_matter_ocean is not None:
density = np.where(shell == len(PREM_BOUNDARIES), float(density_matter_ocean),
density)
return float(density) if scalar_input else density
def _depths_or_zero(source_depth: Optional[float],
detector_depth: Optional[float]) -> tuple[float, float]:
r"""Returns the two depths with None read as zero, i.e. as an endpoint on the surface.
Both parameters are declared Optional, so None has to mean something; the only thing
it can mean is "no depth". Normalizing here rather than at each use keeps the test
that selects the default code path (``source_depth == 0.0 and detector_depth == 0.0``)
from sending a None down the general branch, where it would surface as a TypeError
from ``float(None)`` instead of this package's descriptive ValueError.
.. versionadded:: 1.1.1
"""
return (0.0 if source_depth is None else source_depth,
0.0 if detector_depth is None else detector_depth)
def _validated_endpoint_radii(costhz: float, source_depth: float, detector_depth: float,
source_func_name: str) -> tuple[float, float]:
r"""Returns (r_source, r_detector) [km] for two depths below the surface.
Shared by the three trajectory functions so that one depth cannot be rejected by one
of them and accepted by another. A depth equal to the Earth's radius would put an
endpoint at the center, where the zenith angle no longer names a direction, so the
interval is half open.
The cosine is checked here too, and only here, which means only on the buried branch.
A cosine outside [-1, 1] is not a direction, and the general trajectory formulas take
the square root of :math:`1 - \cos^2\theta_z` and would return NaN for one -- quietly,
since NaN propagates all the way to a probability. The surface branch keeps its own
long-standing behavior of returning a number for such input rather than raising,
because tightening it would change results that already exist.
.. versionadded:: 1.1.1
"""
R = gd.EARTH_RADIUS
if not (-1.0 <= float(costhz) <= 1.0):
raise ValueError('Error in magnus: earth.' + source_func_name + ': costhz is the ' + \
'cosine of the zenith angle, so it must lie in [-1, 1]; got ' + str(costhz) + '.')
for name, value in (('source_depth', source_depth), ('detector_depth', detector_depth)):
if not (0.0 <= float(value) < R):
raise ValueError('Error in magnus: earth.' + source_func_name + ': ' + name + \
' is measured below the surface of the Earth, so it must lie in ' + \
'[0, globaldefs.EARTH_RADIUS = ' + str(R) + ') km; got ' + str(value) + '.')
return R - float(source_depth), R - float(detector_depth)
[docs]
def distance_traveled_inside_earth(costhz: float, source_depth: Optional[float]=0.0,
detector_depth: Optional[float]=0.0) -> float:
r"""Returns the distance traveled by a neutrino inside the Earth,
traveling with a cosine of zenith angle costhz.
Returns the length of the path traveled by a neutrino from its point
of entry into the Earth, through it, until it reaches a detector. The
direction of the neutrino is parametrized by the zenith angle of the
neutrino, **measured at the detector**.
By default the source and the detector both sit on the surface, which
is the geometry every earlier version assumed: the path is the full
chord, and its length is zero for all values of costhz > 0. Burying
either end moves the corresponding endpoint to a smaller radius. The
zenith angle keeps its meaning throughout, since at zero depth the
detector is on the surface; a buried detector sees a downward-going
neutrino (costhz > 0) through its overburden, so the path length is
then positive rather than zero.
Writing :math:`r_{\rm s} = R_\oplus - {}` ``source_depth`` and
:math:`r_{\rm d} = R_\oplus - {}` ``detector_depth``, the impact
parameter of the trajectory is :math:`b = r_{\rm d} \sqrt{1 -
\cos^2\theta_z}` and its length is
.. math::
L = \sqrt{r_{\rm s}^2 - b^2} - r_{\rm d} \cos\theta_z ~.
.. versionadded:: 1.0.0
.. versionchanged:: 1.1.1
Added ``source_depth`` and ``detector_depth``. Their defaults of
zero reproduce the surface-to-surface chord bit for bit, through
the same expression as before.
Parameters
----------
costhz : float
Cosine of the zenith angle of the neutrino, measured at the
detector; must lie in [-1, 1]. The bound is enforced only when an
endpoint is buried: with both ends on the surface, a costhz outside
it returns a number rather than raising.
source_depth : float or None, optional
Depth of the entry point below the surface [km]. None is read as
0.0. Default: 0.0, i.e. the neutrino enters at the surface.
detector_depth : float or None, optional
Depth of the detector below the surface [km]. None is read as 0.0.
Default: 0.0, i.e. the detector sits on the surface.
Returns
-------
float
Path length inside the Earth [km]. Every ``osc_prob_*`` baseline is
in :math:`\text{eV}^{-1}`, so multiply by
:data:`magnus.globaldefs.UNIT_KM` before passing this on; handing the
raw value over returns a converged, unitary, wrong answer.
Raises
------
ValueError
If ``costhz`` is outside [-1, 1] and an endpoint is buried, if either
depth is outside [0, R_earth), if the trajectory never reaches the
source radius, which happens when the source is buried below the
trajectory's closest approach to the center, or if the resulting path
length is negative, which happens when the source sits deeper than the
detector on a downward-going trajectory.
Examples
--------
.. jupyter-execute::
from magnus import earth
for costhz in (-0.2, -0.5, -1.0):
print('costhz = %5.2f -> %8.1f km'
% (costhz, earth.distance_traveled_inside_earth(costhz)))
A detector 2 km underground sees a shorter upward-going path, and a
downward-going one through its overburden:
.. jupyter-execute::
for costhz in (-1.0, 0.5, 1.0):
print('costhz = %5.2f -> %10.4f km'
% (costhz, earth.distance_traveled_inside_earth(
costhz, detector_depth=2.0)))
"""
source_depth, detector_depth = _depths_or_zero(source_depth, detector_depth)
# The default geometry returns through the expression this function has always used,
# rather than through the general one below. The two agree bitwise on every zenith
# angle tested (220001 of them, randomized and gridded), but "agrees on everything
# measured" is a weaker promise than "is the same expression", and every result this
# package has published was computed with this line.
if source_depth == 0.0 and detector_depth == 0.0:
return 0.0 if costhz > 0.0 else -2.0 * gd.EARTH_RADIUS * costhz
r_s, r_d = _validated_endpoint_radii(costhz, source_depth, detector_depth,
'distance_traveled_inside_earth')
# Grouped as (r_s - r_d)(r_s + r_d) + (r_d costhz)^2 rather than the algebraically
# equal r_s^2 - r_d^2 (1 - costhz^2). The second form subtracts two numbers near
# R_earth^2 and loses the leading digits when the two radii are close, which is the
# common case; the first is exact when they are equal.
disc = (r_s - r_d)*(r_s + r_d) + (r_d*costhz)**2
if disc < 0.0:
raise ValueError('Error in magnus: earth.distance_traveled_inside_earth: the ' + \
'trajectory never reaches the source radius. Its closest approach to the ' + \
'center is ' + str(r_d*np.sqrt(1.0-costhz*costhz)) + ' km, which is above ' + \
'the source at ' + str(r_s) + ' km. Move the source outwards, or point the ' + \
'trajectory closer to the vertical.')
length = np.sqrt(disc) - r_d*costhz
if length < 0.0:
raise ValueError('Error in magnus: earth.distance_traveled_inside_earth: a ' + \
'downward-going trajectory (costhz = ' + str(costhz) + ' >= 0) travels ' + \
'outward as it is traced back from the detector, so it can only start at a ' + \
'radius above the detector. The source is at ' + str(r_s) + ' km and the ' + \
'detector at ' + str(r_d) + ' km.')
return float(length)
[docs]
def earth_radial_distance_from_depth(costhz: float, l: Union[float, np.ndarray],
tol: Optional[float]=1.e-8, source_depth: Optional[float]=0.0,
detector_depth: Optional[float]=0.0) -> Union[float, np.ndarray]:
r"""Returns the radial distance measured from the center of the
Earth to a position inside the Earth, given by costhz and l.
A neutrino with direction given by the cosine of the zenith angle,
costhz, travels from l=0 to l=distance_traveled_inside_earth,
computed below. The routine returns the radial distance to the
neutrino when its distance from its point of entry into the Earth is
l. Accepts a single distance or an array of distances; array input
is evaluated in a single vectorized pass.
The two depths move the endpoints of that trajectory, exactly as they
do in :func:`distance_traveled_inside_earth`, and carry the same
defaults: the entry point is on the surface and so is the detector.
Whatever the depths, ``l`` is measured from the entry point, so
``l = 0`` returns the source radius and ``l = L`` the detector's.
.. versionadded:: 1.0.0
.. versionchanged:: 1.1.1
Added ``source_depth`` and ``detector_depth``. Their defaults of
zero reproduce the surface-to-surface chord bit for bit, through
the same expression as before.
Parameters
----------
costhz : float
Cosine of the zenith angle of the neutrino, measured at the
detector.
l : float or np.ndarray
Distance(s) of the neutrino from its point of entry into the
Earth [km].
tol : float, optional
Absolute tolerance [km] by which ``l`` may exceed the distance
traveled inside the Earth before a ValueError is raised;
distances within the tolerance are clamped onto the exit point.
Default: 1e-8.
source_depth : float, optional
Depth of the entry point below the surface [km]. Default: 0.0.
detector_depth : float, optional
Depth of the detector below the surface [km]. Default: 0.0.
Returns
-------
float or np.ndarray
Radial distance to the neutrino [km].
Raises
------
ValueError
If any l exceeds the distance traveled inside the Earth for
this value of costhz by more than the tolerance tol, or if
either depth is rejected by
:func:`distance_traveled_inside_earth`.
Examples
--------
.. jupyter-execute::
from magnus import earth
L = earth.distance_traveled_inside_earth(-0.8, detector_depth=2.0)
for l in (0.0, 0.5*L, L):
print('l = %9.3f km -> r = %8.3f km'
% (l, earth.earth_radial_distance_from_depth(
-0.8, l, detector_depth=2.0)))
"""
source_depth, detector_depth = _depths_or_zero(source_depth, detector_depth)
scalar_input = (np.ndim(l) == 0)
l = np.asarray(l, dtype=float)
d = distance_traveled_inside_earth(costhz, source_depth, detector_depth)
if np.any(l - d > tol):
raise ValueError('Error in magnus: earth_radial_distance_from_depth: value of ' + \
'l cannot be larger than the distance traveled ' + \
'inside Earth for this value of costhz')
# Clamp values of l within tolerance of the exit point onto the exit point
l = np.minimum(l, d)
# As in distance_traveled_inside_earth, the default geometry returns through the
# expression this function has always used. Here the two forms are not bitwise equal
# -- they differ by 1.5e-14 relative, the general one being the more accurate, since
# it is exactly symmetric about the closest approach where this one is symmetric only
# to 8.8e-10 km -- so routing the default through the general form would perturb every
# Earth result the package has ever produced. Improving the default is a change worth
# making on its own evidence, not a side effect of adding a parameter.
if source_depth == 0.0 and detector_depth == 0.0:
r2 = gd.EARTH_RADIUS*gd.EARTH_RADIUS
r2 = r2 + (d-l)**2
r2 = r2 + 2.0*gd.EARTH_RADIUS*(d-l)*costhz
r = np.sqrt(np.abs(r2))
return float(r) if scalar_input else r
_, r_d = _validated_endpoint_radii(costhz, source_depth, detector_depth,
'earth_radial_distance_from_depth')
# Distance from the entry point to the trajectory's closest approach to the center,
# which is where its radius is the impact parameter b. Everything else follows from
# Pythagoras on the right triangle (b, l - l_closest, r).
b = r_d*np.sqrt(1.0 - costhz*costhz)
l_closest = d + r_d*costhz
r = np.sqrt(b*b + (l - l_closest)**2)
return float(r) if scalar_input else r
[docs]
def prem_layer_edges_along_chord(costhz: float, source_depth: Optional[float]=0.0,
detector_depth: Optional[float]=0.0) -> np.ndarray:
r"""Returns the positions along a chord through the Earth at which
the chord crosses the PREM layer boundaries.
A neutrino entering the Earth with direction ``costhz`` travels along
a chord from :math:`l = 0` to
:math:`l =` :func:`distance_traveled_inside_earth` (``costhz``).
The matter density along the chord is piecewise-smooth, with
discontinuities (or kinks) where the chord crosses the boundaries
between PREM shells. This routine returns those crossing positions,
which are useful as mandatory slab edges for the Magnus expansion:
high-order quadrature converges at its nominal order only if the
Hamiltonian is smooth within each slab.
The crossing positions solve :math:`r(l) = r_b` for each boundary
radius :math:`r_b`, which is a quadratic equation in :math:`l`: with
:math:`u = d - l` and :math:`d = -2 R \cos\theta_z`, one has
.. math::
u^2 + 2 R \cos\theta_z\, u + \left(R^2 - r_b^2\right) = 0 .
The two depths move the endpoints of the trajectory, exactly as they
do in :func:`distance_traveled_inside_earth`. Only crossings strictly
inside the trajectory are returned, so burying an endpoint drops the
boundaries the shortened path no longer reaches. An endpoint that
lands exactly on a boundary radius is not a crossing: the density is
smooth on the whole of a path that stops there.
.. versionadded:: 1.0.0
.. versionchanged:: 1.1.1
Added ``source_depth`` and ``detector_depth``. Their defaults of
zero reproduce the surface-to-surface chord bit for bit, through
the same expression as before.
Parameters
----------
costhz : float
Cosine of the zenith angle of the neutrino, measured at the
detector. With both endpoints on the surface, crossings exist
only for costhz < 0; a buried detector also sees a downward-going
trajectory through its overburden.
source_depth : float, optional
Depth of the entry point below the surface [km]. Default: 0.0.
detector_depth : float, optional
Depth of the detector below the surface [km]. Default: 0.0.
Returns
-------
np.ndarray
Sorted crossing positions l [km], each strictly inside (0, d).
Empty if the chord crosses no boundary.
Examples
--------
.. jupyter-execute::
import numpy as np
from magnus import earth
edges = earth.prem_layer_edges_along_chord(-0.8)
d = earth.distance_traveled_inside_earth(-0.8)
print('%d crossings; the first three at %s km'
% (len(edges), np.round(edges[:3], 1)))
print('symmetric about the midpoint:',
np.allclose(edges + edges[::-1], d))
A buried detector loses the crossings its shortened path no longer
reaches. Looking straight up, it keeps all eighteen until it passes
below a boundary itself; looking straight down, its whole overburden
lies inside PREM's outermost shell, so it crosses nothing:
.. jupyter-execute::
for depth in (0.0, 2.0, 20.0):
print('%5.1f km down: %2d crossings looking up, %d looking down'
% (depth,
len(earth.prem_layer_edges_along_chord(-1.0, detector_depth=depth)),
len(earth.prem_layer_edges_along_chord(1.0, detector_depth=depth))))
"""
source_depth, detector_depth = _depths_or_zero(source_depth, detector_depth)
# The default geometry keeps the expression this function has always used, for the
# reason given in distance_traveled_inside_earth. The general branch below reduces to
# it, the two differing only in which end each parameterizes from.
if source_depth == 0.0 and detector_depth == 0.0:
if costhz >= 0.0:
return np.array([])
R = gd.EARTH_RADIUS
d = -2.0*R*costhz # chord length [km]
rmin2 = R*R*(1.0 - costhz*costhz) # (squared) closest approach to the center
crossings = []
for rb in PREM_BOUNDARIES:
disc = rb*rb - rmin2
if disc <= 0.0: # chord never reaches this depth
continue
s = np.sqrt(disc)
for u in (-R*costhz - s, -R*costhz + s):
if 0.0 < u < d:
crossings.append(d - u)
return np.unique(np.array(sorted(crossings)))
_, r_d = _validated_endpoint_radii(costhz, source_depth, detector_depth,
'prem_layer_edges_along_chord')
d = distance_traveled_inside_earth(costhz, source_depth, detector_depth)
if d <= 0.0: # a trajectory of no length crosses nothing
return np.array([])
rmin2 = r_d*r_d*(1.0 - costhz*costhz) # (squared) closest approach to the center
l_closest = d + r_d*costhz # entry point to closest approach [km]
crossings = []
for rb in PREM_BOUNDARIES:
disc = rb*rb - rmin2
if disc <= 0.0: # trajectory never reaches this depth
continue
s = np.sqrt(disc)
for l in (l_closest - s, l_closest + s):
if 0.0 < l < d:
crossings.append(l)
return np.unique(np.array(sorted(crossings)))
[docs]
def dms_to_decimal(degrees: float, minutes: float, seconds: float) -> float:
r"""Converts (degree, minute, second) coordinates to decimal degrees.
A West longitude or a South latitude is negative. The sign is read from
the first non-zero part, and the minutes and seconds are magnitudes that
count away from zero in that direction, so ``(-88, 15, 26)`` and
``(-88, -15, -26)`` both give :math:`-88.257^\circ`. A coordinate within
a degree of the meridian or the equator carries its sign on the minutes,
``(0, -30, 0)``, or on the degrees as ``-0.0``.
.. versionadded:: 1.0.0
.. versionchanged:: 1.1.1
The minutes and seconds now follow the sign of the coordinate. They
used to be added as given, so ``(-88, 15, 26)``, the form the built-in
location table uses, came out :math:`-87.743^\circ` and moved every
West or South site toward zero by up to one degree: the chord from
Fermilab to Homestake was 1207 km instead of 1285 km.
Parameters
----------
degrees : float
Degree part of the coordinate.
minutes : float
Minute part of the coordinate.
seconds : float
Second part of the coordinate.
Returns
-------
float
Coordinate in decimal degrees.
"""
sign = 1.0
for part in (degrees, minutes, seconds):
if part != 0 or np.copysign(1.0, part) < 0:
sign = np.copysign(1.0, part)
break
return sign*(abs(degrees) + abs(minutes) / 60 + abs(seconds) / 3600)
[docs]
def chord_length_inside_earth(lat1_dms: tuple[float, float, float],
lon1_dms: tuple[float, float, float], lat2_dms: tuple[float, float, float],
lon2_dms: tuple[float, float, float]) -> float:
r"""Returns the chord length between two locations on the surface of
the Earth.
Computes the straight-line (chord) distance between two locations on the surface of the
Earth, assumed spherical, using the haversine formula for the central angle between the two
locations and converting it to a chord length.
.. versionadded:: 1.0.0
Parameters
----------
lat1_dms : tuple of float
Latitude of the first location, as (degrees, minutes, seconds).
lon1_dms : tuple of float
Longitude of the first location, as (degrees, minutes, seconds).
lat2_dms : tuple of float
Latitude of the second location, as (degrees, minutes, seconds).
lon2_dms : tuple of float
Longitude of the second location, as (degrees, minutes, seconds).
Returns
-------
float
Chord length between the two locations [km].
Examples
--------
.. jupyter-execute::
from magnus import earth
fermilab = ((41.0, 49.0, 55.0), (-88.0, -15.0, -26.0))
sanford = ((44.0, 21.0, 12.0), (-103.0, -45.0, -5.0))
print('Fermilab to Sanford: %.1f km'
% earth.chord_length_inside_earth(fermilab[0], fermilab[1],
sanford[0], sanford[1]))
"""
# Convert DMS to decimal degrees
lat1 = dms_to_decimal(*lat1_dms)
lon1 = dms_to_decimal(*lon1_dms)
lat2 = dms_to_decimal(*lat2_dms)
lon2 = dms_to_decimal(*lon2_dms)
# Convert decimal degrees to radians
lat1_rad = np.radians(lat1)
lon1_rad = np.radians(lon1)
lat2_rad = np.radians(lat2)
lon2_rad = np.radians(lon2)
# Differences in coordinates
delta_lat = lat2_rad - lat1_rad
delta_lon = lon2_rad - lon1_rad
# Haversine formula to calculate the central angle
a = np.sin(delta_lat / 2)**2 + np.cos(lat1_rad) * np.cos(lat2_rad) * np.sin(delta_lon / 2)**2
central_angle = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
# Straight-line distance (chord length)
distance = 2 * gd.EARTH_RADIUS * np.sin(central_angle / 2)
return distance
[docs]
def costhz_between_points_on_surface(lat1_dms: tuple[float, float, float],
lon1_dms: tuple[float, float, float], lat2_dms: tuple[float, float, float],
lon2_dms: tuple[float, float, float]) -> float:
r"""Returns the zenith angle of the chord between two locations on
the surface of the Earth.
Computes the cosine of the zenith angle at which a neutrino would need to travel in a
straight chord through the Earth to reach the second location from the first (e.g., a source
and a detector both on the surface). Assumes a spherical Earth and a detector on the
surface, not underground, so the returned value is always non-positive. A neutrino
arriving from above, ``costhz > 0``, crosses no part of the Earth's interior, and
``costhz = 0`` grazes the surface horizontally.
Two surface coordinates cannot describe a buried endpoint; to place one, give the zenith
angle directly and pass ``source_depth`` or ``detector_depth`` to
:func:`distance_traveled_inside_earth`.
.. versionadded:: 1.0.0
Parameters
----------
lat1_dms : tuple of float
Latitude of the first location, as (degrees, minutes, seconds).
lon1_dms : tuple of float
Longitude of the first location, as (degrees, minutes, seconds).
lat2_dms : tuple of float
Latitude of the second location, as (degrees, minutes, seconds).
lon2_dms : tuple of float
Longitude of the second location, as (degrees, minutes, seconds).
Returns
-------
float
Cosine of the zenith angle of the chord connecting the two locations.
"""
chord_length = chord_length_inside_earth(lat1_dms, lon1_dms, lat2_dms, lon2_dms) # [km]
return -0.5 * chord_length / gd.EARTH_RADIUS
[docs]
def coordinates_of_named_location(source_func_name: str, loc_name: str) -> np.ndarray:
r"""Returns the coordinates of a predefined location (e.g., a
neutrino detector site).
Looks up ``loc_name`` (case-insensitively, spaces treated as underscores) in the
``loc_coords_dms`` dictionary of predefined locations (neutrino telescopes/detector sites and
a few reference points) and returns its latitude and longitude.
.. versionadded:: 1.0.0
Parameters
----------
source_func_name : str
Name of the calling function, used only to build a more informative error message if
``loc_name`` is not found. The message prefixes it with ``oscprob.``, so a caller
from another module is reported under that name.
loc_name : str
Name of the predefined location (e.g., ``'kamioka'``, ``'south_pole'``). See
``earth.loc_coords_dms`` for the full list.
Returns
-------
np.ndarray
Array of shape ``(2, 3)`` in degrees, arcminutes and arcseconds: row 0 the
latitude, row 1 the longitude. The stored tuples are converted to floats.
Raises
------
ValueError
If ``loc_name`` is not one of the predefined locations.
"""
# The latitude and longitude are each returned in day-minute-second format, (dd, mm, ss)
try:
lat = loc_coords_dms[loc_name.lower().replace(" ", "_")]['lat']
lon = loc_coords_dms[loc_name.lower().replace(" ", "_")]['lon']
except KeyError:
raise ValueError(gd.ERROR_MSG_NO_COLOR + " oscprob." + source_func_name + ": the given name of the" + \
" location (" + loc_name + ") is not one of the predefined named locations" + \
" in Magnus. The available predefined named locations (in" + \
" earth.loc_coords_dms)" + " are: " + str(list(loc_coords_dms.keys())) + ".")
return np.array([lat, lon])
# ---------------------------------------------------------------------------------------
# Composition: the electron fraction Y_e = <Z/A>, by PREM layer.
#
# PREM is a *density* model and carries no composition, so Y_e has to be supplied. The
# library used to assume 0.5 everywhere -- exactly isoscalar matter, which nothing in the
# Earth is -- and that is worth up to a factor of ten in P(nu_mu -> nu_e) on a
# core-crossing chord, because the core is iron.
#
# Each value is <Z/A> of the material:
#
# core 0.4656 pure iron, Z/A = 26/55.845; Ni pulls up, light elements down
# mantle 0.4957 peridotite; O and Si sit at ~0.5, a few per cent of Fe pulls it down
# crust 0.4952 granitic continental crust -- 0.0005 from the mantle, i.e. nothing
# ocean 0.5551 seawater, H2O; hydrogen has Z/A = 1, which is why this is ABOVE 0.5
#
# The crust value is here for explicitness rather than for effect: it differs from the
# mantle by 0.1%, far inside PREM's own density uncertainty. The two splits that move a
# number are core-vs-rest (6%) and ocean-vs-rest (12%).
#
# The ocean is a caveat rather than a fact about a given baseline. PREM's outermost 3 km
# is a global-average ocean; a detector under rock has none, and for a trajectory within
# ~2.3 degrees of horizontal that shell is the *entire* path. Pass
# `electron_fraction_ocean=Y_E_CRUST_PREM` for a land baseline -- PREM cannot know which
# you mean.
r"""float: Module-level constant
Electron fraction :math:`Y_e = \langle Z/A \rangle` of the Earth's core (:math:`r \le
3480` km), taken as pure iron. Units: [Adimensional]
.. versionadded:: 1.0.0
"""
[docs]
Y_E_MANTLE_PREM = 0.4957
r"""float: Module-level constant
Electron fraction of the mantle (:math:`3480 < r \le 6346.6` km), peridotite.
Units: [Adimensional]
.. versionadded:: 1.0.0
"""
[docs]
Y_E_CRUST_PREM = 0.4952
r"""float: Module-level constant
Electron fraction of the crust (:math:`6346.6 < r \le 6368` km), granitic. Within
0.1% of the mantle; separate for explicitness rather than for effect.
Units: [Adimensional]
.. versionadded:: 1.0.0
"""
[docs]
Y_E_OCEAN_PREM = 0.5551
r"""float: Module-level constant
Electron fraction of PREM's ocean layer (:math:`r > 6368` km), seawater. Above 0.5
because hydrogen has :math:`Z/A = 1`. Units: [Adimensional]
.. versionadded:: 1.0.0
"""
# The two radii that separate the four compositions. Both are already PREM shell
# boundaries (PREM_BOUNDARIES[1] and [6], [8]), so this introduces no new constant.
_CORE_MANTLE_BOUNDARY = 3480.0
_MANTLE_CRUST_BOUNDARY = 6346.6
_CRUST_OCEAN_BOUNDARY = 6368.0
[docs]
def electron_fraction_func_prem(
r,
electron_fraction_core=None,
electron_fraction_mantle=None,
electron_fraction_crust=None,
electron_fraction_ocean=None,
):
r"""Electron fraction :math:`Y_e` at one or more radii, by PREM layer.
.. versionadded:: 1.0.0
Parameters
----------
r : float or np.ndarray
Radial distance from the Earth's center [km].
electron_fraction_core : float, optional
:math:`Y_e` for :math:`r \le 3480` km. Default: :data:`Y_E_CORE_PREM`.
electron_fraction_mantle : float, optional
:math:`Y_e` for :math:`3480 < r \le 6346.6` km. Default: :data:`Y_E_MANTLE_PREM`.
electron_fraction_crust : float, optional
:math:`Y_e` for :math:`6346.6 < r \le 6368` km. Default: :data:`Y_E_CRUST_PREM`.
electron_fraction_ocean : float, optional
:math:`Y_e` for :math:`r > 6368` km. Default: :data:`Y_E_OCEAN_PREM`.
Returns
-------
np.ndarray
:math:`Y_e` at each radius, with the shape of ``r``.
"""
core = Y_E_CORE_PREM if electron_fraction_core is None else float(electron_fraction_core)
mantle = (Y_E_MANTLE_PREM if electron_fraction_mantle is None
else float(electron_fraction_mantle))
crust = Y_E_CRUST_PREM if electron_fraction_crust is None else float(electron_fraction_crust)
ocean = Y_E_OCEAN_PREM if electron_fraction_ocean is None else float(electron_fraction_ocean)
rr = np.asarray(r, dtype=float)
out = np.full(rr.shape, mantle, dtype=float)
out = np.where(rr <= _CORE_MANTLE_BOUNDARY, core, out)
out = np.where(rr > _MANTLE_CRUST_BOUNDARY, crust, out)
out = np.where(rr > _CRUST_OCEAN_BOUNDARY, ocean, out)
return out
[docs]
def neutron_to_proton_ratio_from_electron_fraction(electron_fraction):
r"""The neutron-to-proton ratio implied by an electron fraction.
The two are not independent. With charge neutrality :math:`n_p = n_e = Y_e
n_{\rm nucleon}` and :math:`n_n = (1 - Y_e) n_{\rm nucleon}`, so
.. math::
r = \frac{n_n}{n_p} = \frac{1 - Y_e}{Y_e} .
Deriving one from the other is what keeps a medium physical: the two used to be
independent arguments, so setting :math:`Y_e = 0.4656` for an iron core while leaving
:math:`r` at its isoscalar default of 1.0 described matter that cannot exist -- and
silently, since :math:`r` only shows up in the sterile sector.
.. versionadded:: 1.0.0
Parameters
----------
electron_fraction : float or np.ndarray
:math:`Y_e`, in (0, 1]. Unchecked here: the domain is enforced by the caller,
``magnus.oscprob._earth_composition``. A zero returns ``inf`` with a NumPy
divide warning, and a value outside the range returns a negative ratio in
silence.
Returns
-------
np.ndarray
:math:`r = n_n/n_p`, with the shape of the input.
"""
ye = np.asarray(electron_fraction, dtype=float)
return (1.0 - ye)/ye
__all__ = [
'Y_E_CORE_PREM',
'Y_E_MANTLE_PREM',
'Y_E_CRUST_PREM',
'Y_E_OCEAN_PREM',
'electron_fraction_func_prem',
'neutron_to_proton_ratio_from_electron_fraction',
'loc_coords_dms',
'PREM_BOUNDARIES',
'density_matter_func_prem',
'distance_traveled_inside_earth',
'earth_radial_distance_from_depth',
'prem_layer_edges_along_chord',
'dms_to_decimal',
'chord_length_inside_earth',
'costhz_between_points_on_surface',
'coordinates_of_named_location',
]