magnus.earth

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

Attributes

loc_coords_dms

PREM_BOUNDARIES

Y_E_CORE_PREM

Module-level constant

Y_E_MANTLE_PREM

Module-level constant

Y_E_CRUST_PREM

Module-level constant

Y_E_OCEAN_PREM

Module-level constant

Functions

density_matter_func_prem(→ Union[float, numpy.ndarray])

Returns the matter density inside the Earth according to the

distance_traveled_inside_earth(→ float)

Returns the distance traveled by a neutrino inside the Earth,

earth_radial_distance_from_depth(→ Union[float, ...)

Returns the radial distance measured from the center of the

prem_layer_edges_along_chord(→ numpy.ndarray)

Returns the positions along a chord through the Earth at which

dms_to_decimal(→ float)

Converts (degree, minute, second) coordinates to decimal degrees.

chord_length_inside_earth(→ float)

Returns the chord length between two locations on the surface of

costhz_between_points_on_surface(→ float)

Returns the zenith angle of the chord between two locations on

coordinates_of_named_location(→ numpy.ndarray)

Returns the coordinates of a predefined location (e.g., a

electron_fraction_func_prem(r[, ...])

Electron fraction \(Y_e\) at one or more radii, by PREM layer.

neutron_to_proton_ratio_from_electron_fraction(...)

The neutron-to-proton ratio implied by an electron fraction.

Module Contents

magnus.earth.loc_coords_dms[source]
magnus.earth.PREM_BOUNDARIES[source]
magnus.earth.density_matter_func_prem(r: float | numpy.ndarray, tol: float | None = 1e-08, density_matter_ocean: float | None = None) → float | numpy.ndarray[source]

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.

Added in version 1.0.0.

Changed in version 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, \(r > 6368\) km [\(\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 electron_fraction_func_prem() handles. Default: None, i.e. PREM’s own ocean.

Returns:

Matter density [\(\text{g cm}^{-3}\)].

Return type:

float or np.ndarray

Raises:

ValueError – If any radial distance exceeds globaldefs.EARTH_RADIUS by more than the relative tolerance tol.

References

Examples

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)))

r =      0 km ->  13.09 g/cm^3
r =   3000 km ->  10.60 g/cm^3
r =   5000 km ->   4.79 g/cm^3
r =   6371 km ->   1.02 g/cm^3
magnus.earth.distance_traveled_inside_earth(costhz: float, source_depth: float | None = 0.0, detector_depth: float | None = 0.0) → float[source]

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 \(r_{\rm s} = R_\oplus - {}\) source_depth and \(r_{\rm d} = R_\oplus - {}\) detector_depth, the impact parameter of the trajectory is \(b = r_{\rm d} \sqrt{1 - \cos^2\theta_z}\) and its length is

\[L = \sqrt{r_{\rm s}^2 - b^2} - r_{\rm d} \cos\theta_z ~.\]

Added in version 1.0.0.

Changed in version 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:

Path length inside the Earth [km]. Every osc_prob_* baseline is in \(\text{eV}^{-1}\), so multiply by magnus.globaldefs.UNIT_KM before passing this on; handing the raw value over returns a converged, unitary, wrong answer.

Return type:

float

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

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)))

costhz = -0.20 ->   2548.4 km
costhz = -0.50 ->   6371.0 km
costhz = -1.00 ->  12742.0 km

A detector 2 km underground sees a shorter upward-going path, and a downward-going one through its overburden:

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)))

costhz = -1.00 -> 12740.0000 km
costhz =  0.50 ->     3.9981 km
costhz =  1.00 ->     2.0000 km
magnus.earth.earth_radial_distance_from_depth(costhz: float, l: float | numpy.ndarray, tol: float | None = 1e-08, source_depth: float | None = 0.0, detector_depth: float | None = 0.0) → float | numpy.ndarray[source]

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 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.

Added in version 1.0.0.

Changed in version 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:

Radial distance to the neutrino [km].

Return type:

float or np.ndarray

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 distance_traveled_inside_earth().

Examples

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)))

l =     0.000 km -> r = 6371.000 km
l =  5096.450 km -> r = 3821.400 km
l = 10192.900 km -> r = 6369.000 km
magnus.earth.prem_layer_edges_along_chord(costhz: float, source_depth: float | None = 0.0, detector_depth: float | None = 0.0) → numpy.ndarray[source]

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 \(l = 0\) to \(l =\) 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 \(r(l) = r_b\) for each boundary radius \(r_b\), which is a quadratic equation in \(l\): with \(u = d - l\) and \(d = -2 R \cos\theta_z\), one has

\[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 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.

Added in version 1.0.0.

Changed in version 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:

Sorted crossing positions l [km], each strictly inside (0, d). Empty if the chord crosses no boundary.

Return type:

np.ndarray

Examples

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))

14 crossings; the first three at [ 3.8 18.8 30.5] km
symmetric about the midpoint: True

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:

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))))

  0.0 km down: 18 crossings looking up, 0 looking down
  2.0 km down: 18 crossings looking up, 0 looking down
 20.0 km down: 16 crossings looking up, 2 looking down
magnus.earth.dms_to_decimal(degrees: float, minutes: float, seconds: float) → float[source]

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 \(-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.

Added in version 1.0.0.

Changed in version 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 \(-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:

Coordinate in decimal degrees.

Return type:

float

magnus.earth.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[source]

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.

Added in version 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:

Chord length between the two locations [km].

Return type:

float

Examples

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]))

Fermilab to Sanford: 1284.7 km
magnus.earth.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[source]

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 distance_traveled_inside_earth().

Added in version 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:

Cosine of the zenith angle of the chord connecting the two locations.

Return type:

float

magnus.earth.coordinates_of_named_location(source_func_name: str, loc_name: str) → numpy.ndarray[source]

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.

Added in version 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:

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.

Return type:

np.ndarray

Raises:

ValueError – If loc_name is not one of the predefined locations.

magnus.earth.Y_E_CORE_PREM = 0.4656[source]

Module-level constant

Electron fraction \(Y_e = \langle Z/A \rangle\) of the Earth’s core (\(r \le 3480\) km), taken as pure iron. Units: [Adimensional]

Added in version 1.0.0.

Type:

float

magnus.earth.Y_E_MANTLE_PREM = 0.4957[source]

Module-level constant

Electron fraction of the mantle (\(3480 < r \le 6346.6\) km), peridotite. Units: [Adimensional]

Added in version 1.0.0.

Type:

float

magnus.earth.Y_E_CRUST_PREM = 0.4952[source]

Module-level constant

Electron fraction of the crust (\(6346.6 < r \le 6368\) km), granitic. Within 0.1% of the mantle; separate for explicitness rather than for effect. Units: [Adimensional]

Added in version 1.0.0.

Type:

float

magnus.earth.Y_E_OCEAN_PREM = 0.5551[source]

Module-level constant

Electron fraction of PREM’s ocean layer (\(r > 6368\) km), seawater. Above 0.5 because hydrogen has \(Z/A = 1\). Units: [Adimensional]

Added in version 1.0.0.

Type:

float

magnus.earth.electron_fraction_func_prem(r, electron_fraction_core=None, electron_fraction_mantle=None, electron_fraction_crust=None, electron_fraction_ocean=None)[source]

Electron fraction \(Y_e\) at one or more radii, by PREM layer.

Added in version 1.0.0.

Parameters:

  • r (float or np.ndarray) – Radial distance from the Earth’s center [km].

  • electron_fraction_core (float, optional) – \(Y_e\) for \(r \le 3480\) km. Default: Y_E_CORE_PREM.

  • electron_fraction_mantle (float, optional) – \(Y_e\) for \(3480 < r \le 6346.6\) km. Default: Y_E_MANTLE_PREM.

  • electron_fraction_crust (float, optional) – \(Y_e\) for \(6346.6 < r \le 6368\) km. Default: Y_E_CRUST_PREM.

  • electron_fraction_ocean (float, optional) – \(Y_e\) for \(r > 6368\) km. Default: Y_E_OCEAN_PREM.

Returns:

\(Y_e\) at each radius, with the shape of r.

Return type:

np.ndarray

magnus.earth.neutron_to_proton_ratio_from_electron_fraction(electron_fraction)[source]

The neutron-to-proton ratio implied by an electron fraction.

The two are not independent. With charge neutrality \(n_p = n_e = Y_e n_{\rm nucleon}\) and \(n_n = (1 - Y_e) n_{\rm nucleon}\), so

\[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 \(Y_e = 0.4656\) for an iron core while leaving \(r\) at its isoscalar default of 1.0 described matter that cannot exist – and silently, since \(r\) only shows up in the sterile sector.

Added in version 1.0.0.

Parameters:

electron_fraction (float or np.ndarray) – \(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:

\(r = n_n/n_p\), with the shape of the input.

Return type:

np.ndarray