API reference

Generated from the docstrings. Every Examples block is executed when this page is built, so what is shown is what the code returns rather than numbers written beside it.

Physics

Closed-form physics with no terrain in it: atmosphere, Earth chord, tau range and exit probability, shower development, geomagnetic field, Cherenkov footprint. Self-contained and usable on its own.

Closed-form physics the terrain scan needs but cannot measure from a DEM.

Everything here is analytic. The scan supplies geometry – where the tau exits, how far away, through how much local rock – and these functions supply the physics that geometry alone does not determine: how much atmosphere a shower has to develop in, how much Earth a neutrino crossed to get there, how strongly the shower radiates into a given direction, and how large its radio footprint is when it arrives.

Constants that are genuinely uncertain are parameters with stated defaults rather than literals buried in an expression. Where a default encodes a convention or an approximation, the docstring says so.

oroscope.physics.air_density_kgm3(altitude_m: float, sea_level_density: float = 1.225, scale_height_m: float = 8400.0) → float

Exponential atmosphere, rho0 * exp(-h/H).

An isothermal approximation, adequate over the few kilometres of relief a site search spans. It is not a substitute for a real profile at large zenith angles.

Parameters:

altitude_mfloat

Altitude above sea level, in metres.

sea_level_densityfloat, optional

Density at sea level, in kg/m^3.

scale_height_mfloat, optional

Density scale height, in metres.

Returns:

float

Air density in kg/m^3.

Examples

>>> from oroscope import physics
>>> round(physics.air_density_kgm3(0.0), 3)
1.225
>>> round(physics.air_density_kgm3(4000.0), 3)   # a third thinner at Andean altitude
0.761

oroscope.physics.slant_grammage_gcm2(start_altitude_m: float, elevation_deg: float, distance_m: float, sea_level_density: float = 1.225, scale_height_m: float = 8400.0) → float

Atmospheric depth along a slanted path, in g/cm^2.

A shower develops through grammage, not through metres, and air density falls by a third between 2000 m and 4500 m. A site search that compares candidates at different altitudes while measuring path length in metres is comparing unlike things: 20 km at 4000 m is about 1500 g/cm^2, the same 20 km at sea level about 2450 g/cm^2.

Integrating rho0 exp(-(z0 + l sin(theta))/H) dl along the slant path has a closed form, so no numerical integration is needed:

X = rho0 exp(-z0/H) * H/sin(theta) * (1 - exp(-D sin(theta)/H)) / 1

with the horizontal-path limit rho0 exp(-z0/H) * D / cos(theta) as theta -> 0.

Parameters:

start_altitude_mfloat

Altitude of the near end of the path, in metres.

elevation_degfloat

Elevation angle of the path, in degrees, positive upward.

distance_mfloat

Ground distance covered, in metres. Zero or less returns zero.

sea_level_densityfloat, optional

Density at sea level, in kg/m^3.

scale_height_mfloat, optional

Density scale height, in metres.

Returns:

float

Atmospheric depth along the path, in g/cm^2.

Examples

>>> from oroscope import physics
>>> horizontal = physics.slant_grammage_gcm2(4000.0, 0.0, 20000.0)
>>> sea_level = physics.slant_grammage_gcm2(0.0, 0.0, 20000.0)
>>> f"{horizontal:.0f} vs {sea_level:.0f} g/cm^2"
'1522 vs 2450 g/cm^2'

oroscope.physics.shower_maturity(grammage_gcm2: float | ndarray, x_max_gcm2: float = 700.0) → float | ndarray

Path grammage as a fraction of the depth of shower maximum.

Below 1 the shower is still developing when it arrives. What “above 1” means depends on what is being detected, and the two cases are not alike:

Radio. Emission comes from the region around shower maximum and then simply propagates; air is effectively transparent at 50-200 MHz. Being far beyond maximum costs nothing directly, so the criterion is a threshold, not a band. The real trade at greater distance is amplitude against footprint area, which belongs to the footprint term rather than here.

Particles. The charged-particle content peaks at maximum and dies away after, so a particle array such as TAMBO does want to sit near it, and there the criterion genuinely is a band.

This is one more reason criteria have to be per-channel rather than global.

Parameters:

grammage_gcm2float or array_like

Atmospheric depth traversed, in g/cm^2.

x_max_gcm2float, optional

Depth of shower maximum, in g/cm^2.

Returns:

float or ndarray

Grammage as a fraction of shower maximum. Below 1 the shower is still growing.

See also

shower_size_fraction

the particle content itself, rather than this ratio.

Examples

>>> from oroscope import physics
>>> round(physics.shower_maturity(750.0), 3)
1.071

oroscope.physics.shower_maximum_gcm2(energy_pev: float | ndarray, x_max_ref_gcm2: float = 700.0, reference_energy_pev: float = 1000.0, elongation_rate: float = 55.0) → ndarray

Depth of shower maximum at a given primary energy.

X_max(E) = X_max(E_ref) + D * log10(E / E_ref)

Over TAMBO’s 3 PeV to 1 EeV this runs from about 560 to 700 g/cm^2, so the energy dependence is real but mild – the band below is set far more by how much of the profile is being accepted than by where its peak sits.

Parameters:

energy_pevfloat or array_like

Primary energy, in PeV.

x_max_ref_gcm2float, optional

Depth of maximum at the reference energy, in g/cm^2.

reference_energy_pevfloat, optional

Energy at which x_max_ref_gcm2 is quoted, in PeV.

elongation_ratefloat, optional

Deepening per decade of energy, in g/cm^2. About 55 for a hadronic cascade; a purely electromagnetic one is nearer the Heitler value of 85.

Returns:

ndarray

Depth of shower maximum, in g/cm^2.

Examples

>>> from oroscope import physics
>>> [f"{float(physics.shower_maximum_gcm2(e)):.0f}" for e in (3.0, 1000.0)]
['561', '700']

oroscope.physics.shower_size_fraction(grammage_gcm2: float | ndarray, x_max_gcm2: float, lambda_gcm2: float = 70.0) → ndarray

Charged-particle content at depth X, as a fraction of the content at maximum.

Gaisser-Hillas with the first interaction at X_0 = 0:

N(X)/N_max = (X / X_max)^(X_max/lambda) * exp((X_max - X) / lambda)

Zero at and below X = 0. This is what makes a particle array’s criterion a band rather than a threshold: the content rises steeply, peaks, and then dies.

Parameters:

grammage_gcm2float or array_like

Depth at which to evaluate the profile, in g/cm^2.

x_max_gcm2float

Depth of shower maximum, in g/cm^2.

lambda_gcm2float, optional

Gaisser-Hillas interaction length, in g/cm^2, setting how fast the profile rises and falls.

Returns:

ndarray

Particle content as a fraction of the content at maximum, in [0, 1].

Examples

>>> from oroscope import physics
>>> round(float(physics.shower_size_fraction(700.0, 700.0)), 3)   # at maximum
1.0
>>> round(float(physics.shower_size_fraction(172.0, 561.0)), 3)   # a 2 km crossing
0.02

oroscope.physics.grammage_band_from_energy(energy_min_pev: float, energy_max_pev: float, fraction: float = 0.1, lambda_gcm2: float = 70.0, samples: int = 4000, **x_max_kw: float) → tuple[float, float]

Atmospheric-depth band over which a particle array still sees a usable shower.

A particle detector needs the cascade to have grown but not yet died, so the criterion is the depth range where the charged-particle content stays above fraction of its own maximum. The band is taken across the requested energy range: the low edge at the lowest energy, whose maximum is shallowest, and the high edge at the highest, whose maximum is deepest.

At TAMBO’s 3 PeV to 1 EeV and fraction 0.1 this gives roughly 235 to 1300 g/cm^2. That matters for siting: a canyon crossing supplies only what its own width of air contains – about 170 g/cm^2 across 2 km of Colca, and ~390 g/cm^2 across its full 4.5 km rim to rim – so the criterion selects the widest crossings, which is the physically right answer and not one the default (X_max, 4*X_max) band could express.

Parameters:

energy_min_pev, energy_max_pevfloat

Ends of the primary-energy range, in PeV.

fractionfloat, optional

Fraction of peak particle content that still counts as a usable shower. A choice about detector capability rather than a property of the shower, and one of the parameters a result is most sensitive to.

lambda_gcm2float, optional

Gaisser-Hillas interaction length, in g/cm^2.

samplesint, optional

Points on the depth grid searched for the band edges.

**x_max_kw

Passed through to shower_maximum_gcm2().

Returns:

tuple of float

(low_gcm2, high_gcm2), the low edge taken at the lowest energy and the high edge at the highest, so the band spans the whole requested range.

Examples

>>> from oroscope import physics
>>> lo, hi = physics.grammage_band_from_energy(3.0, 1000.0)
>>> f"{lo:.0f} - {hi:.0f} g/cm^2"
'236 - 1287 g/cm^2'

oroscope.physics.earth_chord_m(elevation_deg: float, radius_m: float = 6371000.0) → float

Chord length through the Earth for a ray arriving from below the horizontal.

A ray making angle theta below the local horizontal cuts a chord of 2R sin(theta): zero along the tangent, a full diameter straight down. Directions at or above the horizontal return zero.

This matters because it dwarfs local topography. At -1 degree the chord is about 220 km, at -3 degrees about 670 km, against the tens of km of mountain a DEM can see. The deepest part of even a 670 km chord lies only ~9 km below the surface, so it stays in the crust and a constant density is adequate.

Parameters:

elevation_degfloat

Arrival elevation angle, in degrees. Zero or above returns zero.

radius_mfloat, optional

Earth radius, in metres. The true radius: the neutrino is not refracted.

Returns:

float

Chord length through the Earth, in metres.

Examples

>>> from oroscope import physics
>>> f"{physics.earth_chord_m(-3.0) / 1000:.0f} km"
'667 km'

oroscope.physics.earth_chord_gcm2(elevation_deg: float, radius_m: float = 6371000.0, density_gcm3: float = 2.65) → float

Column depth of the Earth chord, in g/cm^2.

Parameters:

elevation_degfloat

Arrival elevation angle, in degrees.

radius_mfloat, optional

Earth radius, in metres.

density_gcm3float, optional

Crust density, in g/cm^3.

Returns:

float

Column depth of the chord, in g/cm^2.

Examples

>>> from oroscope import physics
>>> f"{physics.earth_chord_gcm2(-3.0):.2e}"
'1.77e+08'

oroscope.physics.neutrino_survival(elevation_deg: float, interaction_length_gcm2: float, radius_m: float = 6371000.0, density_gcm3: float = 2.65, nc_regeneration: bool = False, spectral_index: float = 2.0) → float

Fraction of neutrinos surviving the Earth chord to reach the exit region.

exp(-X_chord / X_int). The interaction length is a parameter because it depends on the cross-section at the energy of interest; around an EeV it is of order 1e8 g/cm^2, which is the same order as the chord at -3 degrees, so the suppression across a +/-3 degree window is substantial rather than marginal.

With nc_regeneration=True the result is multiplied by nc_regeneration_factor(), and clipped at 1 so that regeneration can offset absorption but never invent flux. Off by default, because it is an approximation and every published number here was computed without it.

Parameters:

elevation_degfloat

Arrival elevation angle, in degrees.

interaction_length_gcm2float

Neutrino interaction length at the energy of interest, in g/cm^2. Zero or less disables the attenuation and returns 1.

radius_mfloat, optional

Earth radius, in metres.

density_gcm3float, optional

Crust density, in g/cm^3.

nc_regenerationbool, optional

Apply the leading-order neutral-current regeneration correction.

spectral_indexfloat, optional

Power-law index used by that correction.

Returns:

float

Surviving fraction, in [0, 1].

Examples

>>> from oroscope import physics
>>> lam = physics.neutrino_interaction_length_gcm2(1000.0)
>>> round(physics.neutrino_survival(-1.0, lam), 3)
0.7

Regeneration lifts it, because absorption alone was overstating the suppression:

>>> round(physics.neutrino_survival(-1.0, lam, nc_regeneration=True), 3)
0.778

oroscope.physics.muon_shielding_gcm2(thickness_km: float, density_gcm3: float = 2.65) → float

Column depth corresponding to a rock thickness, for muon rejection.

A mountain in the arrival direction is a muon filter: an air-shower muon from that direction would have to cross the whole thickness, which a few km of rock makes impossible. Anything detected from behind that much rock is therefore not a cosmic-ray muon, which is what lets a surface array claim neutrino purity.

Unlike the production-and-escape band, this is a floor: more rock is always better for background rejection, and only the signal side wants an upper limit.

4 km of standard rock is about 1.06e6 g/cm^2.

Parameters:

thickness_kmfloat

Rock thickness along the arrival direction, in km.

density_gcm3float, optional

Rock density, in g/cm^3.

Returns:

float

Column depth, in g/cm^2.

Examples

>>> from oroscope import physics
>>> f"{physics.muon_shielding_gcm2(physics.DEFAULT_MUON_SHIELDING_KM):.2e}"
'1.06e+06'

Through the constant rather than a literal 4.0: the thickness is a sourced number, and this function was exported while the answer to “how thick?” was not.

oroscope.physics.tau_decay_length_m(energy_pev: float, mass_gev: float = 1.77686, ctau_m: float = 8.703e-05) → float

Lorentz-boosted tau decay length, (E/m) c*tau.

The quantity that decides whether a tau decays inside a canyon crossing or flies through it: 147 m at 3 PeV against 49 km at 1 EeV, over a gap of a few km.

Parameters:

energy_pevfloat

Tau energy, in PeV.

mass_gevfloat, optional

Tau mass, in GeV.

ctau_mfloat, optional

Proper decay length c * tau, in metres.

Returns:

float

Decay length in the laboratory frame, in metres.

Examples

>>> from oroscope import physics
>>> [f"{physics.tau_decay_length_m(e):.0f}" for e in (3.0, 1000.0)]
['147', '48980']

oroscope.physics.cc_cross_section_cm2(energy_pev: float | ndarray) → float | ndarray

Charged-current neutrino-nucleon cross-section.

A power-law fit to the standard parameterisations, good to tens of per cent over 1e8-1e10 GeV and an extrapolation above that, where no data constrain it.

Parameters:

energy_pevfloat or array_like

Neutrino energy, in PeV.

Returns:

float or ndarray

Cross-section, in cm^2.

Examples

>>> from oroscope import physics
>>> f"{physics.cc_cross_section_cm2(1000.0):.2e}"
'1.01e-32'

oroscope.physics.neutrino_interaction_length_gcm2(energy_pev: float | ndarray) → float | ndarray

Column depth over which a neutrino interacts once, 1/(N_A sigma).

Falls from about 3.8e8 g/cm^2 at 100 PeV to 7e7 at 10 EeV. Only charged-current attenuation is counted; neutral-current regeneration would soften it slightly.

Parameters:

energy_pevfloat or array_like

Neutrino energy, in PeV.

Returns:

float or ndarray

Interaction length, in g/cm^2.

Examples

>>> from oroscope import physics
>>> f"{physics.neutrino_interaction_length_gcm2(100.0):.2e}"
'3.76e+08'

oroscope.physics.nc_regeneration_factor(chord_gcm2: float, cc_length_gcm2: float, spectral_index: float = 2.0, nc_to_cc: float = 0.42, inelasticity: float = 0.25) → float

Leading-order enhancement from neutral-current regeneration. An approximation.

A charged-current interaction removes a neutrino from the beam. A neutral-current one does not – it degrades the energy and the neutrino continues. Counting only CC absorption, as neutrino_survival does by default, therefore understates the flux arriving at a given energy: neutrinos that started higher up the spectrum scatter down into the band and partially refill it.

The size of that refilling follows from the spectrum. For a power law Phi(E) ~ E**-gamma, a neutrino observed at E after one NC interaction of inelasticity y began at E' = E/(1-y). The flux there is larger by (1-y)**gamma, and the Jacobian dE'/dE = 1/(1-y) gives one more power, so each NC scatter contributes (1-y)**(gamma-1) relative to the unscattered flux. With tau_nc = (chord/X_cc) * (sigma_NC/sigma_CC) NC interactions expected along the chord, the leading term is:

1 + tau_nc * (1 - y)**(gamma - 1)

What this is not. It is the first term of a series, not a solution of the cascade equations, and it ignores the nu_tau -> tau -> nu_tau chain that makes the Earth genuinely translucent to tau neutrinos at high energy. It is a correction of the right sign and roughly the right size, and it is capped below (see neutrino_survival()) so it cannot manufacture more flux than arrived. Treat a result that depends strongly on it as a result that needs a real transport code.

Parameters:

chord_gcm2float

Column depth along the chord, in g/cm^2.

cc_length_gcm2float

Charged-current interaction length, in g/cm^2.

spectral_indexfloat, optional

gamma of the assumed power-law flux. A steeper spectrum regenerates less: the neutrinos scattering down into the band come from E/(1-y), above it, and a steeper spectrum has less flux there.

nc_to_ccfloat, optional

sigma_NC / sigma_CC.

inelasticityfloat, optional

Mean fraction of energy lost in one NC interaction.

Returns:

float

Multiplicative enhancement, at least 1.

Examples

A chord one CC interaction length deep, on an E^-2 spectrum:

>>> from oroscope import physics
>>> round(physics.nc_regeneration_factor(1.0e8, 1.0e8), 3)
1.315

A steeper spectrum regenerates less, since less flux sits above the band:

>>> round(physics.nc_regeneration_factor(1.0e8, 1.0e8, spectral_index=2.7), 3)
1.258

A short chord regenerates nothing:

>>> round(physics.nc_regeneration_factor(0.0, 1.0e8), 3)
1.0

oroscope.physics.tau_energy_loss_beta(energy_pev: float, reference: float | None = None, reference_energy_pev: float | None = None, index: float | None = None) → float

Energy-loss coefficient beta(E), rising with energy as photonuclear does.

An estimate, not a fit to published tables: see the module comments for how it was arrived at. It is the least certain number in this module, and it moves the production-and-escape optimum in proportion.

Parameters:

energy_pevfloat

Tau energy, in PeV.

referencefloat, optional

Value of beta at reference_energy_pev, in cm^2/g. Defaults to whatever set_tau_energy_loss() last established, and to BETA_REFERENCE_CM2G if it has not been called.

reference_energy_pevfloat, optional

Energy at which reference applies, in PeV.

indexfloat, optional

Power-law index of the energy dependence. Zero recovers a constant beta.

Returns:

float

Energy-loss coefficient, in cm^2/g.

Examples

>>> from oroscope import physics
>>> f"{physics.tau_energy_loss_beta(100.0):.2e}"
'3.79e-07'

An argument overrides the module setting for that one call:

>>> f"{physics.tau_energy_loss_beta(100.0, reference=1.0e-6, index=0.0):.2e}"
'1.00e-06'

oroscope.physics.tau_range_gcm2(energy_pev: float, beta_cm2g: float | None = None, density_gcm3: float = 2.65) → float

Column depth over which a tau’s survival probability falls to 1/e.

Decay and energy loss couple, because losing energy shortens the boosted decay length. With E(X) = E0 exp(-beta X) the decay probability per unit depth is exp(beta X)/X_decay(E0), and integrating the survival gives

S(X) = exp( -(X_loss/X_decay) (exp(X/X_loss) - 1) ), X_loss = 1/beta

so the 1/e point is

R = X_loss * ln(1 + X_decay/X_loss)

Note this grows logarithmically at high energy rather than saturating at 1/beta. An earlier version of this module combined the two lengths harmonically, which saturates and underestimates the range by a factor of 2 at an EeV and 4 at 10 EeV.

Parameters:

energy_pevfloat

Tau energy on entering the rock, in PeV.

beta_cm2gfloat, optional

Energy-loss coefficient, in cm^2/g. Defaults to tau_energy_loss_beta() at this energy.

density_gcm3float, optional

Rock density, in g/cm^3.

Returns:

float

Column depth at which survival falls to 1/e, in g/cm^2.

Examples

>>> from oroscope import physics
>>> f"{physics.tau_range_gcm2(1000.0) / physics.CRUST_DENSITY_GCM3 / 1e5:.1f} km"
'13.7 km'

oroscope.physics.tau_survival(depth_gcm2: float | ndarray, energy_pev: float, beta_cm2g: float | None = None, density_gcm3: float = 2.65) → ndarray

Probability a tau of the given energy crosses depth_gcm2 of rock without decaying.

The double-exponential form above, which falls far more sharply than a simple exponential once the depth exceeds 1/beta: the tau is losing energy, so its decay length shrinks as it goes.

Parameters:

depth_gcm2float or array_like

Rock traversed, in g/cm^2.

energy_pevfloat

Tau energy on entering the rock, in PeV.

beta_cm2gfloat, optional

Energy-loss coefficient, in cm^2/g.

density_gcm3float, optional

Rock density, in g/cm^3.

Returns:

ndarray

Survival probability, in [0, 1].

Examples

>>> from oroscope import physics
>>> round(float(physics.tau_survival(0.0, 1000.0)), 3)
1.0

oroscope.physics.tau_exit_probability(column_depth_gcm2: float | ndarray, energy_pev: float, beta_cm2g: float | None = None, density_gcm3: float = 2.65, inelasticity: float = 0.2, samples: int = 2000) → float | ndarray

Relative probability that a traversing neutrino yields a tau that escapes.

The neutrino interacts at depth x with probability exp(-x/lambda) dx/lambda; the tau, carrying (1 - y) of the energy, must then cross the remaining X - x and survive:

P(X) = integral_0^X (dx/lambda) exp(-x/lambda) S(X - x)

Evaluated numerically, because S is a double exponential and the integral has no clean closed form. Relative, not absolute: normalisation needs the trigger response (see aperture).

Parameters:

column_depth_gcm2float or array_like

Rock along the arrival direction, in g/cm^2.

energy_pevfloat

Neutrino energy, in PeV.

beta_cm2gfloat, optional

Tau energy-loss coefficient, in cm^2/g.

density_gcm3float, optional

Rock density, in g/cm^3.

inelasticityfloat, optional

Mean inelasticity of the charged-current interaction; the tau carries away 1 - y of the neutrino energy.

samplesint, optional

Points used for the numerical integration over interaction depth, spaced logarithmically in the remaining depth X - x.

Returns:

float or ndarray

Relative exit probability, matching the shape of column_depth_gcm2.

Notes

The grid is in ``u = X - x``, not in ``x``, and that is the whole of the accuracy. Only interactions within a few tau ranges of the far surface contribute anything – everything produced deeper is absorbed – so the integrand is a spike against the far end of a range that can be five decades wide. Sampled uniformly in x, as this was, the spacing outruns the spike and the trapezoid rule reports the area of something it never resolved. Measured at 3 PeV and X = 10^9 g/cm^2:

grid	P(X)
uniform, 2000 pts	8.884e-05
uniform, 20,000	1.328e-05
uniform, 200,000	1.103e-05
uniform, 2,000,000	1.100e-05
log in u, 2000	1.1004e-05

So the old default was 8x the converged value, and the substitution reaches that value with a thousandth of the points. Worse than the magnitude, it inverted the trend. P(X) has a genuine maximum – production_escape_optimum_gcm2() exists to find it, and dP/dX = (S(X) - P(X))/lambda is positive wherever S(X) exceeds P(X) – so it rises to that peak and falls beyond. At 3 PeV the peak sits near 3.1e5 g/cm^2. What the unresolved grid did was invent a second, spurious rise far past the peak, at the edge of the grid, where the converged curve is falling.

The error was confined to low energy against deep rock – at 100 PeV and above the worst case over the same grid was 3%, and below 10^7 g/cm^2 it was exact everywhere – which is why nothing else showed it. But depth_band_from_energy() takes its low edge at the lowest energy asked for, and TAMBO’s configured range starts at 3 PeV, so the band inherited the whole of it: (1.18e8, 2.89e8) over 3 PeV - 1 EeV where the converged answer is (2.17e4, 1.15e8). The published 1 EeV optimum of 5.7e6 g/cm^2 lies inside the second and 20x below the first. No published number moved: every config leaves depth_band_gcm2 null, so no search has ever called this (roadmap 6.44, 6.56).

Examples

>>> from oroscope import physics
>>> p = physics.tau_exit_probability(1.0e6, 1000.0)
>>> 0.0 <= p <= 1.0
True

Converged where it used to be eight times too large:

>>> round(physics.tau_exit_probability(1.0e9, 3.0) * 1.0e5, 3)
1.1

and falling with depth, as more rock must:

>>> p = physics.tau_exit_probability([1.0e8, 3.0e8, 1.0e9, 3.0e9], 3.0)
>>> bool((p[1:] < p[:-1]).all())
True

oroscope.physics.set_tau_energy_loss(reference: float | None = None, reference_energy_pev: float | None = None, index: float | None = None) → dict

Adopts a different beta, in one place, for every function that uses it.

beta is the least certain number in this module – an estimate from mass scaling, in the range (0.4-1.0)e-6 cm^2/g – and the collaboration’s own value should replace it when there is one. Until now that meant editing the source, which is not something a user of an installed package can reasonably do, and which leaves no record of what was used.

This does not change any site-search result. beta enters tau range and survival – production and escape through rock – and the search does not model those: it uses the decay length L = (E/m_tau) c*tau, which is kinematics and carries no beta. So this affects tau_range_gcm2(), tau_survival() and tau_exit_probability(), which the notebooks and the physics page use, and nothing in a search. See ROADMAP 6.38.

Parameters:

referencefloat, optional

beta at reference_energy_pev, in cm^2/g. Left unchanged when None.

reference_energy_pevfloat, optional

Energy at which reference applies, in PeV.

indexfloat, optional

Power-law index of the energy dependence. Zero gives a constant beta.

Returns:

dict

The settings now in force, as tau_energy_loss_settings() returns.

Examples

A collaboration value, adopted once:

>>> from oroscope import physics
>>> _ = physics.set_tau_energy_loss(reference=0.8e-6, index=0.0)
>>> f"{physics.tau_energy_loss_beta(100.0):.2e}"
'8.00e-07'

A constant beta means the same value at every energy:

>>> f"{physics.tau_energy_loss_beta(10000.0):.2e}"
'8.00e-07'

And back to the shipped estimate, which does rise with energy:

>>> _ = physics.restore_tau_energy_loss()
>>> f"{physics.tau_energy_loss_beta(100.0):.2e}"
'3.79e-07'

oroscope.physics.restore_tau_energy_loss() → dict

Restores the shipped estimate, undoing set_tau_energy_loss().

Returns:

dict

The settings now in force.

Examples

>>> from oroscope import physics
>>> physics.restore_tau_energy_loss() == {
...     "reference": physics.BETA_REFERENCE_CM2G,
...     "reference_energy_pev": physics.BETA_REFERENCE_ENERGY_PEV,
...     "index": physics.BETA_ENERGY_INDEX}
True

oroscope.physics.tau_energy_loss_settings() → dict

The beta parameters currently in force.

Returns:

dict

{"reference", "reference_energy_pev", "index"}, in cm^2/g and PeV.

Examples

>>> from oroscope import physics
>>> physics.tau_energy_loss_settings()["reference"]
6e-07

oroscope.physics.production_escape_optimum_gcm2(energy_pev: float, beta_cm2g: float | None = None, density_gcm3: float = 2.65, inelasticity: float = 0.2, samples: int = 400) → float

Column depth maximising tau_exit_probability().

Found on a log grid, since the corrected exit probability has no closed-form peak. Rises with energy and then flattens – about 12 km of standard rock at 100 PeV, 20 km at 1 EeV, 23 km at 10 EeV – as the logarithmic growth of the tau range is tempered by beta rising.

Parameters:

energy_pevfloat

Neutrino energy, in PeV.

beta_cm2gfloat, optional

Tau energy-loss coefficient, in cm^2/g.

density_gcm3float, optional

Rock density, in g/cm^3.

inelasticityfloat, optional

Mean charged-current inelasticity.

samplesint, optional

Points on the log-spaced depth grid searched for the peak.

Returns:

float

Column depth maximising the exit probability, in g/cm^2.

Examples

>>> from oroscope import physics
>>> km = physics.production_escape_optimum_gcm2(100.0) / physics.CRUST_DENSITY_GCM3 / 1e5
>>> 8.0 < km < 16.0
True

oroscope.physics.depth_band_from_energy(energy_min_pev: float, energy_max_pev: float, fraction: float = 0.5, beta_cm2g: float | None = None, density_gcm3: float = 2.65, inelasticity: float = 0.2, samples: int = 400) → tuple[float, float]

Column-depth band where the tau exit probability stays above fraction of peak.

Very wide – roughly 5e5 to 2.6e8 g/cm^2 at half maximum across 100 PeV to 10 EeV, some two and a half decades – which is itself the result: column depth is an intrinsically weak discriminant and the criterion should not pretend otherwise.

Parameters:

energy_min_pev, energy_max_pevfloat

Ends of the neutrino energy range, in PeV.

fractionfloat, optional

Fraction of peak exit probability defining the band edges.

beta_cm2gfloat, optional

Tau energy-loss coefficient, in cm^2/g.

density_gcm3float, optional

Rock density, in g/cm^3.

inelasticityfloat, optional

Mean charged-current inelasticity.

samplesint, optional

Points on the log-spaced depth grid.

Returns:

tuple of float

(low_gcm2, high_gcm2), the low edge taken at the lowest energy and the high edge at the highest, so the band spans the whole requested range.

Examples

>>> from oroscope import physics
>>> lo, hi = physics.depth_band_from_energy(100.0, 10000.0)
>>> lo < hi
True

oroscope.physics.earth_absorption_cutoff_deg(energy_pev: float, fraction: float = 0.5, radius_m: float = 6371000.0, density_gcm3: float = 2.65, **kw: float) → float | None

Elevation below which the Earth chord itself exceeds the useful column depth.

The chord is 2R sin(theta), hundreds of km for degrees below the horizontal, so steep arrival directions carry far more matter than the optimum wants and the neutrino is absorbed before reaching the exit region. Setting the chord equal to the upper band edge gives the elevation at which acceptance has fallen to fraction of peak.

The result narrows sharply with energy – about -4.5 degrees at 100 PeV, -2.1 at 1 EeV, -1.0 at 10 EeV – so the effective arrival window is not a fixed +/-3 degrees but an energy-dependent one whose lower edge climbs toward the horizon.

Parameters:

energy_pevfloat

Neutrino energy, in PeV.

fractionfloat, optional

Fraction of peak acceptance defining the cut.

radius_mfloat, optional

Earth radius, in metres.

density_gcm3float, optional

Crust density, in g/cm^3.

**kw

Passed through to tau_exit_probability().

Returns:

float or None

A negative elevation angle in degrees, or None when the cut lies below the horizon entirely.

Examples

>>> from oroscope import physics
>>> cut = physics.earth_absorption_cutoff_deg(1000.0)
>>> -6.0 < cut < 0.0
True

oroscope.physics.spectrum_weighted_decay_probability(distance_m, energy_min_pev, energy_max_pev, spectral_index=2.0, shower_development_m=0.0, samples=96, index_samples=129, weight_by='flux', response=None)

Probability the tau decays inside the usable gap, folded over a power-law spectrum.

A tau of energy \(E\) decays within a usable path \(u\) with probability \(1 - e^{-u/L(E)}\), and \(L = (E/m_\tau)c\tau\) runs over three decades across a single experiment’s reach. Evaluating that at one representative energy is therefore not an approximation but a choice of answer: measured on a real canyon search, the reported capacity ran from 10878 detector positions at 3 PeV to zero at 100 PeV. Weighting by the flux instead gives a number that is a property of the terrain and the spectrum rather than of the energy someone picked.

\[P(u) = \frac{\int E^{-\gamma}\left(1 - e^{-u/L(E)}\right)\,{\rm d}E} {\int E^{-\gamma}\,{\rm d}E}\]

integrated on a log-spaced grid, since the range spans decades.

What weights the average is selectable. An event rate is \(\int \Phi(E) A(E) P(E)\,{\rm d}E\), and weighting by the flux alone is only one of three defensible choices:

"flux"

\(w(E) = E^{-\gamma}\). The default, and what every published number here was computed with. Says: of the neutrinos that arrive, what fraction decays usefully?

"acceptance"

\(w(E) = A(E)\), per unit energy. Says: over the energies this detector actually responds to, what fraction decays usefully?

This does not remove the spectral assumption; it replaces it. Weighting by \(A(E)\,{\rm d}E\) is weighting by \(A(E)\) against a flat differential flux, which is \(\gamma = 0\) – harder than any astrophysical spectrum, so the average leans on the top of the range rather than sitting free of it. Measured against TAMBO’s rising \(A(E)\) over 3 PeV - 1 EeV, 97.4% of the weight lands above 100 PeV, against 31.3% for flux. That is a large part of why weighting by acceptance alone emptied a search (roadmap 6.42), alongside the shape of the inferred response itself. Read it as a different assumption, useful for asking what the instrument’s own band favours, and not as the assumption-free option.

"flux_times_acceptance"

\(w(E) = E^{-\gamma} A(E)\), the event-rate integrand itself, and the right one when both the spectrum and a response table are trusted.

A(E) is not something this module can supply; pass a callable, most usefully one recovered from a published curve by aperture.infer_response(), which divides that curve by the geometric model and leaves everything else.

Note a steep spectrum weights low energies heavily, where the tau decays readily – so a soft spectrum drives \(P\) toward 1 and the term stops discriminating. That is the physics rather than a defect, but it means the spectral index deserves the same scrutiny as any other assumption – including the \(\gamma = 0\) that "acceptance" quietly substitutes for it.

Parameters:

distance_mfloat or array_like

Distance from the exit point to the detector, in metres.

energy_min_pev, energy_max_pevfloat

Ends of the tau energy range, in PeV.

spectral_indexfloat, optional

\(\gamma\) in \({\rm d}N/{\rm d}E \propto E^{-\gamma}\).

shower_development_mfloat, optional

Path the shower needs after the decay, in metres. Subtracted from the gap, so a target closer than this yields nothing usable.

samplesint, optional

Points on the log-spaced energy grid.

index_samplesint, optional

Points across the spectral-index range, when one is given. Generous by default because the index integral is folded into a weight vector computed once, so a fine grid costs nothing per candidate. Ignored for a single index.

weight_by{‘flux’, ‘acceptance’, ‘flux_times_acceptance’}, optional

What weights the average. See the discussion above. 'acceptance' ignores spectral_index entirely, which is the point of it.

responsecallable, optional

A(E) in PeV, required by the two acceptance weightings. A response that is zero everywhere on the grid leaves nothing to average and raises.

Returns:

ndarray

Flux-weighted decay probability, in [0, 1], matching the shape of distance_m.

Raises:

ValueError

If the energy range is inverted, or spectral_index is neither one value nor a pair.

See also

arrival_scan.decay_probability

the single-energy form, for one baseline window.

Examples

>>> from oroscope import physics
>>> p = physics.spectrum_weighted_decay_probability(3000.0, 3.0, 1000.0)
>>> round(float(p), 3)
0.954

A harder spectrum puts more weight at high energy, where the tau outruns the gap:

>>> soft = physics.spectrum_weighted_decay_probability(3000.0, 3.0, 1000.0, 2.7)
>>> hard = physics.spectrum_weighted_decay_probability(3000.0, 3.0, 1000.0, 1.5)
>>> bool(hard < soft)
True

Or marginalise over the index rather than choosing one, which lands between the extremes it spans:

>>> spread = physics.spectrum_weighted_decay_probability(
...     3000.0, 3.0, 1000.0, (1.5, 2.7))
>>> bool(hard < spread < soft)
True

oroscope.physics.geomagnetic_latitude_deg(latitude_deg: float, longitude_deg: float, pole_lat_deg=80.7, pole_lon_deg=-72.7)

Latitude in the centered-dipole frame.

Parameters:

latitude_deg, longitude_degfloat

Geographic coordinates of the site, in degrees.

pole_lat_deg, pole_lon_degfloat, optional

Geographic coordinates of the north geomagnetic pole, in degrees.

Returns:

float

Magnetic latitude, in degrees.

Examples

>>> from oroscope import physics
>>> round(physics.geomagnetic_latitude_deg(-16.4, -71.5), 1)
-7.1

oroscope.physics.centered_dipole_inclination(latitude_deg: float, longitude_deg: float, **kw: float) → float

Inclination from a centered dipole, tan(I) = 2 tan(magnetic latitude).

Useful when no IGRF lookup is available: it captures the dominant behaviour, that inclination passes through zero at the magnetic equator and steepens away from it. At Arequipa it gives about -14 degrees, the site being roughly 7 degrees south of the magnetic equator.

Note this approximation is only worth using for inclination. The dipole declination at Arequipa is about -0.2 degrees against an IGRF value of -6.9, so the non-dipole terms dominate there and a dipole declination would be misleading.

Parameters:

latitude_deg, longitude_degfloat

Geographic coordinates of the site, in degrees.

**kw

Passed through to geomagnetic_latitude_deg().

Returns:

float

Inclination in degrees, positive downward.

Examples

>>> from oroscope import physics
>>> round(physics.centered_dipole_inclination(-16.4, -71.5), 1)
-14.0

oroscope.physics.default_field_for_site(latitude_deg: float, longitude_deg: float, declination_deg: float | None = None, inclination_deg: float | None = None) → tuple[float, float]

Geomagnetic field for an arbitrary site, falling back sensibly.

Inclination is computed from the site’s own coordinates with the dipole model, so moving the search to another location gets the right inclination automatically – it is the quantity that varies most across Peru, from about -5 degrees near Lima to -14 near Arequipa.

Declination cannot be had that way: the dipole gives about -0.2 degrees at Arequipa against an IGRF -6.9, so non-dipole terms dominate. It therefore falls back to the Arequipa IGRF value, which is right for the prototype region and approximately right for the rest of southern Peru. Supply the IGRF declination for anywhere else.

Parameters:

latitude_deg, longitude_degfloat

Geographic coordinates of the site, in degrees.

declination_degfloat, optional

Declination in degrees east of north. Supply the IGRF value for the site; the fallback is Arequipa’s.

inclination_degfloat, optional

Inclination in degrees, positive downward. Derived from the site’s coordinates when omitted.

Returns:

tuple of float

(declination_deg, inclination_deg).

Examples

>>> from oroscope import physics
>>> dec, inc = physics.default_field_for_site(-16.4, -71.5)
>>> f"{dec:.1f} {inc:.1f}"
'-6.9 -14.0'

oroscope.physics.set_declination_model(model)

Supplies a declination model, so declination can follow the site as inclination does.

Inclination is computed from the site’s coordinates with a centred dipole and is good enough. Declination is not: the dipole gives about -0.2 degrees at Arequipa against an IGRF -6.9, because the non-dipole terms dominate it. So declination has always fallen back to a single constant – Arequipa’s IGRF value – wherever the DEM happened to be, which is right for the prototype region and wrong anywhere else.

No IGRF model is shipped, deliberately. IGRF is a spherical-harmonic expansion with a couple of hundred coefficients per epoch, and coefficients typed from memory would produce declinations that look entirely plausible and are wrong – the exact failure this project keeps finding. What is provided instead is the socket: pass any callable model(latitude_deg, longitude_deg) -> declination_deg, from ppigrf/pyIGRF, from a NOAA grid via declination_from_grid(), or from the collaboration’s own numbers.

Parameters:

modelcallable or None

model(latitude_deg, longitude_deg) returning degrees east of north. None restores the constant fallback.

Returns:

callable or None

The model now in force.

Examples

>>> from oroscope import physics
>>> _ = physics.set_declination_model(lambda lat, lon: -6.0 + 0.1 * (lon + 71.5))
>>> dec, inc = physics.default_field_for_site(-16.4, -71.5)
>>> f"{dec:.2f}"
'-6.00'

It follows the site, which the constant never did:

>>> dec_east, _ = physics.default_field_for_site(-16.4, -66.5)
>>> f"{dec_east:.2f}"
'-5.50'
>>> _ = physics.set_declination_model(None)
>>> physics.default_field_for_site(-16.4, -66.5)[0]
-6.9

oroscope.physics.declination_model() → object

The declination model in force, or None when the constant fallback is in use.

Examples

>>> from oroscope import physics
>>> physics.declination_model() is None
True

oroscope.physics.declination_from_grid(latitudes, longitudes, declinations)

Builds a declination model by bilinear interpolation over a supplied grid.

The practical way to get real IGRF values in without a spherical-harmonic implementation: export a declination grid covering the DEM (NOAA’s geomagnetic calculator will produce one), and hand the three arrays here.

Parameters:

latitudesarray_like

Grid latitudes, ascending, in degrees.

longitudesarray_like

Grid longitudes, ascending, in degrees.

declinationsarray_like

(len(latitudes), len(longitudes)) of declination in degrees east of north.

Returns:

callable

model(latitude_deg, longitude_deg) -> float, suitable for set_declination_model(). Queries outside the grid are clamped to its edge, which is the right behaviour for a DEM that pokes slightly past its corners and a poor one for a grid that does not cover the region at all – so cover it.

Examples

>>> import numpy as np
>>> from oroscope import physics
>>> lats, lons = np.array([-18.0, -14.0]), np.array([-74.0, -70.0])
>>> dec = np.array([[-5.0, -7.0], [-6.0, -8.0]])
>>> model = physics.declination_from_grid(lats, lons, dec)
>>> f"{model(-18.0, -74.0):.2f}"
'-5.00'

Halfway in both directions is the average of the four corners:

>>> f"{model(-16.0, -72.0):.2f}"
'-6.50'

oroscope.physics.geomagnetic_unit_vector(declination_deg: float, inclination_deg: float) → tuple[float, float, float]

Unit geomagnetic field in local East-North-Up coordinates.

Declination is measured from geographic north, positive eastward; inclination is positive downward, so it enters the Up component with a minus sign.

Parameters:

declination_degfloat

Declination, in degrees east of north.

inclination_degfloat

Inclination, in degrees, positive downward.

Returns:

tuple of float

(east, north, up) components of the unit field vector.

Examples

>>> from oroscope import physics
>>> e, n, u = physics.geomagnetic_unit_vector(0.0, 0.0)
>>> f"{e:.1f} {n:.1f} {u:.1f}"                  # due north, horizontal
'0.0 1.0 -0.0'

oroscope.physics.geomagnetic_sin_alpha(azimuth_deg: float, elevation_deg: float, field_unit_vector: tuple[float, float, float]) → float

sin(alpha) between a shower axis and the geomagnetic field.

Radio emission from an air shower is dominantly geomagnetic, with amplitude proportional to |v x B|, so a shower travelling along the field radiates very little of it. Peru lies near the magnetic equator, where the field is close to horizontal and roughly northward: north-south showers are therefore strongly suppressed and east-west ones near maximal. The azimuth of a target matters, not merely whether one exists.

The sign of the axis is irrelevant – |(-v) x B| = |v x B| – so the same value applies whether the direction is taken as arrival or propagation.

Parameters:

azimuth_degfloat

Shower azimuth, in degrees clockwise from north.

elevation_degfloat

Shower elevation angle, in degrees.

field_unit_vectortuple of float

Field direction as (east, north, up), from geomagnetic_unit_vector().

Returns:

float

sin(alpha), in [0, 1]. Zero for a shower along the field.

Examples

>>> from oroscope import physics
>>> B = physics.geomagnetic_unit_vector(0.0, 0.0)      # horizontal, northward
>>> round(physics.geomagnetic_sin_alpha(0.0, 0.0, B), 3)     # along the field
0.0
>>> round(physics.geomagnetic_sin_alpha(90.0, 0.0, B), 3)    # across it
1.0

oroscope.physics.refractivity(altitude_m: float, sea_level_value: float = 0.00029, scale_height_m: float = 8400.0) → float

n - 1 at altitude, falling with density.

Parameters:

altitude_mfloat

Altitude above sea level, in metres.

sea_level_valuefloat, optional

Refractivity at sea level.

scale_height_mfloat, optional

Density scale height, in metres.

Returns:

float

Refractivity n - 1.

Examples

>>> from oroscope import physics
>>> f"{physics.refractivity(4000.0):.2e}"
'1.80e-04'

oroscope.physics.cherenkov_angle_rad(altitude_m: float, **kw: float) → float

Cherenkov angle in air, sqrt(2(n-1)) for small angles.

About 1.4 degrees at sea level and 1.1 degrees at 4000 m: the cone narrows with altitude because the air is thinner.

Parameters:

altitude_mfloat

Altitude of the emission point, in metres.

**kw

Passed through to refractivity().

Returns:

float

Cherenkov angle, in radians.

Examples

>>> import math
>>> from oroscope import physics
>>> f"{math.degrees(physics.cherenkov_angle_rad(4000.0)):.2f} deg"
'1.09 deg'

oroscope.physics.cherenkov_footprint_radius_m(altitude_m: float, distance_m: float, **kw: float) → float

Radius of the radio footprint on the ground, D * theta_C.

The consequence for layout is counter-intuitive: a higher site has a smaller footprint, so it needs a denser array for the same trigger efficiency. A 1 km grid under-samples a footprint of a few hundred metres either way, which is why counted antennas are a cost proxy rather than an effective area.

Parameters:

altitude_mfloat

Altitude of the emission point, in metres.

distance_mfloat

Distance from emission point to the ground, in metres.

**kw

Passed through to refractivity().

Returns:

float

Footprint radius, in metres.

Examples

>>> from oroscope import physics
>>> f"{physics.cherenkov_footprint_radius_m(4000.0, 20000.0):.0f} m"
'380 m'

oroscope.physics.footprint_sampling(spacing_m: float, altitude_m: float, distance_m: float, **kw: float) → float

Antennas per footprint diameter: 2 r / spacing.

Below 1 the array does not resolve the footprint and triggering relies on a single antenna happening to fall inside it.

Parameters:

spacing_mfloat

Detector spacing, in metres. Zero or less returns 0.

altitude_mfloat

Altitude of the emission point, in metres.

distance_mfloat

Distance from emission point to the ground, in metres.

**kw

Passed through to refractivity().

Returns:

float

Antennas spanning one footprint diameter.

Examples

>>> from oroscope import physics
>>> round(physics.footprint_sampling(1000.0, 4000.0, 20000.0), 2)
0.76

The arrival scan

The scan kernel: profile walking, column depth, Fresnel clearance and radio-noise line-of-sight. Compiled with Numba; see The physics for what it computes.

Arrival-direction scanning: what a candidate site can actually see, and what lies behind it.

The question a site search must answer is not “is there a tall mountain out there” but “from which arrival directions does a backward ray from this pixel enter rock, and how much rock does it cross”. This module answers that directly.

Tracing backward from a candidate along an arrival direction (azimuth phi, elevation angle theta measured from horizontal):

• rays above the local horizon escape to the sky and see no matter;

• rays below it strike terrain. That first intersection is the tau exit point, the distance to it is the decay baseline, and the path length beyond it that runs under the surface is the column depth.

Different experiments accept different arrival directions, so the elevation window is a parameter rather than a constant. GRAND accepts neutrinos within roughly +/-3 deg of the horizon and cosmic rays from above it; TAMBO looks across a canyon. The cosmic-ray case inverts the test – terrain in the accepted directions is an obstruction rather than a target – which require_terrain=False expresses without a second code path.

Algorithm

For a fixed azimuth, one walk outward yields every elevation bin at once. Writing the elevation angle of the terrain at ground distance d as

theta_terrain(d) = atan( (z(d) - d^2/2R - z0) / d )

then a ray at angle theta first meets terrain at the smallest d where theta_terrain(d) >= theta. Since the running maximum of theta_terrain only increases, each new maximum claims a contiguous band of elevation bins, so first-intersection distances for all bins are filled in a single pass.

Column depth follows from the same samples: the ray at angle theta is underground wherever theta_terrain(d) > theta, so binning theta_terrain and taking an inclusive suffix sum gives the underground path length for every bin at once. Rays crossing several ridges accumulate all of the rock they traverse, not just the first chord.

This costs one profile walk per (candidate, azimuth) regardless of how finely the elevation window is sampled.

oroscope.arrival_scan.scan(candidates, elevation, map_grid, *, elev_min_deg=-3.0, elev_max_deg=3.0, n_elev_bins=12, n_azimuths=9, half_width_deg=60.0, use_aspect=True, step_m=None, max_range_m=80000.0, min_dist_km=0.0, max_dist_km=80.0, min_depth_gcm2=0.0, require_terrain=True, min_target_slope_deg=None, max_target_slope_deg=None, rock_density=2650.0, earth_radius_m=6371000.0, radio_earth_radius_m=8500000.0, frequency_mhz=None, shower_offset_m=3000.0, antenna_height_m=2.0, near_field_m=500.0, bilinear=True, geomag_declination_deg=None, geomag_inclination_deg=None)

Convenience wrapper over scan_candidates() with defaults for GRAND neutrinos.

Sampling defaults to one DEM pixel, since a coarser step can miss a ridge entirely.

Parameters:

candidatesndarray

(N, 3) array of [row, col, aspect_deg], as produced by the topographic screen.

elevationndarray

The DEM, as a 2-D array. Converted to float32 if it is neither float32 nor float64.

map_gridMapGrid

Angular and metric pixel sizes of the DEM.

elev_min_deg, elev_max_degfloat, optional

Edges of the accepted arrival window, in degrees.

n_elev_binsint, optional

Bins across that window. Nearly free: one walk serves every bin, so cost scales with azimuths rather than with this.

n_azimuthsint, optional

Azimuths scanned per candidate. This is what sets the cost.

half_width_degfloat, optional

Half-width of a forward arc about each candidate’s aspect. None sweeps the full circle.

use_aspectbool, optional

Treat half_width_deg offsets as relative to each candidate’s aspect rather than as absolute bearings.

step_mfloat, optional

Sampling step along the profile, in metres. Defaults to one DEM pixel, since a coarser step can miss a ridge entirely.

max_range_mfloat, optional

How far to walk, in metres.

min_dist_km, max_dist_kmfloat, optional

Accepted range to the first intersection – the decay-baseline window.

min_depth_gcm2float, optional

Column depth a direction must have to count.

require_terrainbool, optional

True selects directions striking rock (neutrino channels); False selects directions escaping to the sky (cosmic-ray channels), where terrain is an obstruction and the depth and distance criteria do not apply.

min_target_slope_deg, max_target_slope_degfloat, optional

Bounds on the struck terrain’s slope along the arrival azimuth. Unset by default, which asks only that rock is present – true almost everywhere in mountainous terrain.

rock_densityfloat, optional

Density used to turn path length into column depth, in kg/m^3.

earth_radius_mfloat, optional

True Earth radius, for the particle geometry.

radio_earth_radius_mfloat, optional

Inflated radius for the refracted radio path. Used only by the Fresnel term.

frequency_mhzfloat, optional

Radio band for the Fresnel clearance measurement. None skips that pass entirely.

shower_offset_mfloat, optional

Path the shower needs after the tau decays, in metres. The far endpoint of the Fresnel measurement.

antenna_height_mfloat, optional

Height of the receiver above ground, in metres. Without it every path scores near zero, since a ground-level receiver always has terrain in its own first Fresnel zone.

near_field_mfloat, optional

Stretch of path skipped by the Fresnel measurement, in metres.

bilinearbool, optional

Interpolate the terrain profile between pixel centres. Costs about 1.44x and removes an asymmetric half-pixel bias.

geomag_declination_deg, geomag_inclination_degfloat, optional

Field direction. Both must be given for geomagnetic weighting to apply; guessing a field would be worse than declining to weight at all.

Returns:

dict of ndarray

One entry per candidate for each of: cells, solid_angle_sr, mean_distance_m, max_depth_gcm2, mean_depth_gcm2, horizon_deg, best_clearance_ratio, geomag_solid_angle_sr, path_grammage_gcm2, earth_chord_gcm2 and target_slope_deg.

See also

scoring.score_candidates

turns these observables into a comparable score.

oroscope.arrival_scan.rfi_exposure(candidates: ndarray, elevation: ndarray, map_grid, zones_rowcol_weight, step_m: float | None = None, earth_radius_m: float = 8500000.0) → ndarray

Line-of-sight-weighted radio noise exposure for each candidate.

Parameters:

candidatesndarray

(N, 3) array of [row, col, aspect_deg].

elevationndarray

The DEM.

map_gridMapGrid

Angular and metric pixel sizes.

zones_rowcol_weightiterable

Noise sources as (row, col, weight) in pixel coordinates; weight is normally the zone’s radius or a population proxy.

step_mfloat, optional

Sampling step along the sight line, in metres. Defaults to one pixel.

earth_radius_mfloat, optional

Radius for the curvature drop. The radio one, since this is a radio path.

Returns:

ndarray

Exposure in weight per metre squared, one entry per candidate; smaller is quieter. Sources hidden behind terrain contribute nothing, which a plain distance-based exclusion zone cannot express.

oroscope.arrival_scan.earth_radius_for_k(k_factor: float) → float

Effective Earth radius for a refraction k-factor.

k = 1 is true geometry; k = 4/3 is the standard radio convention and gives the 8500 km the searcher has always used. The choice is not negligible: over an 80 km path the apparent drop is 376 m at k = 4/3 against 502 m at k = 1, a difference comparable to the Fresnel clearance itself.

Parameters:

k_factorfloat

Refraction factor. 1 is true geometry; 4/3 is the standard radio convention.

Returns:

float

Effective Earth radius, in metres.

Examples

>>> from oroscope import arrival_scan
>>> f"{arrival_scan.earth_radius_for_k(4/3) / 1e3:.0f} km"
'8495 km'

oroscope.arrival_scan.azimuth_fan(n_azimuths: int, half_width_deg: float | None = None) → ndarray

Azimuths to scan, as offsets from each candidate’s aspect.

half_width_deg restricts the fan to a forward arc, which is what a slope-mounted array sees; leave it None for a full 360 degree sweep, appropriate when the array orientation does not constrain the acceptance.

Parameters:

n_azimuthsint

Number of azimuths to scan. This is what sets the cost of a search: one profile walk per (candidate, azimuth), with the elevation binning nearly free.

half_width_degfloat, optional

Half-width of a forward arc about the aspect, in degrees. None gives a full sweep.

Returns:

ndarray

Azimuth offsets, in degrees.

Examples

Samples sit at cell centres, so each stands for the same arc and the arcs tile the fan exactly. The wedge used to be endpoint-inclusive – linspace(-hw, hw, n) – which places two samples on the fan’s edges and gives them the same weight as the interior ones. The solid angle then only came out right when every azimuth accepted: for azimuth_fan(9, 60) the true arcs are 7.5, 15x7, 7.5 degrees against a uniform 13.333, so a candidate open only at the two edges reported 1.78x the sky it saw, and one open everywhere but the edges reported 0.89x. Centring the samples makes span / n the exact arc for every one of them, which is the same convention the full sweep already used.

With an odd n_azimuths the centre sample still lands exactly on the aspect.

>>> from oroscope import arrival_scan
>>> arrival_scan.azimuth_fan(4, None)
array([  0.,  90., 180., 270.])
>>> arrival_scan.azimuth_fan(3, 60.0)
array([-40.,   0.,  40.])

oroscope.arrival_scan.balanced_order(n_candidates: int, n_threads: int, block: int = 256) → ndarray | None

Candidate ordering that balances thread load without destroying locality.

Numba’s prange schedules statically, giving each thread one contiguous slice of the index range. Candidates leave the topographic screen in spatial order, so that slice is a contiguous patch of map – and walk cost varies enormously across the map, since rays near an edge terminate early while interior ones run the full range. The result is that some threads finish long before others: measured scaling was 2.4x on 12 cores against 4-5x for randomly scattered candidates.

Shuffling fixes the balance but destroys cache locality, and measured barely better overall. Dealing blocks of neighbouring candidates round-robin keeps locality inside a block while spreading each thread’s slice across the whole map.

Parameters:

n_candidatesint

Number of candidates to be scanned.

n_threadsint

Threads the scan will run on.

blockint, optional

Candidates per block. Large enough to keep locality inside a block, small enough that dealing them spreads each thread across the map.

Returns:

ndarray or None

An index array, or None when reordering cannot help – one thread, or too few candidates for the deal to be worth its own cost.

Examples

>>> from oroscope import arrival_scan
>>> arrival_scan.balanced_order(100, 1) is None       # nothing to balance
True
>>> order = arrival_scan.balanced_order(20000, 8)
>>> sorted(order.tolist()) == list(range(20000))      # a permutation, nothing lost
True

oroscope.arrival_scan.tau_decay_length_m(energy_pev: float) → float

Lorentz-boosted tau decay length, (E/m) * c*tau.

This is what sets the scale of the useful detector-to-exit-point distance, and it is exactly analytic — no simulation input required. Worth noting that it reproduces the published numbers on both sides: 1-100 PeV gives 49 m to 4.9 km, matching TAMBO’s quoted 50 m - 5 km range, while the searcher’s inherited 10-80 km GRAND window corresponds to 0.2-1.6 EeV.

Parameters:

energy_pevfloat

Tau energy, in PeV.

Returns:

float

Decay length in the laboratory frame, in metres.

Examples

>>> from oroscope import arrival_scan
>>> f"{arrival_scan.tau_decay_length_m(1.0):.0f} m"
'49 m'

oroscope.arrival_scan.energy_pev_for_decay_length(distance_m: float) → float

Inverse of tau_decay_length_m(), for reporting what a distance implies.

Parameters:

distance_mfloat

A decay length, in metres.

Returns:

float

The tau energy having that decay length, in PeV.

Examples

>>> from oroscope import arrival_scan
>>> round(arrival_scan.energy_pev_for_decay_length(49000.0))
1000

oroscope.arrival_scan.decay_probability(min_dist_m: float, max_dist_m: float, energy_pev: float) → float

Probability the tau decays inside the accepted baseline window.

exp(-d_min/L) - exp(-d_max/L). One of the few factors that needs no acceptance table, and the one that couples most strongly to site geometry.

Parameters:

min_dist_m, max_dist_mfloat

Ends of the accepted baseline window, in metres.

energy_pevfloat

Tau energy, in PeV.

Returns:

float

Probability of decaying inside the window, in [0, 1].

Examples

>>> from oroscope import arrival_scan
>>> round(arrival_scan.decay_probability(0.0, 3000.0, 3.0), 3)
1.0
>>> round(arrival_scan.decay_probability(0.0, 3000.0, 1000.0), 3)
0.059

oroscope.arrival_scan.distance_window_from_energy(energy_min_pev: float, energy_max_pev: float, shower_development_m: float = 3000.0) → tuple[float, float]

A decay-baseline window implied by an energy range.

The tau must decay before reaching the detector and the shower then needs room to develop, so the window runs from about one decay length at the low end to one at the high end plus the shower length. Ref. [2] quotes 3-10 km of shower development.

This is a stated convention rather than a derivation: it fixes the scale correctly, but the useful window also depends on acceptance details this tool does not model. Callers can always set the distances directly.

Parameters:

energy_min_pev, energy_max_pevfloat

Ends of the tau energy range, in PeV.

shower_development_mfloat, optional

Path the shower needs after the tau decays, in metres.

Returns:

tuple of float

(min_dist_m, max_dist_m).

Examples

>>> from oroscope import arrival_scan
>>> lo, hi = arrival_scan.distance_window_from_energy(1.0, 100.0)
>>> f"{lo:.0f} m to {hi / 1000:.1f} km"
'49 m to 7.9 km'

Scoring

Score shapes — band, saturating, ramp — and their composition into a single figure of merit. Read the warning in Assumptions and limitations about thresholding a product before choosing min_score.

Scoring: turning geometric observables into comparable site quality.

Every criterion returns a value in [0, 1] with a documented shape, so criteria can be combined and so a site’s weakness can be attributed to a named component rather than disappearing into a single opaque number.

Why scores rather than cuts. Measurement on real terrain showed the binary geometric test carries almost no discriminating power: in the Andes the whole of a +/-3 degree arrival window sits below the local horizon, so nearly every direction strikes rock and a hit/no-hit criterion selects most of the map. What separates sites is how much rock, at what distance, over how much solid angle – all continuous quantities.

What these scores are and are not. They rank sites against each other for one experiment and energy band. They are not apertures. Any factor that is energy-dependent but site-independent cancels in a ranking, which is precisely why useful relative scoring is possible without the differential acceptance table that absolute apertures would need (see the roadmap, section 4.10).

oroscope.scoring.band_score(x: float | ndarray, lo: float, hi: float, soft_lo: float | None = None, soft_hi: float | None = None) → ndarray

A plateau of 1 between lo and hi, falling linearly to 0 outside it.

The shape criteria of this kind want: acceptable over a range, degrading either side, rather than a cliff at an arbitrary threshold. soft_lo/soft_hi set the width of the falling flanks and default to a quarter of the band width, so a value must be well outside the band before it scores zero.

Column depth is the motivating case: the tau must be produced, which needs rock, and must escape, which limits how much – so its score is a band with an optimum, not a floor.

Parameters:

xfloat or array_like

Value or values to score.

lo, hifloat

Edges of the plateau. Swapped automatically if given the wrong way round.

soft_lo, soft_hifloat, optional

Widths of the falling flanks below lo and above hi. Default to a quarter of the band width each. Zero gives a hard cliff.

Returns:

ndarray

Scores in [0, 1].

Examples

>>> from oroscope import scoring
>>> float(scoring.band_score(5.0, 0.0, 10.0))          # inside the plateau
1.0
>>> float(scoring.band_score(-5.0, 0.0, 10.0))         # well below it
0.0

oroscope.scoring.saturating_score(x: float | ndarray, half_value: float) → ndarray

x / (x + half_value): rises from 0, reaching 0.5 at half_value.

For quantities where more is better with diminishing returns and no natural maximum – accepted solid angle being the case in hand.

Parameters:

xfloat or array_like

Value or values to score. Negative values are clipped to zero.

half_valuefloat

Value scoring 0.5. Choosing it far below the range actually observed saturates the term so that it stops discriminating, which is a real failure mode: a GRAND-scale value flattens completely against a canyon, which sees several times the sky. Note this function takes a bare number and does not care what it measures – score_candidates() feeds it the fraction of available sky for exactly that reason, since an absolute in steradians silently encodes the azimuth fan and the arrival window.

Returns:

ndarray

Scores in [0, 1).

Examples

>>> from oroscope import scoring
>>> float(scoring.saturating_score(0.05, 0.05))
0.5
>>> round(float(scoring.saturating_score(0.7, 0.05)), 3)   # saturated
0.933

oroscope.scoring.ramp_score(x: float | ndarray, zero_at: float, one_at: float) → ndarray

Linear ramp from 0 at zero_at to 1 at one_at; either order.

Parameters:

xfloat or array_like

Value or values to score.

zero_atfloat

Value scoring 0.

one_atfloat

Value scoring 1. May be below zero_at, giving a falling ramp.

Returns:

ndarray

Scores in [0, 1].

Examples

>>> from oroscope import scoring
>>> [float(scoring.ramp_score(v, 0.0, 10.0)) for v in (-1.0, 5.0, 99.0)]
[0.0, 0.5, 1.0]

oroscope.scoring.compose(components: dict[str, ndarray], mode: str = 'product', weights: dict[str, float] | None = None) → ndarray

Combines named component scores into one value in [0, 1].

Parameters:

componentsdict

Component name to array of scores in [0, 1].

mode{‘product’, ‘mean’, ‘min’}, optional

How to combine them. product is unforgiving – one bad component sinks the site – mean lets a strong component compensate, and min reports the weakest link.

weightsdict, optional

Per-component weights, used by product as exponents and by mean as linear weights. A weight of 0 excludes a component, in every mode. min has no meaningful notion of a relative weight – the smallest component is the smallest however it is scaled – so 0 is the only value it reads.

Returns:

ndarray

The composed score, in [0, 1].

Raises:

ValueError

If no components are given, the mode is unknown, a weight is negative, the weights for mean do not sum to a positive value, or every component has been excluded by a zero weight.

Warns:

UserWarning

When a weight names a component this composition does not contain. The name is a real one – a typo is rejected earlier, by parse_score_weights() – so this means the component is switched off in this run and the weight does nothing.

Notes

A product of several components each in [0, 1] concentrates near zero, so a threshold on the result sits on a cliff. Measured on a real search, a score cut of 0.0, 0.35 and 0.5 gave 65268, 10437 and 0 detector positions. Prefer ranking sites over thresholding a product.

A weight naming a component that is not here used to be dropped in silence, by an if n in w filter with nothing behind it. Combined with an unvalidated spelling at the command line, that made --score_weights geomag=0 – for geomagnetic – a request the tool accepted, ignored, and never mentioned again: the component the user had switched off ran at full weight and every number in the run moved.

Examples

>>> import numpy as np
>>> from oroscope import scoring
>>> parts = {"a": np.array([0.5]), "b": np.array([0.5])}
>>> float(scoring.compose(parts, "product")[0])
0.25
>>> float(scoring.compose(parts, "mean")[0])
0.5

A zero weight excludes a component under min as it does under product:

>>> float(scoring.compose(parts, "min", {"a": 0.0})[0])
0.5

oroscope.scoring.score_candidates(observables: dict[str, ndarray], config: dict | None = None, distance_window_m: tuple[float, float] | None = None) → tuple[ndarray, dict[str, ndarray]]

Scores candidates from their scan observables.

Components:

• depth band on column depth: enough rock to interact, not so much that

  the tau cannot escape;

• distance band on the exit-point distance, defaulting to the configured

  decay-baseline window;

• solid_angle saturating in accepted solid angle, since more acceptance is

  better with diminishing returns;

• decay probability the tau decays in the gap with room left for a shower.

  Present when either an energy range (decay_energy_min_pev and decay_energy_max_pev, folded over the spectrum – preferred) or a single decay_energy_pev is supplied;

• clearance present only when a Fresnel frequency was configured.

Parameters:

observablesdict

Per-candidate arrays from arrival_scan.scan(). cells, solid_angle_sr, mean_distance_m and max_depth_gcm2 are required; the rest enable their components when present.

configdict, optional

Overrides for DEFAULT_SCORE_CONFIG. None values are ignored, so a sparse configuration leaves the remaining defaults alone.

distance_window_mtuple of float, optional

Fallback distance band, used when the configuration does not set one. Normally the decay-baseline window the scan was run with.

Returns:

totalndarray

The composed score per candidate. Zero wherever no direction was accepted, whatever the components say.

componentsdict

The individual component scores, by name, so a site’s weakness can be attributed rather than disappearing into one number.

Examples

>>> import numpy as np
>>> from oroscope import scoring
>>> obs = {"cells": np.array([4, 0]),
...        "solid_angle_sr": np.array([0.4, 0.4]),
...        "mean_distance_m": np.array([2.0e4, 2.0e4]),
...        "max_depth_gcm2": np.array([1.0e6, 1.0e6])}
>>> total, parts = scoring.score_candidates(obs)
>>> float(total[1])                      # no accepted direction scores zero
0.0
>>> sorted(parts)
['depth', 'solid_angle']

oroscope.scoring.summarize_scores(scores: ndarray, components: dict[str, ndarray]) → dict

Distribution summary of a score set, for reporting and per-site records.

Storing a distribution rather than a single number is deliberate: it lets a site’s quality be re-examined later without re-running the terrain analysis.

Parameters:

scoresndarray

Composed scores, one per candidate.

componentsdict

The per-component scores behind them, from score_candidates().

Returns:

dict

Count, mean, median and 90th percentile of the total, and the mean of each named component so a weak site can be attributed rather than merely ranked.

Aperture

Aperture estimate, tabulated response, and inference of a response curve from a published one.

Aperture estimation, and validation against what is independently known.

The tool’s own geometry gives a site’s usable area and the solid angle over which it sees rock. Turning that into an aperture in m^2 sr needs a detection response, and the part of that response which depends on tau production and escape is exactly what no available table supplies (see the roadmap, section 4.10). This module therefore separates cleanly into:

• the geometric aperture, computed here and fully determined by terrain,

• the analytic decay factor, also computed here and free of free parameters,

• a pluggable response, defaulting to unity and replaceable by a table.

Under that split the absolute normalisation is unknown but the shape in energy and the ranking between sites are not. Both are testable, and are tested.

On validating against published apertures. Ref. [1] Fig. 25 and ref. [2] Fig. 3 are integral quantities over a whole array, all geometries and one site, so they cannot be applied per pixel. They can anchor the normalisation once supplied as data. Reading numbers off a published figure by eye is not a measurement, so this module provides the machinery to compare against a supplied curve rather than a transcription of one. What is validated here are the physical invariants the estimate must satisfy regardless of normalisation.

oroscope.aperture.unit_response(energy_pev: float | ndarray) → ndarray

Default detection response: energy-independent and equal to one.

A placeholder that makes the normalisation explicit rather than hidden. Replace it with a TabulatedResponse once a real acceptance table is available.

Parameters:

energy_pevfloat or array_like

Energies at which to evaluate the response, in PeV.

Returns:

ndarray

Ones, matching the shape of the input.

Examples

>>> from oroscope import aperture
>>> aperture.unit_response([1.0, 10.0]).tolist()
[1.0, 1.0]

class oroscope.aperture.TabulatedResponse(energy_pev, response)

Detection response interpolated from a supplied table.

Log-log interpolation, since both response and energy span decades. Supply a two-column CSV of energy in PeV against relative response, or the arrays directly.

Outside the tabulated range the response is returned as zero rather than extrapolated: beyond the table the response is unknown, and a silent extrapolation over decades of energy would be an invention.

Parameters:

energy_pevarray_like

Tabulated energies, in PeV. Sorted internally, so any order will do.

responsearray_like

Relative response at each energy, same shape.

Raises:

ValueError

If the two are not matching one-dimensional arrays.

Examples

>>> from oroscope import aperture
>>> r = aperture.TabulatedResponse([1.0, 10.0, 100.0], [0.1, 1.0, 0.5])
>>> round(float(r(10.0)), 3)
1.0
>>> float(r(1000.0))                 # outside the table: unknown, not extrapolated
0.0

classmethod from_csv(path: str) → TabulatedResponse

Builds a response from a two-column CSV of energy in PeV against response.

Parameters:

pathstr

Path to the CSV. # introduces a comment.

Returns:

TabulatedResponse

oroscope.aperture.geometric_aperture_m2sr(area_km2: float | ndarray, solid_angle_sr: float | ndarray) → ndarray

Area times accepted solid angle, in m^2 sr.

The purely geometric part of an aperture: no physics beyond the terrain.

Parameters:

area_km2float or array_like

Usable area, in km^2.

solid_angle_srfloat or array_like

Accepted solid angle, in steradians.

Returns:

ndarray

Geometric aperture, in m^2 sr.

Examples

>>> from oroscope import aperture
>>> f"{float(aperture.geometric_aperture_m2sr(100.0, 0.5)):.2e}"
'5.00e+07'

oroscope.aperture.aperture_vs_energy(area_km2: float, solid_angle_sr: float, min_dist_m: float, max_dist_m: float, energies_pev, response=None)

Aperture as a function of energy for one site.

A(E) = area * Omega * P_decay(E) * response(E)

P_decay is the probability the tau decays inside the accepted baseline window, exp(-d_min/L) - exp(-d_max/L) with L = (E/m_tau) c*tau. It is exact and carries the whole geometric energy dependence: short baselines favour low energies, long baselines high ones.

Parameters:

area_km2float

Usable area of the site, in km^2.

solid_angle_srfloat

Accepted solid angle, in steradians.

min_dist_m, max_dist_mfloat

Ends of the accepted decay-baseline window, in metres.

energies_pevarray_like

Energies at which to evaluate, in PeV.

responsecallable, optional

Detection response as a function of energy. Defaults to unit_response(), which leaves the normalisation explicit.

Returns:

ndarray

Aperture in m^2 sr, one entry per energy.

Examples

>>> import numpy as np
>>> from oroscope import aperture
>>> a = aperture.aperture_vs_energy(100.0, 0.5, 1.0e4, 8.0e4, [1.0, 100.0])
>>> bool(a[1] > a[0])          # a long baseline favours higher energies
True

oroscope.aperture.peak_energy_pev(min_dist_m: float, max_dist_m: float, energies_pev: ndarray | None = None) → float

Energy at which the decay factor peaks for a given baseline window.

A site’s geometry therefore predicts which energies it is best suited to, which is a check that can be made without any normalisation.

Parameters:

min_dist_m, max_dist_mfloat

Ends of the accepted decay-baseline window, in metres.

energies_pevndarray, optional

Energies to search, in PeV. Defaults to a wide log-spaced grid.

Returns:

float

Energy at which the decay factor peaks, in PeV.

Examples

>>> from oroscope import aperture
>>> near = aperture.peak_energy_pev(1.0e3, 5.0e3)
>>> far = aperture.peak_energy_pev(1.0e4, 8.0e4)
>>> bool(far > near)           # longer baselines peak at higher energy
True

oroscope.aperture.infer_response(published_energy_pev: ndarray, published_value: ndarray, area_km2: float, solid_angle_sr: float, min_dist_m: float, max_dist_m: float, min_model_fraction: float = 0.001) → tuple[ndarray, ndarray]

Response function implied by a published curve, given our geometric model.

This is the useful thing to do with an integral aperture when no differential acceptance table exists. Our model supplies the two factors terrain and kinematics determine – geometric aperture and the analytic tau decay probability – so dividing a published curve by them leaves everything else:

response(E) = published(E) / (area * Omega * P_decay(E))

What remains is the neutrino interaction and tau exit probability, the trigger efficiency, and any normalisation the published configuration carries. Its shape in energy is then usable as a weight for a site of the same experiment, which is exactly the piece section 4.10 had to leave out.

The caveat is the same one that applies to the published curves themselves: they are integral over one array and one site, so the inferred response inherits that site’s geometry. It is a better weight than a flat response, not a substitute for a differential table.

Where the decay probability is negligible the division is ill-conditioned and the ratio explodes – at 0.35 PeV over a canyon baseline it is of order 1e-8 – so energies whose model value falls below min_model_fraction of its own peak are excluded rather than allowed to dominate the normalisation.

Parameters:

published_energy_pevarray_like

Energies of the published curve, in PeV.

published_valuearray_like

Published aperture at each energy, in m^2 sr.

area_km2float

Usable area of the configuration the curve describes, in km^2.

solid_angle_srfloat

Accepted solid angle of that configuration, in steradians.

min_dist_m, max_dist_mfloat

Ends of its accepted decay-baseline window, in metres.

min_model_fractionfloat, optional

Energies whose model aperture falls below this fraction of its own peak are excluded, since the division there is ill-conditioned.

Returns:

energies_pevndarray

The well-conditioned subset of the input energies.

responsendarray

Inferred response over that range, normalised to 1 at its maximum.

oroscope.aperture.load_curve_csv(path: str) → tuple[ndarray, ndarray]

Loads a two-column digitized curve, returning (energy_pev, value).

The files under data/ store energy in GeV, as the published axes do, and this converts to PeV on the way in.

Parameters:

pathstr

Path to a two-column CSV of energy in GeV against value.

Returns:

energy_pevndarray

Energies, converted to PeV.

valuendarray

The second column, unchanged.

oroscope.aperture.summarize_sites(site_details: list[dict], min_dist_m: float, max_dist_m: float, energies_pev: ndarray, response: ndarray | None = None) → dict

Aperture-versus-energy for every site that carries scan observables, plus the total.

Uses each site’s median accepted solid angle, so a site is credited with the acceptance typical of its pixels rather than of its best one.

Parameters:

site_detailslist of dict

Site records from the pipeline. Sites without an arrival_scan entry are skipped.

min_dist_m, max_dist_mfloat

Ends of the accepted decay-baseline window, in metres.

energies_pevarray_like

Energy grid, in PeV.

responsecallable, optional

Detection response. Defaults to unit_response().

Returns:

dict

energies_pev, a sites mapping of site id to aperture curve, and total, their sum.

oroscope.aperture.array_scale_factor(target_units, published, target_spacing_km=None, target_grid_type=None)

How much larger this array is than the one a published curve was simulated for.

A published aperture or effective area belongs to a specific array at a specific site. Oroscope changes both, and only one of them can be corrected for by arithmetic. This is that one: the array size.

Aperture scales with instrumented ground, not with detector count as such, so the factor is

(N_target * s_target^2) / (N_published * s_published^2)

which collapses to the ratio of counts when the two spacings agree. Both terms are carried because they must be: doubling the count at fixed spacing doubles the ground and roughly doubles the aperture, while doubling it at fixed ground only makes the array denser, and a denser array past the point where it already samples the Cherenkov footprint adds very little. Scaling a densified array by its count alone would inflate the answer by exactly the factor by which it was densified.

Parameters:

target_unitsint

Detectors oroscope fits on this ground.

publisheddict or str

An entry of PUBLISHED_ARRAYS, or its key.

target_spacing_kmfloat, optional

Lattice spacing oroscope used. Defaults to the published spacing, which makes the factor a plain ratio of counts.

target_grid_type{‘hex’, ‘square’}, optional

Lattice oroscope used. Ground per detector is spacing^2 times sin60 for a triangular lattice and 1 for a square one, so the factor cancels only when both lattices match. Defaults to the published array’s.

Returns:

float

Multiplier to apply to the published curve.

Raises:

ValueError

If the published entry is unknown or either count is not positive.

Notes

This corrects the array and not the site. The published simulation carries its own terrain – Colca’s walls for TAMBO, a prototypical site for GRAND – with its own distribution of column depth, arrival elevation and target distance, and no operation on an integral curve can remove them. A scaled curve therefore reads: “what this many detectors would have achieved on the ground the simulation assumed”, not “what they will achieve on this ground”. See Assumptions and limitations.

Examples

>>> from oroscope import aperture
>>> round(aperture.array_scale_factor(10000, "tambo_aperture_fig3"), 3)
2.0

Same detector count, spread twice as far apart, is four times the ground:

>>> round(aperture.array_scale_factor(5000, "tambo_aperture_fig3",
...                                   target_spacing_km=0.30), 3)
4.0

oroscope.aperture.scale_published_curve(values, target_units, published, target_spacing_km=None, target_grid_type=None)

A published curve rescaled to the array oroscope actually found room for.

Parameters:

valuesarray_like

The published curve, in whatever units it came in – m^2 sr for an aperture, cm^2 for an effective area. The scaling is dimensionless, so the units survive.

target_unitsint

Detectors oroscope fits on this ground.

publisheddict or str

As array_scale_factor().

target_spacing_kmfloat, optional

As array_scale_factor().

target_grid_type{‘hex’, ‘square’}, optional

As array_scale_factor().

Returns:

ndarray

The curve, scaled.

Examples

Half the published detector count is half the aperture, at the same spacing:

>>> from oroscope import aperture
>>> published = [10.0, 100.0, 1000.0]
>>> aperture.scale_published_curve(published, 2500, "tambo_aperture_fig3")
array([  5.,  50., 500.])

oroscope.aperture.absolute_from_published(results, curve_path, published, target_spacing_km=None, target_grid_type=None)

A run’s absolute aperture, by scaling a published curve to the array it found room for.

This is post-processing, exactly as roadmap §4.10 step 4 intends: the search stores what terrain determines, and folding a published curve against it needs no re-run.

Parameters:

resultsdict

A results dictionary, or one loaded from a run’s JSON.

curve_pathstr

Two-column digitized curve, as in data/.

publisheddict or str

The array that curve was simulated for – an entry of PUBLISHED_ARRAYS or its key.

target_spacing_kmfloat, optional

Spacing this run used. Read from the run’s own spacing_km when omitted, which is what makes the answer honest: reading it from anywhere else invites the density error array_scale_factor() exists to prevent.

target_grid_type{‘hex’, ‘square’}, optional

Lattice this run used. Read from the run’s own grid_type when omitted.

Returns:

dict

{"energies_pev", "aperture", "units", "scale_factor", "detectors", "published", "caveat"}. aperture carries the published curve’s own units.

Notes

This corrects the array, not the site, and the returned caveat says so in the artefact itself rather than only here. The published simulation carries its own terrain – its column depths, its target distances, its trigger geometry – and no operation on an integral curve separates those from it. Read the result as what this many detectors would have achieved on the ground that simulation assumed.

Examples

Paths are resolved by the caller, so this example builds one from the package’s own location rather than from the working directory – the test job runs from tests/.

This resolves only in a source checkout. data/ sits beside src/ in the repository and is not shipped in the wheel, so from an installed site-packages/oroscope/aperture.py the three dirname calls land outside site-packages entirely. That is not a limitation of the function – it takes whatever path it is given – but the example below cannot be pasted into a fresh install unchanged. Point it at your own copy of the curve, or at a clone.

>>> import os
>>> from oroscope import aperture
>>> repo = os.path.dirname(os.path.dirname(os.path.dirname(aperture.__file__)))
>>> curve = os.path.join(repo, "data", "tambo_aperture_fig3.csv")
>>> results = {"results": {"total_capacity": 10000},
...            "parameters": {"antenna_spacing_km": 0.15}}
>>> out = aperture.absolute_from_published(results, curve, "tambo_aperture_fig3")
>>> out["units"], round(out["scale_factor"], 3), out["detectors"]
('m^2 sr', 2.0, 10000)

The pipeline

Screening, morphology, capacity and outputs, plus the command-line entry point.

The search itself: from a digital elevation model to sites with detector capacity.

Six stages, each streamed so that a DEM larger than memory is not a special case. Terrain is screened by slope, aspect, altitude and exclusion zones; the survivors are scanned over a fan of arrival directions (arrival_scan) and scored against per-experiment criteria (scoring); the accepted mask is cleaned morphologically, labelled into sites and packed with a detector lattice; and the result is written as GeoTIFF, world file, KML, PNG and JSON, with a selection funnel, a provenance record and a plain-language summary (explain).

Three things are worth knowing before reading further.

GRAND and TAMBO are configurations, not code paths. Adding an experiment means writing a JSON file. Nothing that shapes a result is hard-coded; every criterion is a parameter of find_grand_regions_interactive(), and the command line, the configuration file and the library all reach the same function.

The funnel is the diagnostic. Every filter records how many pixels survived it, so a search that returns little or nothing names the constraint responsible rather than leaving it to be guessed. It is printed, stored in the results JSON, and read back by explain.binding_constraint().

Geometry comes from the file. Pixel size and the north-west corner are read from the DEM’s own GeoTIFF tags, because an origin typed by hand does not fail when it is wrong – it silently georeferences every output to the wrong ground. A supplied origin that disagrees with the file is reported rather than honoured in silence.

The pipeline returns its results dictionary, so a caller does not have to find and re-read the file it was just handed the path to.

oroscope.site_searcher.find_grand_regions_interactive(dem_path, cell_size_deg=None, target_antennas=10000, rfi_zones=None, origin_lat=None, origin_lon=None, min_width_km=2.0, min_altitude=None, max_altitude=None, antenna_spacing_km=1.0, min_dist_km=10.0, max_dist_km=80.0, road_map_path=None, max_road_dist_km=20.0, grid_type='hex', generate_kml=False, search_mode='distributed', min_sub_array_size=500, min_aspect_deg=None, max_aspect_deg=None, min_slope_deg=3.0, max_slope_deg=25.0, region_name=None, downsample_factor=4, run_output_dir='.', output_image_format='png', tile_size=2048, resume=False, resume_dir=None, num_cores=-1, candidate_stride=5, slope_baseline_m=None, energy_min_pev=None, energy_max_pev=None, n_azimuths=9, azimuth_half_width_deg=60.0, elev_min_deg=-3.0, elev_max_deg=3.0, n_elev_bins=12, min_column_depth_gcm2=0.0, require_terrain=True, min_target_slope_deg=None, max_target_slope_deg=None, max_range_km=None, score_percentile=None, decay_weight_by='flux', decay_response_csv=None, settlements='auto', roads_geojson=None, stop_at_target=False, decay_energy_pev=None, decay_energy_min_pev=None, decay_energy_max_pev=None, decay_spectral_index=None, shower_development_m=3000.0, gap_close_km=None, fresnel_frequency_mhz=None, refraction_k=None, antenna_height_m=2.0, fresnel_near_field_m=500.0, exclude_near_field=True, depth_band_gcm2=None, score_composition='product', score_weights=None, distance_band_m=None, solid_angle_half_sr=None, solid_angle_half_fraction=None, clearance_full_at=None, min_score=0.0, geomag_declination_deg=None, geomag_inclination_deg=None, use_geomagnetic=True, grammage_mode='radio', grammage_band_gcm2=None, grammage_maturity_gcm2=None, grammage_band_fraction=None, shower_elongation_rate_gcm2=None, shower_lambda_gcm2=None, muon_shielding_km=None, bilinear_sampling=True, nu_interaction_length_gcm2=None, max_memory_gb=None, explain=True)

The main orchestrator. Now decoupled from logic, it sets up the environment, calls the pipeline helpers in sequence, and manages memory cleanup and checkpointing.

The map resolution (cell_size_deg) is read from the DEM’s own georeferencing tags unless the caller overrides it; every metric conversion downstream derives from it.

Parameters:

dem_pathstr

Path to the input elevation GeoTIFF.

cell_size_degfloat, optional

Pixel size in degrees, overriding the GeoTIFF’s own tag.

target_antennasint, optional

Capacity wanted from a single site.

rfi_zoneslist or str, optional

Exclusion zones, or the name of a bundled set.

origin_lat, origin_lonfloat, optional

Coordinates of the DEM’s north-west corner, in degrees. Read from the file’s own ModelTiepointTag when omitted, which is the recommended use; a supplied value that disagrees with the tag by more than ~100 m is reported rather than silently honoured.

min_width_kmfloat, optional

Narrowest feature to keep, in km. 0 disables pruning, which is what a strip-shaped array along a canyon wall needs.

min_altitude, max_altitudefloat, optional

Altitude bounds, in metres.

antenna_spacing_kmfloat, optional

Detector spacing, in km.

min_dist_km, max_dist_kmfloat, optional

Accepted range to the first terrain intersection, in km – the decay-baseline window. Also sets how far the profile is walked.

road_map_pathstr, optional

Aligned GeoTIFF of distance-to-road values.

max_road_dist_kmfloat, optional

Maximum allowed distance from a road, in km.

grid_type{‘square’, ‘hex’}, optional

Lattice the detectors are placed on.

generate_kmlbool, optional

Also write a Google Earth .kml.

search_mode{‘single’, ‘distributed’}, optional

Whether one site must hold the whole array.

min_sub_array_sizeint, optional

Capacity a sub-array must reach in distributed mode.

min_aspect_deg, max_aspect_degfloat, optional

Required facing directions, in degrees clockwise from north.

min_slope_deg, max_slope_degfloat, optional

Slope band the detector site must fall in, in degrees. This is the near wall, the ground the array stands on.

region_namestr, optional

Human-readable label for the outputs.

downsample_factorint, optional

Factor at which sites are labelled and areas measured. Above 1, a feature only a few pixels wide loses area it keeps detectors on.

run_output_dirstr, optional

Directory to write into.

output_image_formatstr, optional

Extension for the overview map, such as png or pdf.

tile_sizeint, optional

Side of the square tile held in RAM at once.

resumebool, optional

Reuse a previous run’s scan buffer instead of recomputing it.

resume_dirstr, optional

Directory holding that buffer.

num_coresint, optional

Threads for the scan. -1 uses all of them.

candidate_strideint, optional

Keeps every Nth screened pixel. Measured to be unbiased; see get_candidates_chunked().

slope_baseline_mfloat, optional

Ground distance over which slope is measured, in metres. Slope is scale-dependent, so this is an explicit choice rather than an accident of the DEM’s resolution.

energy_min_pev, energy_max_pevfloat, optional

Tau energy range. When given it overrides the distance window with one derived from the decay length, and in particle mode also sets the shower band.

n_azimuthsint, optional

Azimuths scanned per candidate. This is what sets the cost of a run.

azimuth_half_width_degfloat, optional

Half-width of the fan about each candidate’s aspect, in degrees.

elev_min_deg, elev_max_degfloat, optional

Edges of the accepted arrival window, in degrees.

n_elev_binsint, optional

Bins across that window. Nearly free: one walk serves them all.

min_column_depth_gcm2float, optional

Column depth a direction must have to count, in g/cm^2.

require_terrainbool, optional

True selects directions striking rock; False selects directions escaping to the sky, which is the cosmic-ray channel.

min_target_slope_deg, max_target_slope_degfloat, optional

Bounds on the struck terrain’s slope along the arrival azimuth. This is the far wall. Unset by default, which asks only that rock is present – true almost everywhere in mountainous terrain.

max_range_kmfloat, optional

How far each profile is walked, in km. Defaults to max_dist_km. Worth setting larger for a short-range search: column depth accumulates over the whole walk, so tying the two makes the reported depth a property of where the walk stopped rather than of the target’s thickness.

score_percentilefloat, optional

Keep this percentage of viable candidates, ranked by score, instead of cutting at an absolute min_score. Preferred, and for the same reason: a rank is scale-free, so it does not move when the composition or the number of components changes.

decay_weight_by{‘flux’, ‘acceptance’, ‘flux_times_acceptance’}, optional

What weights the spectrum-folded decay probability. 'flux' asks what fraction of arriving neutrinos decays usefully and is the default; 'acceptance' asks the same over the energies the detector responds to, with no assumed spectrum; 'flux_times_acceptance' is the event-rate integrand. The latter two require decay_response_csv.

decay_response_csvstr or callable, optional

Detection response A(E) for the acceptance weightings: a path to a two-column CSV of energy in PeV against relative response, or a callable. aperture.infer_response() recovers one from a published integral curve.

settlementsstr or sequence, optional

Named places to mark on the map. 'auto' uses whichever curated list has points inside the DEM; a path to a places GeoJSON from oroscope-fetch-roads marks those; 'none' marks nothing. See resolve_settlements().

roads_geojsonstr, optional

Road geometry to draw on the map, from oroscope-fetch-roads. Context only: distinct from road_map_path, which screens candidates, this changes no result.

stop_at_targetbool, optional

In distributed mode, stop selecting sites once target_antennas is reached. Sites are ranked by capacity, so this reports the best sites for the array actually wanted rather than every patch of qualifying ground.

decay_energy_pevfloat, optional

Single energy at which to score the probability the tau decays in the gap. Superseded by the range below and kept for asking what one energy would have said: measured, the answer ran from 10878 detector positions at 3 PeV to zero at 100, so one number does not stand in for a spectrum.

decay_energy_min_pev, decay_energy_max_pevfloat, optional

Tau energy range over which to fold the decay probability against the flux. The defensible form, and what makes the result a property of the terrain rather than of the energy someone picked.

decay_spectral_indexfloat, optional

Gamma in dN/dE ~ E^-gamma for that folding (default 2.0). A softer spectrum weights low energies, where the tau decays readily, so it drives the term toward 1 – an assumption deserving the same scrutiny as any other.

shower_development_mfloat, optional

Path the shower needs after the tau decays, in metres.

gap_close_kmfloat, optional

Size of the morphological closing element, in km. Defaults to antenna_spacing_km. Closing more than doubles reported area on real terrain, so it is worth setting deliberately; 0 disables it.

fresnel_frequency_mhzfloat, optional

Radio band for the Fresnel clearance measurement. None skips that pass, as a particle experiment wants.

refraction_kfloat, optional

Refraction factor for the radio path only. The particle geometry always uses the true Earth radius.

antenna_height_mfloat, optional

Receiver height above ground, in metres.

fresnel_near_field_mfloat, optional

Stretch of path the clearance measurement skips, in metres.

exclude_near_fieldbool, optional

Apply that near-field cut-off.

depth_band_gcm2tuple of float, optional

Column-depth band scoring 1, in g/cm^2.

score_composition{‘product’, ‘mean’, ‘min’}, optional

How components combine.

score_weightsdict, optional

Per-component weights.

distance_band_mtuple of float, optional

Exit-distance band scoring 1, in metres. Defaults to the decay window.

solid_angle_half_srfloat, optional

Accepted solid angle scoring 0.5, in steradians. The 0.05 default is GRAND-scale and saturates against a canyon’s much larger acceptance.

solid_angle_half_fractionfloat, optional

Fraction of the sky the azimuth fan and arrival window could accept that scores 0.5. Dimensionless, so it does not move when either of those does – which the steradian form could not manage.

clearance_full_atfloat, optional

Fresnel clearance ratio scoring 1.

min_scorefloat, optional

Score a candidate must reach. Note a product composition concentrates near zero, so any threshold in the middle sits on a cliff.

geomag_declination_deg, geomag_inclination_degfloat, optional

Field direction, in degrees. Supply the IGRF values for the site.

use_geomagneticbool, optional

Weight directions by |v x B|. Radio only; particles do not care.

grammage_mode{‘radio’, ‘particle’}, optional

Whether atmospheric depth is scored as a maturity threshold or as a band.

grammage_band_gcm2tuple of float, optional

Explicit shower band in particle mode. Setting it disables grammage_band_fraction.

grammage_maturity_gcm2float, optional

Depth at which the radio maturity ramp reaches 1, in g/cm^2.

grammage_band_fractionfloat, optional

Fraction of peak particle content still counted as a usable shower, when the band is derived from an energy range.

shower_elongation_rate_gcm2, shower_lambda_gcm2float, optional

Shower-profile parameters: how much deeper maximum sits per decade of energy, and the Gaisser-Hillas interaction length.

muon_shielding_kmfloat, optional

Rock overburden required for muon rejection, in km.

bilinear_samplingbool, optional

Interpolate the terrain profile between pixel centres. Costs about 1.44x and removes an asymmetric half-pixel bias.

nu_interaction_length_gcm2float, optional

Neutrino interaction length, enabling the Earth-chord attenuation term.

max_memory_gbfloat, optional

Ceiling on this process’s address space, in GiB. None uses 80% of what the system reports available, so a search that outgrows the machine fails with MemoryError instead of inviting the OOM killer to choose a victim; 0 disables the cap. See preflight_memory().

explainbool, optional

Print the plain-language account of the run – what was found, which constraint set the size of the answer, and which numbers are assumptions – and save it as explanation.txt beside the results. On by default. The text is also in the returned dictionary under "explanation", and can be regenerated from any results file with explain.explain_results().

Returns:

dict

The run’s results: parameters, results (sites, capacity), funnel, regions, timings_sec, aperture, provenance, explanation and the paths of the files written. The same content as the results JSON, so a caller no longer has to find and re-read the file this just wrote. A run that finds no candidate at all still returns its funnel, which is the case where the funnel matters most.

Notes

Writes GeoTIFF, world file, PNG, optional KML, a results JSON and a provenance record into run_output_dir, and prints a selection funnel. When a search returns nothing, the funnel is the first place to look: the constraint responsible is the line where the survivor count collapses.

oroscope.site_searcher.main()

Command-line entry point: parses arguments, reconciles them against the config file and the fallbacks, validates, and runs one search.

Kept as a function rather than a bare __main__ block so the console script declared in pyproject.toml has something to call, and so the argument handling can be exercised from a test without spawning a subprocess.

oroscope.site_searcher.resolve_grid_geometry(dem_path, origin_lat, cell_size_deg=None, origin_lon=None)

Determines the sampling geometry used by the whole pipeline.

Resolution priority: explicit user value > GeoTIFF ModelPixelScaleTag > 1 arc-second.

A geographic raster has pixels that are square in degrees but not in metres: a degree of longitude shrinks with the cosine of the latitude, so a 1 arc-second pixel spans roughly 30.7 m north-south but only ~29.5 m east-west at 17 degrees south. The longitude scale is evaluated at the DEM’s centre latitude so the residual error from ignoring its north-south variation is spread evenly over the map rather than accumulating towards one edge.

Parameters:

dem_pathstr

Path to the input elevation GeoTIFF.

origin_latfloat

Latitude of the DEM’s northern edge, in degrees, for the latitude-dependent east-west scaling.

cell_size_degfloat, optional

Explicit pixel size in degrees, overriding whatever the file says.

origin_lonfloat, optional

Longitude of the DEM’s western edge, in degrees. Supplied only so the result can carry center_lon; nothing about the pixel sizes depends on it. Omitted, center_lon is None.

Returns:

MapGrid

Angular pixel size, both metric pixel sizes, the centre latitude used for the longitude scaling, and where the resolution value came from – recorded so a run’s provenance says whether the resolution was detected or asserted.

Examples

>>> from oroscope import site_searcher as ss
>>> grid = ss.resolve_grid_geometry("nonexistent.tif", -15.6, cell_size_deg=1/3600)
>>> f"{grid.cell_size_y:.1f} m x {grid.cell_size_x:.1f} m"
'30.7 m x 29.8 m'

oroscope.site_searcher.read_dem_geometry(dem_path)

Reads the angular pixel size and row count of a GeoTIFF DEM from its header.

Standard geographic (EPSG:4326) DEMs such as SRTMGL1 or AW3D30 store the pixel size in degrees, which is what the georeferenced outputs (.tfw, .kml) require. Only the header is touched, so this stays cheap on multi-gigabyte files.

Parameters:

dem_pathstr

Path to the input elevation GeoTIFF.

Returns:

tuple

Pixel size in degrees, the number of rows and the number of columns. Any of them is None when the file or the tag cannot be read, which is not an error: the caller falls back to an explicit value or to 1 arc-second.

oroscope.site_searcher.read_dem_origin(dem_path)

North-west corner of a GeoTIFF, from its ModelTiepointTag.

Standard geographic DEMs carry their own corner, so asking a user to type it is asking for a mistake that nothing catches: an origin that disagrees with the file does not fail, it silently georeferences every output to the wrong ground. Reading it removes the most error-prone input the tool has.

Parameters:

dem_pathstr

Path to the input elevation GeoTIFF.

Returns:

tuple

(latitude, longitude) of the north-west corner in degrees, or (None, None) when the file or the tag cannot be read – which is not an error, since a caller may supply the origin explicitly.

oroscope.site_searcher.resolve_origin(dem_path, origin_lat=None, origin_lon=None, tolerance_deg=0.001)

Settles the DEM’s origin, preferring the file and checking anything supplied.

Two failure modes, and the second is the dangerous one. An origin nobody supplied used to be a fatal error even though the file knows it. And an origin supplied wrongly was accepted in silence, mis-georeferencing every output – the GeoTIFF, the world file, the KML and every coordinate in the results – while the search itself ran perfectly and looked right.

So the tag wins when nothing is given, and disagreement past tolerance_deg is reported loudly rather than resolved quietly. 1e-3 degrees is about 100 m, which is a few pixels: closer than that is rounding in a config file, further is a mistake.

Parameters:

dem_pathstr

Path to the input elevation GeoTIFF.

origin_lat, origin_lonfloat, optional

Origin as supplied by the user, if any.

tolerance_degfloat, optional

Disagreement beyond which the supplied value is called out.

Returns:

tuple

(latitude, longitude, source), where source describes where the value came from and is recorded in the run’s provenance.

Examples

>>> from oroscope import site_searcher as ss
>>> lat, lon, source = ss.resolve_origin("nonexistent.tif", -15.3, -72.4)
>>> (lat, lon, source)
(-15.3, -72.4, 'supplied (DEM carries no tiepoint)')

oroscope.site_searcher.build_elevation_cache(dem_path, npy_path, block_rows=2048)

Converts a DEM to the memory-mapped float32 cache, without ever holding it in RAM.

The obvious tiff.imread(path).astype(np.float32) materialises the whole DEM and then a second full copy of it — which defeats the point of the out-of-core design the rest of the pipeline is built around, and fails outright on the multi-gigabyte DEMs this tool is meant to handle. Instead the page is decoded straight into a native-dtype file, then converted a block of rows at a time. Peak memory is one block, whatever the size of the DEM.

float32 with NaN is kept rather than the DEM’s own integer dtype: NaN propagates through the gradient and comparison chain in the screening stage, so nodata is excluded without a sentinel test in every kernel. That costs twice the disk of an int16 cache and buys correctness that would otherwise have to be re-established in half a dozen places.

Parameters:

dem_pathstr

Path to the input GeoTIFF.

npy_pathstr

Path to write the float32 memory-mapped cache to.

block_rowsint, optional

Rows converted at a time. Peak memory is one block, whatever the DEM’s size.

oroscope.site_searcher.load_dem_and_init_buffers(dem_path, temp_dir, resume=False, resume_dir=None)

Step 1 Pipeline: Converts TIF to memory-mapped NPY for rapid random access and initializes the ping-pong buffers for later morphology steps. If resume is True and resume_dir is provided, it attempts to load an existing ray-tracing buffer.

Parameters:

dem_pathstr

Path to the input elevation GeoTIFF.

temp_dirstr

Directory for the working buffers.

resumebool, optional

Reuse a previous run’s scan buffer instead of recomputing it.

resume_dirstr, optional

Directory holding that buffer.

Returns:

elevationndarray

Memory-mapped DEM.

rows, colsint

Array dimensions.

path_A, path_Bstr

Paths to the boolean ping-pong buffers.

buf_andarray

Open memory map of buffer A.

is_resumingbool

True if a previous scan buffer was loaded successfully.

oroscope.site_searcher.terrain_gradients(elevation_block: ndarray, cell_size_y: float, cell_size_x: float, smooth_y: int = 0, smooth_x: int = 0) → tuple[ndarray, ndarray]

Smoothed partial derivatives of the surface, the raw material for slope and aspect.

Smoothing before differentiating gives the average gradient over the window, which is what “slope at 1 km scale” means physically. Callers must supply a block with a halo of at least max(smooth)//2 + 1 and crop the result, otherwise the window reaches past the block edge.

Kept separate from terrain_derivatives() because the screening stage wants the gradients themselves: a slope band can be tested without ever forming the angle (see slope_band_gradient_sq()), and aspect is needed only at the few pixels that survive.

Parameters:

elevation_blockndarray

Elevation tile, including a halo of at least max(smooth)//2 + 1.

cell_size_y, cell_size_xfloat

Ground size of one pixel on each axis, in metres. They differ on a geographic grid, which is why they are separate.

smooth_y, smooth_xint, optional

Smoothing window in pixels, from slope_baseline_pixels().

Returns:

tuple of ndarray

(d/dy, d/dx), in metres per metre.

oroscope.site_searcher.terrain_derivatives(elevation_block: ndarray, cell_size_y: float, cell_size_x: float, smooth_y: int = 0, smooth_x: int = 0) → tuple[ndarray, ndarray]

Slope and aspect over a stated measurement baseline.

Parameters:

elevation_blockndarray

Elevation tile, including a halo. See terrain_gradients().

cell_size_y, cell_size_xfloat

Ground size of one pixel on each axis, in metres.

smooth_y, smooth_xint, optional

Smoothing window in pixels.

Returns:

tuple of ndarray

Slope in degrees, and aspect in degrees clockwise from north.

See also

slope_band_gradient_sq

tests a slope band without forming the angle at all, which is what the screening stage uses.

oroscope.site_searcher.slope_band_gradient_sq(min_slope_deg: float | None, max_slope_deg: float | None) → tuple[float | None, float | None]

The slope band restated as bounds on the squared gradient magnitude.

slope = atan(|grad|) rises monotonically with the gradient magnitude, so

min <= atan(sqrt(g)) <= max <=> tan(min)^2 <= g <= tan(max)^2

which tests the same pixels without a sqrt or an arctan. Bounds at or beyond the vertical, and non-positive lower bounds, are returned as None meaning “unbounded”: tan is singular at 90 degrees and every real gradient satisfies them anyway.

Parameters:

min_slope_deg, max_slope_degfloat or None

Edges of the accepted slope band, in degrees.

Returns:

tuple

Lower and upper bounds on dx^2 + dy^2. Either may be None, meaning unbounded on that side.

Examples

>>> from oroscope import site_searcher as ss
>>> lo, hi = ss.slope_band_gradient_sq(3.0, 25.0)
>>> f"{lo:.4f} {hi:.4f}"
'0.0027 0.2174'
>>> ss.slope_band_gradient_sq(0.0, 90.0)      # both edges degenerate
(None, None)

oroscope.site_searcher.slope_baseline_pixels(map_grid, slope_baseline_m: float | None) → tuple[int, int]

Converts a slope measurement baseline in metres to a per-axis window in pixels.

Slope is scale-dependent: on real Andean terrain the median slope falls from ~17.8 deg measured over the DEM’s native ~61 m to ~10.8 deg over 1 km, and the fraction passing a 3-25 deg band rises from 60% to 78%. Which of those is “the” slope depends on the footprint being deployed, so the baseline is an explicit parameter rather than an accident of the DEM’s resolution.

Parameters:

map_gridMapGrid

Angular and metric pixel sizes of the DEM.

slope_baseline_mfloat or None

Ground distance over which slope is measured, in metres. None or 0 uses the DEM’s native resolution, which on 30 m data is dominated by DEM noise.

Returns:

tuple of int

Smoothing window as (rows, columns) in pixels. (0, 0) when no baseline is requested, meaning the native gradient.

oroscope.site_searcher.get_candidates_chunked(elevation, map_grid, rfi_zones, origin_lat, origin_lon, min_alt=None, max_alt=None, min_aspect_deg=None, max_aspect_deg=None, road_map_path=None, max_road_dist_km=None, min_slope_deg=3.0, max_slope_deg=25.0, tile_size=2048, candidate_stride=5, slope_baseline_m=None, funnel=None)

Step 2 Pipeline: Memory-efficient topographic screening. Iterates over the large DEM in chunks (tiles) to find pixels that meet the primary geometrical criteria (slope, aspect, altitude) and logistics constraints (RFI distance, road distance) prior to running ray-tracing.

Parameters:

• elevation (ndarray): Full DEM array (usually memory-mapped).

• map_grid (MapGrid): Angular and metric pixel sizes of the DEM.

• rfi_zones (list): List of configured exclusion zones (circles/polygons).

• origin_lat, origin_lon (float): Reference coordinates for converting km to pixels.

• min_alt, max_alt (float): Elevation restrictions.

• min_aspect_deg, max_aspect_deg (float): Required facing directions for slopes.

• road_map_path (str): Path to an aligned TIFF containing distance-to-road values.

• max_road_dist_km (float): Maximum allowed distance from a road.

• min_slope_deg, max_slope_deg (float): Required steepness limits for detector slopes.

• tile_size (int): Size of the square chunk to process in RAM at one time.

• candidate_stride (int): Keeps every Nth surviving pixel before ray-tracing. Higher values trade spatial sampling density for speed; 1 keeps every candidate.

• slope_baseline_m (float): Ground distance over which slope is measured. None uses the DEM’s native resolution, which on 30 m data is dominated by DEM noise.

• funnel (Funnel): Optional accounting object recording per-filter survivor counts.

Parameters:

elevationndarray

Full DEM, usually memory-mapped.

map_gridMapGrid

Angular and metric pixel sizes.

rfi_zoneslist or None

Exclusion zones as ('circle', lat, lon, radius_km, name) or ('poly', [(lat, lon), ...], name).

origin_lat, origin_lonfloat

North-west corner of the DEM, in degrees, for converting zones to pixels.

min_alt, max_altfloat, optional

Altitude bounds, in metres.

min_aspect_deg, max_aspect_degfloat, optional

Required facing directions, in degrees clockwise from north. Wraps through 360 when the lower bound exceeds the upper.

road_map_pathstr, optional

Aligned GeoTIFF of distance-to-road values.

max_road_dist_kmfloat, optional

Maximum allowed distance from a road, in km.

min_slope_deg, max_slope_degfloat, optional

Slope band, in degrees. Tested on the squared gradient, so neither a square root nor an arctangent is formed over the tile.

tile_sizeint, optional

Side of the square chunk processed in RAM at once.

candidate_strideint, optional

Keeps every Nth surviving pixel. Measured to be unbiased: acceptance is identical at strides 1 and 5, and the stride-corrected area matches the stride-1 truth to 0.05%.

slope_baseline_mfloat, optional

Ground distance over which slope is measured, in metres. None uses the DEM’s native resolution, which on 30 m data is dominated by DEM noise.

funnelFunnel, optional

Accounting object recording per-filter survivor counts.

Returns:

ndarray

(N, 3) array of surviving pixels as [row, col, aspect_deg], ready for run_arrival_scan().

oroscope.site_searcher.run_arrival_scan(candidates_arr, elevation, map_grid, buf_a, scan_params, score_config=None, min_score=0.0, rfi_zones_px=None, score_percentile=None, funnel=None)

Step 3 alternative: scan arrival directions instead of casting one ray per pixel.

Marks a candidate as valid when at least one accepted (azimuth, elevation) direction strikes rock within the decay-baseline window with enough column depth. See arrival_scan.py for the geometry.

Parameters:

candidates_arrndarray

(N, 3) array of [row, col, aspect_deg].

elevationndarray

The DEM.

map_gridMapGrid

Angular and metric pixel sizes.

buf_andarray

Open memory map to mark accepted pixels in.

scan_paramsdict

Keyword arguments for arrival_scan.scan().

score_configdict, optional

Overrides for scoring.DEFAULT_SCORE_CONFIG.

min_scorefloat, optional

Absolute score a candidate must reach. Used only when score_percentile is not given. The default composition is a product, whose distribution piles up near zero, so any threshold in the middle sits on a cliff.

rfi_zones_pxsequence, optional

Radio-noise sources in pixel coordinates, enabling the exposure observable.

score_percentilefloat, optional

Keep this percentage of viable candidates, by score. Rank-based and so scale-free: preferred over min_score for exactly the reason above.

funnelFunnel, optional

Records the two stages this function decides: directions accepted, the candidates the geometry accepted, and – only when a score cut is in force – how many of those the cut kept.

Returns:

n_hitsint

Number of accepted candidates, after the score cut when one applies.

observablesdict

Per-candidate arrays, including the scores and their named components, kept for per-site aggregation.

Notes

directions accepted counts the geometry alone. It previously carried the post-cut count, identical to the score row beside it, which made the score stage invisible: binding_constraint() compares each stage against the one before, so a stage that keeps exactly 100% can never be named the binding constraint however much it removed. A run emptied by min_score was therefore blamed on the arrival geometry, and the summary told the reader to widen the arrival and distance windows – advice that could not help.

Stored results written before this carry the post-cut count under directions accepted. Only runs with a cut are affected: min_score 0 accepts every viable candidate, so every GRAND run in the store already held the geometric number and is unchanged.

oroscope.site_searcher.summarize_observables_by_site(labeled, downsample_factor, candidates_arr, observables, site_ids)

Aggregates per-candidate scan observables over each labelled site.

Storing the distributions rather than a single score is deliberate: absolute apertures can then be obtained later by folding these against an acceptance table, without re-running the terrain analysis (roadmap 4.10).

Parameters:

labeledndarray

Downsampled labelled site map.

downsample_factorint

Factor relating full-resolution candidate coordinates to labeled.

candidates_arrndarray

(N, 3) array of [row, col, aspect_deg].

observablesdict

Per-candidate arrays from run_arrival_scan().

site_idssequence of int

Sites to summarise.

Returns:

dict

Site id to summary statistics – mean, median and 90th percentile of each observable – over that site’s accepted candidates. Empty when there are no accepted candidates at all.

oroscope.site_searcher.clean_shape_artifacts(path_A, path_B, rows, cols, cell_size_y, cell_size_x, antenna_spacing_km, min_width_km, tile_size, gap_close_km=None)

Step 4 Pipeline: Prunes spatial artifacts to ensure solid, block-like arrays. Applies closing to fill gaps and opening to prune unusable tendrils.

The structuring elements are sized per axis so that they cover the requested ground distance in both directions rather than only north-south.

min_width_km = 0 degenerates the opening to a 1x1 element, i.e. an identity that only carries the closed map back into path_A. That is deliberate: a “block-like array” is a GRAND assumption, and an experiment deployed along a canyon wall is a strip a few hundred metres wide and tens of kilometres long, which the opening would delete outright.

Parameters:

path_A, path_Bstr

The two ping-pong buffers. The result is left in path_A.

rows, colsint

Array dimensions.

cell_size_y, cell_size_xfloat

Ground size of one pixel on each axis, in metres.

antenna_spacing_kmfloat

Detector spacing, used as the default closing scale.

min_width_kmfloat

Narrowest feature to keep. 0 disables pruning, which is what a strip-shaped array needs.

tile_sizeint

Side of the square tile held in RAM at once.

gap_close_kmfloat, optional

Size of the closing element, in km. Defaults to antenna_spacing_km.

Returns:

tuple of int

Set-pixel counts after closing and after pruning.

oroscope.site_searcher.apply_morphology_pingpong(source_path, dest_path, shape, dtype, operation_func, structure, desc='Processing', tile_size=2048)

Applies image morphology operations (closing/opening) on a massive memory-mapped array without loading the whole array into RAM. It reads from one file and writes to another (“ping-pong”).

Parameters:

source_path, dest_pathstr

Paths to the two .npy buffers, read and written respectively.

shapetuple of int

Shape of the arrays.

dtypedtype

Element type of the destination.

operation_funccallable

Morphological operation, applied tile by tile.

structurendarray

Structuring element.

descstr, optional

Label for the progress bar.

tile_sizeint, optional

Side of the square tile held in RAM at once.

Returns:

int

Set pixels in the result, counted while writing so the funnel accounting costs nothing extra.

oroscope.site_searcher.separable_closing(chunk, structure)

Binary closing with a rectangular structuring element, done separably.

A rectangle of ones factorises into a column and a row, so dilation or erosion by (h, w) is dilation by (h, 1) followed by (1, w). That turns an O(N h w) operation into O(N (h + w)) – about 10x for the 33x33 element a 1 km antenna spacing implies – and the result is bit-identical, not an approximation.

Parameters:

chunkndarray

Boolean tile to operate on.

structurendarray

Rectangle of ones. Its two side lengths are what the operation factorises into.

Returns:

ndarray

The closed tile.

oroscope.site_searcher.separable_opening(chunk, structure)

Binary opening with a rectangular element, separably.

Prunes features narrower than the element. See separable_closing() for why the factorisation is exact rather than an approximation.

Parameters:

chunkndarray

Boolean tile to operate on.

structurendarray

Rectangle of ones.

Returns:

ndarray

The opened tile.

oroscope.site_searcher.analyze_sites_and_capacity(path_A, elevation, rows, cols, cell_size_y, cell_size_x, downsample_factor, search_mode, target_antennas, min_sub_array_size, antenna_spacing_km, grid_type, funnel=None, origin_lat=None, origin_lon=None, cell_size_deg=None, candidates_arr=None, observables=None, stop_at_target=False)

Step 5 Pipeline: Isolates unique sites and measures their capacity mathematically. Uses SciPy labeling to find continuous regions and simulates physical grid placement.

Returns: - small_final (ndarray): Downsampled binary mask of the validated sites. - labeled_viz (ndarray): Multi-integer labeled array for color coding visualizations. - site_details (list): Dictionaries containing metadata about each valid site found. - cumulative_capacity (int): Sum of all antennas fitting in valid sites. - count (int): Total number of independent valid sites found. - region_stats (dict): Region-level accounting for the funnel report.

Parameters:

path_Astr

Buffer holding the cleaned, full-resolution site mask.

elevationndarray

The DEM, for the per-site mean aspect.

rows, colsint

Full-resolution dimensions.

cell_size_y, cell_size_xfloat

Ground size of one pixel on each axis, in metres.

downsample_factorint

Factor at which labelling and area are computed. Note area is measured on the downsampled map while capacity is measured at full resolution, so a feature only a few pixels wide loses area it keeps detectors on.

search_modestr

single or distributed, deciding which capacity threshold applies.

target_antennasint

Capacity wanted from a single site.

min_sub_array_sizeint

Capacity a sub-array must reach in distributed mode.

antenna_spacing_kmfloat

Detector spacing, in km.

grid_typestr

square or hex.

funnelFunnel, optional

Accounting object recording survivor counts.

origin_lat, origin_lonfloat, optional

The DEM’s north-west corner. With cell_size_deg, each site record gains its centre coordinates and bounding box, so a reader can find the ground without opening the raster.

cell_size_degfloat, optional

Pixel size in degrees, at full resolution.

candidates_arrndarray, optional

Candidates, for folding scan observables into each site’s record.

observablesdict, optional

Their per-candidate observables.

stop_at_targetbool, optional

In distributed mode, stop selecting sites once target_antennas is reached. Sites are sorted by capacity, so this takes the best ones and reports the array actually wanted rather than every patch of qualifying ground.

Returns:

small_finalndarray

Downsampled binary mask of the validated sites.

labeled_vizndarray

Site labels for colour coding, sized from the label count so that selecting more than 255 sites does not overflow.

site_detailslist of dict

Per-site metadata for every site that cleared the thresholds, sorted by capacity, each carrying a selected flag. With stop_at_target the list is longer than the selection: cumulative_capacity, count and the exported mask cover the selected sites only, so anything totalling this list must filter on selected or it will over-report both area and site count.

cumulative_capacityint

Total capacity across the selected sites.

countint

Number of sites selected.

region_statsdict

Region-level accounting for the funnel report.

oroscope.site_searcher.count_grid_capacity(mask_chunk, cell_size_y, cell_size_x, spacing_m, grid_type_code)

Counts the detectors that fit on the validated terrain at a given ground spacing.

Detector positions are laid out in metres on the ground and only then looked up in the pixel grid. The earlier version did the reverse: it converted the spacing to an integer number of pixels and stepped the array by that stride, which truncated three separate times — int() on the row stride, on the column stride, and on the hexagonal row pitch int(spacing_r * sin60). Every truncation shortens the spacing, so detectors ended up closer together than asked for and the count came out high: +7.4% at GRAND’s 1 km, and +58% at TAMBO’s 100 m, where only about three pixels span one separation on a 30 m DEM and the hex pitch collapsed from 2.6 to 2 pixels. Placing points in continuous coordinates has no stride to truncate, so the count follows the requested geometry at any spacing.

The layout is anchored at the chunk’s own corner rather than fitted to it, so this is a capacity estimate for an arbitrarily-placed array, not the best packing achievable by sliding the grid around. That is the same convention as before.

A spacing finer than the DEM’s own pixels is permitted and yields several detectors per pixel. That is the honest continuum limit — capacity is usable area divided by area per detector — but note the terrain mask cannot resolve whether those sub-pixel positions really are usable.

Parameters:

mask_chunkndarray

2D boolean array; True marks valid terrain.

cell_size_y, cell_size_xfloat

Ground size of one pixel, in metres. They differ on a geographic grid, which is why an equal ground spacing is a different number of pixels on each axis.

spacing_mfloat

Distance between neighbouring detectors, in metres. Zero or less returns 0.

grid_type_codeint

0 for a square grid, 1 for a hexagonal (triangular) one.

Returns:

int

Detectors fitting inside the valid terrain.

Examples

>>> import numpy as np
>>> from oroscope import site_searcher as ss
>>> mask = np.ones((100, 100), dtype=bool)          # 3 km square of 30 m pixels
>>> ss.count_grid_capacity(mask, 30.0, 30.0, 1000.0, 1)
12
>>> ss.count_grid_capacity(mask, 30.0, 30.0, 0.0, 1)   # degenerate spacing
0

oroscope.site_searcher.create_world_file(tif_filename, top_left_lat, top_left_lon, cell_size_deg)

Creates an ESRI World File (.tfw) which accompanies a standard TIFF image, allowing GIS software (like QGIS or ArcGIS) to project it correctly on a map.

Parameters:

tif_filenamestr

Path to the raster the world file accompanies. The .tfw is written beside it with a matching stem.

top_left_lat, top_left_lonfloat

Coordinates of the raster’s north-west corner, in degrees.

cell_size_degfloat

Pixel size in degrees, after any downsampling.

oroscope.site_searcher.generate_kml_file(mask, elevation, filename, origin_lat, origin_lon, cell_size_deg, downsample=1)

Generates a Google Earth compatible KML file representing the valid site polygons. It extracts polygon contours from the binary mask using Matplotlib’s contour tool.

Parameters: - mask (ndarray): Binary mask indicating valid deployment sites. - filename (str): Output path for the KML file. - origin_lat, origin_lon, cell_size_deg: Used to convert array pixel indices to GPS coordinates.

Parameters:

maskndarray

Boolean site mask.

elevationndarray

The DEM, used to place the contours in height.

filenamestr

Path to write the .kml to.

origin_lat, origin_lonfloat

North-west corner of the mask, in degrees.

cell_size_degfloat

Pixel size in degrees.

downsampleint, optional

Factor by which mask is already downsampled relative to the DEM.

oroscope.site_searcher.generate_visualizations_and_outputs(dem_path, elevation, small_final, labeled_viz, site_details, count, cumulative_capacity, origin_lat, origin_lon, map_grid, downsample_factor, generate_kml, run_output_dir, output_image_format, rfi_zones, search_mode, grid_type, antenna_spacing_km, min_altitude, max_altitude, region_name, final_params, run_info=None, settlements='auto', roads_geojson=None)

Step 6 Pipeline: Formats and exports all scientific products including geo-registered TIFs, KML models, an annotated map graphic, and a serialized JSON summary of the run parameters and results to the designated unified output directory.

Parameters:

dem_pathstr

Path to the DEM, used for the output stem.

elevationndarray

The DEM, as the map background.

small_finalndarray

Downsampled binary mask of the selected sites.

labeled_vizndarray

Site labels for colour coding.

site_detailslist of dict

Per-site records.

countint

Number of selected sites.

cumulative_capacityint

Total detector capacity across them.

origin_lat, origin_lonfloat

North-west corner of the DEM, in degrees.

map_gridMapGrid

Resolved grid geometry.

downsample_factorint

Factor relating small_final to the DEM.

generate_kmlbool

Also write a Google Earth .kml.

run_output_dirstr

Directory to write into.

output_image_formatstr

Extension for the overview map, such as png or pdf.

rfi_zonessequence

Exclusion zones, drawn on the map.

search_modestr

single or distributed.

grid_typestr

square or hex.

antenna_spacing_kmfloat

Detector spacing, in km.

min_altitude, max_altitudefloat or None

Altitude bounds applied, for the annotation.

region_namestr

Human-readable region label.

final_paramsdict

Resolved parameters, serialised into the results JSON.

run_infodict, optional

Timings, funnel and provenance to record alongside the results.

settlementsstr or sequence, optional

Named places to mark on the map. See resolve_settlements().

roads_geojsonstr, optional

Road geometry to draw as context, from oroscope-fetch-roads.

Returns:

generated_fileslist of str

Absolute paths of everything written.

out_datadict

The results as serialised into the JSON. Returned as well as written so the caller does not have to find and re-read the file it was just handed the path to, which is what every caller was doing.

oroscope.site_searcher.collect_provenance(dem_path, map_grid)

Captures everything needed to reproduce a run: code version, input identity, environment and invocation. Written alongside the scientific outputs.

Parameters:

dem_pathstr

Path to the DEM, whose sha256 is recorded.

map_gridMapGrid

Resolved grid geometry, recorded along with where its resolution came from.

Returns:

dict

Git commit and dirty flag, DEM path, size and checksum, resolved grid geometry, third-party package versions, and the module-level physics state – enough to say what produced a result months later.

oroscope.site_searcher.validate_parameters(params)

Pre-flight validation checks to enforce ‘Fail Fast’ mechanisms. Verifies the existence of critical files and the physical logic of search bounds before engaging the memory-heavy processing loops.

Parameters:

paramsdict

Fully resolved parameters, after the config, fallback and command line have been reconciled.

Raises:

SystemExit

If any check fails. Every problem is collected and reported at once rather than one per run, since the expensive stages come afterwards.

oroscope.site_searcher.parse_score_weights(value)

Normalises per-component score weights from either input form.

A config file is JSON, so it can carry a mapping directly. The command line cannot, so it takes shower=2,solid_angle=1 instead. Both end up as a dict, and anything unnamed keeps weight 1.

Parameters:

valuestr, dict or None

Either name=value pairs separated by commas, or a mapping, or None.

Returns:

dict or None

Component name to weight, or None when nothing was supplied, which leaves the composition unweighted.

Raises:

SystemExit

If a pair lacks =, its value is not a number, or it names something that is not a score component.

Notes

The names are checked. They were not, and compose() then dropped anything it did not recognise with an if n in w filter, so a misspelling was accepted, ignored and never reported. --score_weights geomag=0 – one character short of geomagnetic – parsed cleanly, and the component the user meant to switch off ran at full weight through the whole search: measured on a two-component product, 0.18 where the correct spelling gives 0.9. Nothing anywhere said so. Weights are the one input whose failure leaves no trace in the output, so they are the one input worth rejecting outright.

Examples

>>> from oroscope import site_searcher as ss
>>> ss.parse_score_weights("shower=2,depth=0.5") == {"shower": 2.0, "depth": 0.5}
True
>>> ss.parse_score_weights(None) is None
True

oroscope.site_searcher.explicitly_passed(parser, argv=None)

The set of options the user actually typed, as opposed to argparse’s defaults.

argparse gives no way to distinguish --candidate_stride 5 from the default of 5, which is why the configuration merge used to prefer a config file over the command line: with no way to tell a typed flag from an untyped one, honouring the command line would have let every default silently overwrite the config.

Re-parsing with every default suppressed answers the question directly — with SUPPRESS, argparse only sets an attribute for an option that actually appeared. The defaults are restored afterwards, so the original args is untouched.

Parameters:

parserargparse.ArgumentParser

The parser to interrogate. Left exactly as it was found.

argvlist of str, optional

Arguments to parse. Defaults to sys.argv.

Returns:

set of str

Destinations of the options that actually appeared on the command line.

oroscope.site_searcher.is_point_in_poly(x, y, poly_verts)

Determines if a given 2D point lies inside a polygon using the Ray-Casting algorithm. Optimized for Numba execution.

Parameters:

x, yfloat

Coordinates of the test point.

poly_vertsndarray

(M, 2) array of polygon vertices as (x, y).

Returns:

bool

True if the point is inside the polygon.

oroscope.site_searcher.apply_poly_mask_numba(valid_rows, valid_cols, poly_verts, mask_out)

Parallelized application of the polygon ray-casting check across an array of coordinates. Used for excluding regions defined by arbitrary polygonal RFI zones.

Parameters:

valid_rows, valid_colsndarray

Row and column coordinates of the points to check.

poly_vertsndarray

(M, 2) array of polygon vertices.

mask_outndarray

Boolean array modified in place; entries inside the polygon are cleared. Only ever clears bits, so several polygons can be applied in sequence.

class oroscope.site_searcher.Funnel

Records how many pixels survive each successive stage of the search.

A search that returns nothing gives the user no clue which constraint was responsible. Every stage reports the count of pixels that passed it and all preceding stages, so the table reads as a funnel from the raw DEM down to the selected sites. Counts accumulate across tiles.

add(name, count)

Adds to a stage’s running total, creating the stage on first use.

Parameters:

namestr

Stage label, as it will appear in the funnel table.

countint

Survivors to add. Stages accumulate across tiles, so this is called once per tile per filter.

get(name, default=0)

as_dict()

render()

Formats the funnel as a table: count, share of the DEM, share of the previous stage.

class oroscope.site_searcher.MapGrid(cell_size_deg, cell_size_y, cell_size_x, center_lat, center_lon, source)

Resolved pixel geometry of a DEM.

Parameters:

cell_size_degfloat

Angular pixel size, in degrees. A geographic DEM steps by the same angle on both axes, which is why this is a single number while the metric sizes are two.

cell_size_yfloat

North-south ground size of one pixel, in metres.

cell_size_xfloat

East-west ground size of one pixel, in metres. Smaller than cell_size_y away from the equator, by the cosine of the latitude.

center_latfloat

Latitude at which cell_size_x was evaluated, in degrees.

center_lonfloat or None

Longitude of the DEM’s middle, in degrees, or None when no origin longitude was supplied. Paired with center_lat it is a point actually on the DEM, which is what a field model should be asked about; the geomagnetic field used to be resolved at this latitude and the west edge longitude.

sourcestr

Where the resolution came from – detected from the GeoTIFF, supplied explicitly, or defaulted. Recorded so a run’s provenance says which.

cell_size_deg

Alias for field number 0

cell_size_x

Alias for field number 2

cell_size_y

Alias for field number 1

center_lat

Alias for field number 3

center_lon

Alias for field number 4

source

Alias for field number 5

oroscope.site_searcher.find_results_json(run_dir)

Locates a run’s results JSON, under either the current or the legacy prefix.

Outputs used to be named grand_search_results_* whatever the experiment. The prefix is now oroscope_results_, and both are accepted so that runs made before the rename still load – a reader that could not open last week’s output would make the rename cost more than it saves.

Parameters:

run_dirstr

A run’s output directory.

Returns:

str or None

Path to the results JSON, or None if the directory holds none.

oroscope.site_searcher.default_config(preset='default')

The full set of knobs with their default values, as a plain dictionary.

Every key the pipeline understands appears here, which is the point: a template with holes in it silently falls back for whatever it omits, and the fallbacks are the least visible input the tool has.

Parameters:

presetstr, optional

"default", or "lima"/"arequipa" to fill in that region’s origin, RFI zones, name and DEM filename.

Returns:

dict

Configuration, ready to serialise or to pass to the pipeline.

Raises:

ValueError

If the preset is not one of CONFIG_PRESETS.

Examples

>>> from oroscope import site_searcher as ss
>>> cfg = ss.default_config("arequipa")
>>> cfg["rfi_zones"], cfg["min_slope_deg"], cfg["explain"]
('arequipa', 3.0, True)

oroscope.site_searcher.generate_config(path, preset='default')

Writes a configuration template to path, creating its directory if needed.

What --generate_config does, available to anyone driving the pipeline in a loop or generating a family of runs.

Parameters:

pathstr

Destination JSON file.

presetstr, optional

As default_config().

Returns:

dict

The configuration that was written.

Notes

A generated config names every key, so it also overrides every fallback. That is intended – but note the command line still wins over both, which it did not always do.

oroscope.site_searcher.load_config(path)

Reads a configuration JSON, resolving its relative paths against its own directory.

Parameters:

pathstr

Path to the file.

Returns:

dict

Its contents, or an empty dictionary when the file does not exist – which is how the command line has always treated a missing --config_path, and matching it here keeps one behaviour rather than two.

Path-valued keys (dem_path, road_map_path, resume_dir) come back absolute, resolved against the directory holding the configuration rather than the working directory. See resolve_config_paths().

Raises:

json.JSONDecodeError

If the file exists but is not valid JSON. Unlike a missing file, that is a mistake worth failing on.

oroscope.site_searcher.estimate_peak_memory_gb(rows, cols, downsample_factor=1, candidate_stride=5, survival_fraction=0.6, n_observables=12, n_scoring_arrays=24)

Rough estimate of the anonymous memory one search will need, in GiB.

Only the allocations that can exhaust RAM are counted. The DEM itself is memory-mapped and file-backed, so the kernel can evict it under pressure and it is excluded deliberately – counting it would make every large search look impossible when the streaming design exists precisely so that it is not.

This is an estimate and says so. survival_fraction in particular is the fraction of pixels passing the topographic screen, which is terrain-dependent and not known until the screen has run; 0.6 is typical of Andean terrain at a 3-25 degree band. It is meant to catch the order-of-magnitude mistake – a full DEM at downsample_factor: 1 – rather than to predict a number.

Parameters:

rows, colsint

DEM dimensions in pixels.

downsample_factorint, optional

Factor at which sites are labelled and areas measured. Scales the labelling arrays as its inverse square, and nothing else – see the note above.

candidate_strideint, optional

Keeps every Nth screened pixel. Scales the dominant term directly.

survival_fractionfloat, optional

Fraction of pixels expected to pass the topographic screen.

n_observablesint, optional

Per-candidate arrays the scan returns.

n_scoring_arraysint, optional

Further per-candidate arrays live at the peak, inside compose: the score components, the float64 copy made of each, and the temporaries.

Returns:

float

Estimated peak anonymous memory, in GiB.

Notes

The peak is in the scoring, not in the scan. This counted only the arrays arrival_scan.scan returns, which is not where the high-water mark is: by the time scoring.compose() runs, the scan’s arrays are still live, a score component has been built alongside them for each criterion, compose has clipped a float64 copy of every component, and the composition and the scoring intermediates need several more. About three times the scan’s own count is live at once, all of it n_cand long.

Under-counting that term is not academic: it advertised 2.32 GiB for the full Arequipa DEM, which then peaked at 5.68 GiB measured RSS and died against its own cap 23 minutes in. n_scoring_arrays is calibrated on that run – 15.1M candidates, 7 components – where the anonymous share of the peak implies ~36 live per-candidate arrays against the 12 this modelled.

Re-measured after the audit, the same run peaks at 6.59 GiB resident and 7.80 GiB virtual, against an estimate of 5.08. So this remains optimistic by ~1.5 GiB on the run it was calibrated against, and the calibration above is left as it was rather than quietly re-fitted: n_scoring_arrays moves the pre-flight for every region, and changing it to chase one number is how the estimate came to be sized against the cheaper of two configurations in the first place. What matters for a caller is the other column – what this function estimates is anonymous memory, and what --max-memory-gb caps is address space, 1.2 GiB larger here. Do not set one from the other.

Note also which knob moves it. downsample_factor scales only the labelling and gradient terms, because candidates are taken on the native grid; at full-DEM scale the per-candidate terms dominate, so going from 1 to 4 cuts the estimate by about 1.4x rather than the 16x the inverse-square scaling suggests. To move the dominant term, raise candidate_stride or crop the DEM.

Examples

>>> from oroscope import site_searcher as ss
>>> round(ss.estimate_peak_memory_gb(1981, 3061, downsample_factor=1), 2)
0.77
>>> round(ss.estimate_peak_memory_gb(10204, 12603, downsample_factor=4), 2)
5.08

Downsampling is the weaker of the two levers at this scale, and striding the stronger, because the candidates are taken on the native grid either way:

>>> round(ss.estimate_peak_memory_gb(10204, 12603, downsample_factor=1), 2)
7.21
>>> round(ss.estimate_peak_memory_gb(10204, 12603, downsample_factor=4,
...                                  candidate_stride=10), 2)
2.83

oroscope.site_searcher.estimate_visualisation_memory_gb(rows, cols, downsample_factor=1, combine=False)

Memory the map costs, in GiB, which is not what the search costs.

estimate_peak_memory_gb() models the candidate and scoring arrays and nothing else, because those are what a search allocates. The map is a separate peak that lands on top of them at the very end, and it is the stage that has actually been failing: three runs in one session finished their searches and then died drawing the picture, having written the JSON and the GeoTIFF first. A pre-flight that models only the search will size a cap that cannot survive the run.

The map renders at viz_ds = downsample_factor * 2, so its raster is rows/(2d) x cols/(2d). Measured on this machine over 0.35 to 5.6 Mpx of viz raster – shading through LightSource.shade and saving at 150 dpi:

viz raster	peak RSS above idle
700 x 500 (0.35 Mpx)	126 MB
1400 x 1000 (1.4 Mpx)	263 MB
1981 x 1441 (2.85 Mpx)	539 MB
2800 x 2000 (5.6 Mpx)	959 MB

which is ~190 bytes per viz pixel once matplotlib’s own ~130 MB of import and canvas is taken out. Both terms are carried here.

The combination is a different, larger stage, and was not modelled at all. combine_experiments renders at the mask’s own resolution rather than at viz_ds, so its raster is four times the search map’s, and it draws the two experiments and their overlap on top of the relief. Sized by the search map’s numbers it looked cheap, ran last, and is what actually kept failing. Measured over 2.0 to 12.5 Mpx of raster, with the overlay composited in float32:

combine raster	peak RSS
1276 x 1576 (2.0 Mpx)	0.46 GiB
1806 x 1740 (3.1 Mpx)	0.63 GiB
2551 x 3151 (8.0 Mpx)	1.31 GiB
3200 x 3900 (12.5 Mpx)	1.95 GiB

which fits 182 MiB + 152.2 bytes per pixel to within 1.2% at every point. Four points rather than the one the search estimator is calibrated on.

Parameters:

rows, colsint

Full-resolution DEM dimensions.

downsample_factorint, optional

As the pipeline’s. The search map uses twice it; the combination uses it as-is.

combinebool, optional

Estimate the combination overlay instead of a search’s own map. It renders at four times the pixels and costs about 2.7x as much at DEM scale.

Returns:

float

Estimated peak for that one figure, in GiB.

Examples

>>> from oroscope import site_searcher as ss
>>> round(ss.estimate_visualisation_memory_gb(3961, 2881, 1), 2)
0.63

The same DEM’s combination, which the search’s own estimate does not cover:

>>> round(ss.estimate_visualisation_memory_gb(3961, 2881, 1, combine=True), 2)
1.8

oroscope.site_searcher.apply_memory_cap(max_memory_gb)

Caps this process’s address space, so a runaway fails instead of taking the machine.

Without a cap, an over-large search does not fail: it grows until the kernel’s OOM killer chooses a victim, which may well be something else the user cares about. A MemoryError inside this process is a far better outcome than a dead editor, and it names the run that caused it.

Parameters:

max_memory_gbfloat or None

Ceiling in GiB. None or non-positive leaves the limit alone.

Returns:

bool

Whether a cap was applied. False on platforms without RLIMIT_AS.

Notes

RLIMIT_AS limits virtual address space, which is larger than resident memory: numba and BLAS reserve address space they never touch. Set it generously – a little above physical RAM is usually right – or it will refuse runs that would have fit.

oroscope.site_searcher.available_memory_gb()

Memory the system can give us right now, in GiB, or None if it cannot be told.

Reads MemAvailable from /proc/meminfo, which accounts for reclaimable page cache and so is the figure that matters; free alone understates it badly on a machine that has been running a while.

Returns:

float or None

Available memory in GiB, or None on a platform without /proc/meminfo.

oroscope.site_searcher.preflight_memory(dem_path, downsample_factor=1, candidate_stride=5, max_memory_gb=None, quiet=False, refuse=False, combine=False)

Estimates what a search will need, says so, and caps the address space.

This ran only inside main(), so a library user – a sweep, a notebook, a service – got neither the warning nor the cap unless they knew to ask for both. That is exactly the caller most likely to need them: a ten-point sweep once reached 6.9 GB and was killed by the OOM killer, taking other work with it.

Parameters:

dem_pathstr

DEM whose dimensions set the estimate. An unreadable file skips the estimate but not the cap.

downsample_factorint, optional

As the pipeline’s. Dominates the estimate: the labelling arrays scale as its inverse square.

candidate_strideint, optional

As the pipeline’s.

max_memory_gbfloat, optional

Ceiling in GiB. None uses 80% of what the system reports available, so the cap bites before the kernel does; 0 disables capping entirely.

quietbool, optional

Suppress the printed report, keeping the cap and the returned numbers.

refusebool, optional

Raise MemoryError instead of warning when the estimate exceeds REFUSE_FRACTION of what is available. Default False, which only warns — but warning is what it did while three runs died anyway, so a caller that can afford to stop should pass True.

combinebool, optional

Also account for a combination overlay after the searches. It renders at four times the search map’s pixels and was not modelled at all, so a region could clear this pre-flight and then die at the last step — which is what happened, twice. The stages are sequential rather than simultaneous, so the number judged is max(search + map, combine) rather than their sum.

Returns:

dict

{"estimate_gb", "search_gb", "visualisation_gb", "combine_gb", "available_gb", "cap_gb", "capped", "cap_exceeds_available"}. estimate_gb is the largest stage the run will pass through, and is None when the DEM could not be measured. combine_gb is None unless combine was asked for.

Raises:

MemoryError

When refuse is set and the estimate exceeds REFUSE_FRACTION of the available memory.

Notes

The estimate is rough by construction – it assumes a survival fraction the topographic screen has not yet measured – so an over-large search is warned about rather than refused.

A cap above what is available is not a cap. RLIMIT_AS protects the machine only if the process hits it before the kernel runs out of memory to give; set above the available figure, the OOM killer arrives first and the limit never fires. This is easy to do by accident, because the two pieces of advice around it pull opposite ways: the estimate counts only anonymous memory, so on a large DEM the cap has to clear it by the size of the memory-mapped file (roadmap 6.46), and raising it for that reason can quietly carry it past the available figure. It happened here – a 339 Mpx search capped at 13.0 GiB on a machine with 8.0 GiB available, which took the machine down. When the two constraints cannot both be met, the configuration does not fit and the answer is a coarser one, not a bigger number. Warned about rather than refused, for the same reason as the estimate: the available figure moves while the run is being set up.

oroscope.site_searcher.emit_explanation(results, run_output_dir=None, print_it=True)

Composes the run’s plain-language summary, prints it, and saves it beside the run.

A thin wrapper over explain.explain_results(): the words themselves are that function’s business, so they can be regenerated from an old results file without this one. What is added here is the placement – last on the console, so it is what a reader is left with, and in explanation.txt so the run can be handed on without the terminal it was run in.

Failures are reported and swallowed. A summary that cannot be written is not a reason to lose a search that already succeeded.

Parameters:

resultsdict

The run’s results. Gains an "explanation" key.

run_output_dirstr, optional

Directory to write explanation.txt into. Omitted writes no file.

print_itbool, optional

Whether to print. False still composes and stores the text.

Returns:

str or None

The summary, or None if it could not be composed.

oroscope.site_searcher.resolve_config_paths(config, config_dir, quiet=False)

Makes a configuration’s relative paths absolute, against the config’s own directory.

A configuration that says "dem_path": "../input/dem/colca.tif" is describing where the DEM sits relative to itself, which is the only thing it can know. The pipeline resolved it against the working directory instead, so the bundled configs ran only from src/ and produced a FileNotFoundError anywhere else – the long-standing “you must cd src first” wart.

A path that does not resolve against the configuration’s directory but does resolve against the working directory is left alone, with a warning: that is the old behaviour, and silently breaking someone’s working setup to fix a wart is a poor trade. Absolute paths are untouched.

Parameters:

configdict

A configuration mapping. Not modified; a copy is returned.

config_dirstr

Directory holding the configuration file.

quietbool, optional

Suppress the warning about a working-directory-relative fallback.

Returns:

dict

A copy with the path keys resolved.

Examples

The repository’s own layout is why this is safe to change: config/ and src/ are both one level below the root, so ../input/dem/colca.tif names the same file read either way, and no shipped configuration has to change.

>>> import os
>>> from oroscope import site_searcher as ss
>>> a = os.path.normpath(os.path.join("/repo/config", "../input/dem/colca.tif"))
>>> b = os.path.normpath(os.path.join("/repo/src", "../input/dem/colca.tif"))
>>> a == b == "/repo/input/dem/colca.tif"
True

An absolute path is left as it is:

>>> cfg = ss.resolve_config_paths({"dem_path": "/data/x.tif"}, "/repo/config")
>>> cfg["dem_path"]
'/data/x.tif'

oroscope.site_searcher.stride_gap_m(candidate_stride, cell_size_y_m)

Distance between kept candidates, in metres, after striding.

candidate_stride subsamples the list of surviving pixels rather than the map, so the gap it leaves is a stride’s worth of pixels along a scanline.

Parameters:

candidate_strideint

Keeps every Nth surviving pixel.

cell_size_y_mfloat

Metric pixel size, N-S.

Returns:

float

Gap between kept candidates, in metres.

oroscope.site_searcher.closing_element_m(gap_close_km, antenna_spacing_km)

Size of the morphological closing element in metres, defaulting as the pipeline does.

Parameters:

gap_close_kmfloat or None

Closing element in km. None defaults to antenna_spacing_km.

antenna_spacing_kmfloat

Detector spacing, which the closing element defaults to.

Returns:

float

Closing element size, in metres.

oroscope.site_searcher.warn_stride_outruns_closing(candidate_stride, cell_size_y_m, gap_close_km, antenna_spacing_km, quiet=False)

Warns when the closing element is too small to bridge the gaps striding leaves.

Striding is unbiased in acceptance – measured at both scales, 58.414% against 58.415% for GRAND and 75.750% against 75.736% for TAMBO – so it is tempting to treat it as free. It is not. Accepted pixels are marked one in candidate_stride, and the mask is then closed morphologically before areas are measured. If the closing element is smaller than the gap the stride leaves, the mask never reconnects: it stays a scatter of isolated pixels, small regions fall below the size and capacity thresholds, and the reported area collapses.

Measured at Colca against a 154 m stride-5 gap: at TAMBO’s published 150 m element, 203.0 km² reported against 307.2 km² at stride 1, a 1.51x under-report. At the 100 m element this warning was written for – three pixels against a five-pixel gap – the same comparison was 83.6 against 396.9 km², 4.75x. Acceptance is identical to three decimal places either way. The same run at GRAND’s 1 km element, 32 px against the same gap, is unaffected, which is why this went unnoticed.

The 150 m figure is the one to quote, and the drop from 4.75x to 1.51x is the whole argument for heeding this warning: two pixels of element is the difference between the two sides of a cliff.

The rule is simply that the element must outrun the gap. Raise gap_close_km, lower candidate_stride, or accept the area as a lower bound and say so.

Parameters:

candidate_strideint

Keeps every Nth surviving pixel.

cell_size_y_mfloat

Metric pixel size, N-S.

gap_close_kmfloat or None

Closing element in km. None defaults to antenna_spacing_km.

antenna_spacing_kmfloat

Detector spacing, which the closing element defaults to.

quietbool, optional

Suppress the printed warning, keeping the returned verdict.

Returns:

dict or None

{"gap_m", "element_m", "ratio"} when the element cannot bridge the gap, and None when it can.

Examples

GRAND’s 1 km element easily bridges a stride-5 gap at 30 m pixels:

>>> from oroscope import site_searcher as ss
>>> ss.warn_stride_outruns_closing(5, 30.72, None, 1.0, quiet=True) is None
True

TAMBO’s element does not – not at the old 100 m and not at the published 150 m either, which is why striding costs it so much area:

>>> r = ss.warn_stride_outruns_closing(5, 30.72, None, 0.15, quiet=True)
>>> round(r["gap_m"]), round(r["element_m"]), round(r["ratio"], 2)
(154, 150, 1.02)

>>> r = ss.warn_stride_outruns_closing(5, 30.72, None, 0.1, quiet=True)
>>> round(r["ratio"], 2)
1.54

Closing disabled entirely is not this failure, so it does not warn:

>>> ss.warn_stride_outruns_closing(5, 30.72, 0.0, 0.1, quiet=True) is None
True

oroscope.site_searcher.add_scale_bar(ax, km_per_x_unit, fraction=0.22, colour='black')

Draws a kilometre scale bar on a map, and returns the length it chose.

A map axis labelled in degrees or in pixels does not tell a reader how far anything is, and neither unit converts to distance without knowing where on the Earth it sits: a degree of longitude at Arequipa is 4% shorter than a degree of latitude, and a pixel is whatever the DEM says it is.

Taking km_per_x_unit rather than a latitude keeps one function usable by both maps this project writes – the search map, whose axes are pixels, and the combination overlay, whose axes are degrees.

Parameters:

axmatplotlib.axes.Axes

Axes to draw on. Its limits must already be final; the bar is placed relative to them.

km_per_x_unitfloat

Kilometres per unit of the x axis. For degrees of longitude that is 111.32 * cos(latitude); for pixels it is the metric pixel size over 1000.

fractionfloat, optional

Roughly what fraction of the map width the bar should span, before rounding to a human number.

colourstr, optional

Bar colour. The default reads on both terrain and shaded relief.

Returns:

float

Length of the bar drawn, in km. Always 1, 2 or 5 times a power of ten.

Examples

A two-degree map at 16 degrees south is about 214 km wide, so it gets a 50 km bar:

>>> import matplotlib; matplotlib.use("Agg")
>>> import matplotlib.pyplot as plt, numpy as np
>>> from oroscope import site_searcher as ss
>>> fig, ax = plt.subplots()
>>> _ = ax.set_xlim(-73.0, -71.0); _ = ax.set_ylim(-17.0, -15.0)
>>> ss.add_scale_bar(ax, 111.32 * np.cos(np.radians(-16.0)))
50.0

The same function on a pixel axis, 3061 pixels of 30 m:

>>> _ = ax.set_xlim(0, 3061); _ = ax.set_ylim(1981, 0)
>>> ss.add_scale_bar(ax, 0.030)
20.0
>>> plt.close(fig)

oroscope.site_searcher.altitude_limits(elevation, low_percentile=0.5, high_percentile=99.8)

Altitude range for a colour scale, from the DEM rather than from a constant.

A fixed 0-6000 m scale spends most of its range on altitudes a given DEM does not contain: the Colca crop runs 1500-6300 m, so half the colour bar described ground that is not in the picture and the relief that is there got half the contrast it could have had.

Percentiles rather than the extremes, because one spurious pixel – a nodata sentinel, a spike – otherwise sets the whole scale. Water is excluded from the upper end and pinned to the bottom, so a coastal DEM does not spend a third of its range on ocean.

Parameters:

elevationndarray

Elevation in metres. NaN is ignored.

low_percentile, high_percentilefloat, optional

Percentiles of the land pixels to clip to.

Returns:

tuple of float

(vmin, vmax) in metres, rounded outward to a round number.

Examples

>>> import numpy as np
>>> from oroscope import site_searcher as ss
>>> z = np.linspace(1500.0, 6300.0, 1000)
>>> ss.altitude_limits(z)
(1500.0, 6300.0)

Ocean does not drag the floor down with it:

>>> z = np.concatenate([np.zeros(500), np.linspace(2000.0, 5000.0, 500)])
>>> ss.altitude_limits(z)
(0.0, 5000.0)

oroscope.site_searcher.add_north_arrow(ax, x=0.965, y=0.955, size=0.055)

Draws a north arrow in axes coordinates.

Both maps this project writes are north-up, so the arrow is a convention rather than information – but a map without one asks the reader to assume, and a map in a talk gets read by people who did not make it.

Parameters:

axmatplotlib.axes.Axes

Axes to draw on.

x, yfloat, optional

Position of the arrow’s tip, in axes coordinates.

sizefloat, optional

Length of the arrow, as a fraction of the axes height.

oroscope.site_searcher.attach_colorbar(fig, ax, mappable, label, width='2.6%', pad=0.12, **kwargs)

Adds a colour bar whose height matches the plot panel exactly.

fig.colorbar(..., fraction=...) sizes the bar as a fraction of the figure, so on a map whose aspect is set by its data – which every map here is – the bar overshoots the panel top and bottom by however much the axes shrank to fit. Taking the space out of the axes’ own divider instead ties the two together whatever the aspect turns out to be.

Parameters:

figmatplotlib.figure.Figure

Figure the axes belong to.

axmatplotlib.axes.Axes

Axes to take the space from.

mappablematplotlib.cm.ScalarMappable

What the bar describes.

labelstr

Axis label for the bar.

widthstr, optional

Bar width, as a percentage of the axes width.

padfloat, optional

Gap between panel and bar, in inches.

**kwargs

Passed to fig.colorbar – extend, for instance.

Returns:

matplotlib.colorbar.Colorbar

oroscope.site_searcher.resolve_settlements(value, bounds=None)

Settlements to mark on a map: a preset name, an explicit list, or "auto".

"auto" picks whichever curated list has points inside the map, which is what a reader wants without having to say so, and marks nothing when neither does rather than guessing at a region the project has no coordinates for.

No coordinates are invented here. The lists are the named places already curated as RFI zones; adding more means supplying them, not asking this function to remember them.

Parameters:

valuestr, list, or None

"auto", a preset name, "none"/None, or a list of (latitude, longitude, name).

boundstuple, optional

(south, north, west, east) in degrees, used by "auto".

Returns:

list

(latitude, longitude, name) triples, possibly empty.

Examples

>>> from oroscope import site_searcher as ss
>>> [n for _, _, n in ss.resolve_settlements("arequipa")][:2]
['Arequipa', 'Majes']

Auto picks the list with points on the map, and nothing when none has any:

>>> arequipa = (-17.5, -14.5, -73.6, -70.3)
>>> len(ss.resolve_settlements("auto", arequipa))
5
>>> ss.resolve_settlements("auto", (40.0, 42.0, 0.0, 2.0))
[]
>>> ss.resolve_settlements(None)
[]

oroscope.site_searcher.add_settlements(ax, settlements, to_axes, fontsize=9, max_labels=6, max_markers=25)

Marks settlements, with their names, wherever they fall inside the axes.

Parameters:

axmatplotlib.axes.Axes

Axes to draw on. Its limits must already be final.

settlementssequence

(latitude, longitude, name) triples.

to_axescallable

to_axes(latitude, longitude) -> (x, y) in the axes’ own units, so this serves both the pixel-coordinate search map and the degree-coordinate overlay.

fontsizeint, optional

Fallback label size, for entries carrying no place class.

max_labelsint, optional

How many to name.

max_markersint, optional

How many to mark at all. Entries arrive most-important-first – cities, then towns, then villages, by population within each – so both caps take the ones worth showing. Over the full Arequipa DEM there are 1,268 places, and drawing every one buries the result under a thousand white dots; over one canyon there are sixty, and naming them all is not a map either.

Returns:

int

How many were drawn. Points outside the axes are skipped rather than clipped, so a label cannot appear at the edge pointing at nothing.

oroscope.site_searcher.add_roads(ax, roads, to_axes, colour='#00B33C', alpha=0.9, scale=1.0)

Draws road geometry, as context rather than as a result.

Access is the question a site count cannot answer: a canyon wall that takes a hundred detectors means something different with a road along the rim than with the nearest track forty kilometres away. Roads are drawn in a neutral dark line rather than a colour, because on these maps colour belongs to the categories and a road is not one of them – it is what the reader measures a category against.

Uses one LineCollection per class rather than a call per road. The Arequipa DEM carries 8,780 roads and 367,000 vertices; drawn individually that is slow enough to notice, and matplotlib renders the collection in one pass.

Parameters:

axmatplotlib.axes.Axes

Axes to draw on. Its limits must already be final.

roadsdict or None

As fetch_roads.load_roads() returns. None draws nothing, so a missing road file leaves a map without roads rather than without a map.

to_axescallable

to_axes(latitude, longitude) -> (x, y) in the axes’ own units, so this serves both the pixel-coordinate search map and the degree-coordinate overlay.

colourstr, optional

Line colour.

alphafloat, optional

Line opacity.

scalefloat, optional

Multiplies every width, for maps drawn at a different size.

Returns:

int

How many roads were drawn.

Explaining a run

Turning a results dictionary into an account of what was found and why: which constraint set the size of the answer, what held the surviving sites back, and which of the numbers are assumptions. Pure, so an old results file can be explained months later with no DEM and no pipeline.

Turning a results file into an account of what was found and why.

Everything this module says is already in the results JSON. What is missing there is the story: which constraint did the work, which sites survived and what weakened them, and which of the numbers on the page are choices rather than measurements. A reader who assembles that themselves gets it wrong in predictable ways – most often by reading the reported area as physics-accepted area, which it is not (see Assumptions and limitations, and the 2.35x measured at Colca).

The entry point is explain_results(), which takes the results dictionary and returns a string. It runs nothing, opens nothing and needs no DEM, so the pipeline, the library and a test can all use the same words, and the summary can be regenerated from an old run’s JSON long after the run.

Deliberately plain text: these summaries are meant to be handed to other people, and ANSI colour does not survive a paste into an email.

oroscope.explain.explain_results(results, provenance=None)

Writes a human-readable account of one search: what was found, and why.

Everything it reports is already in the results dictionary. The value added is the reading: which stage of the funnel actually set the size of the answer, what held the surviving sites back, and which of the numbers are assumptions. Those are the three things a reader gets wrong without help, and a run that is going to be handed to someone else needs all three on the same page as the result.

Pure: it opens no files, runs nothing, and needs no DEM. An old run’s JSON can be explained months later, and the pipeline, the library and the tests all get the same words.

Parameters:

resultsdict

A results dictionary as returned by site_searcher.find_grand_regions_interactive() and written to the run’s results JSON. Missing sections are tolerated: each is reported as absent rather than raising, so a partial or older file still explains.

provenancedict, optional

The matching provenance.json contents – git commit, DEM checksum, package versions. Adds the reproducibility block when given.

Returns:

str

The summary, as plain text with no ANSI colour, ready to print or to save beside the results.

Examples

>>> from oroscope import explain
>>> text = explain.explain_results({
...     "funnel": {"DEM pixels": 1000, "slope 3-25 deg": 900,
...                "directions accepted": 12},
...     "results": {"total_sites": 0, "total_capacity": 0, "sites": []}})
>>> "WHERE THE CANDIDATES WENT" in text
True
>>> "directions accepted" in text
True

oroscope.explain.explain_combination(report, runs=None)

Explains an overlay of two or more searches: who can share ground with whom, and why.

The combined report gives a joint area and a Jaccard index. Neither says why the number is what it is, and the reason is usually not about neutrinos at all: a pixel has one slope, and every experiment deployed on it must accept that slope. Measured at Colca, the entire co-location result follows from GRAND’s 3-25 degree deployable band against a canyon’s ~40 degree walls.

So where the runs’ own parameters are available, this compares their screening bands and names the one that limits the sharing.

Parameters:

reportdict

A combined_report.json, as combine_experiments writes it.

runsdict, optional

Label to results dictionary, for the runs being combined. Supplying them adds the constraint comparison; without them only the areas are explained.

Returns:

str

The summary, as plain text.

Examples

>>> from oroscope import explain
>>> text = explain.explain_combination({
...     "runs": [{"label": "A", "area_km2": 100.0, "pixels": 10,
...               "area_in_joint_km2": 20.0, "fraction_of_own_area_in_joint": 0.2,
...               "reported_sites": 1, "reported_capacity": 500}],
...     "joint": {"area_km2": 20.0}, "union": {"area_km2": 100.0},
...     "joint_requires": ["A"], "pairwise_overlap": {}})
>>> "WHERE THESE EXPERIMENTS CAN SHARE GROUND" in text
True

oroscope.explain.binding_constraint(funnel)

Finds the funnel stage that removed the largest share of what reached it.

This is the single most useful sentence a summary can offer, and it matters most when a search returns little or nothing: the stage where the survivor count collapses is the constraint responsible, and every other explanation is a guess.

Two stages are excluded by construction. kept by stride is a deliberate subsample whose acceptance is unbiased, so calling it a constraint would name the same answer on nearly every run; after gap closing adds pixels rather than removing them.

Parameters:

funneldict

Ordered stage-name to survivor-count mapping, as written to results["funnel"].

Returns:

dict or None

{"stage", "survivors", "before", "kept_fraction", "knob", "fatal"}, or None when the funnel has fewer than two stages. fatal is True when the stage left nothing at all.

Examples

>>> from oroscope import explain
>>> f = {"DEM pixels": 1000, "slope 3-25 deg": 800, "directions accepted": 40}
>>> b = explain.binding_constraint(f)
>>> b["stage"], round(b["kept_fraction"], 3), b["fatal"]
('directions accepted', 0.05, False)
>>> explain.binding_constraint({"DEM pixels": 1000}) is None
True

oroscope.explain.weakest_component(arrival_scan, statistic='p50')

Names the score component that held a site back, and its value.

The score is a product of components each in [0, 1] and each named, so a low total can be attributed rather than merely reported. The lowest median component is the one to look at first: under a product composition it bounds the total from above.

Parameters:

arrival_scandict

A site’s arrival_scan record.

statisticstr, optional

Which per-site statistic to compare, "mean", "p50" or "p90".

Returns:

tuple or None

(name, value) for the lowest-scoring component, or None when the record carries no components – which is the case for runs made before they were stored.

Examples

>>> from oroscope import explain
>>> rec = {"score_p50": 0.2, "score_decay_p50": 0.9, "score_shower_p50": 0.22}
>>> explain.weakest_component(rec)
('shower', 0.22)
>>> explain.weakest_component({"score_p50": 0.2}) is None
True

oroscope.explain.site_strengths(arrival_scan, statistic='p50', threshold=0.75)

Why a site is good: the criteria it satisfies well, and the measurement behind each.

The mirror of weakest_component(), and the more useful half when a site has been selected. “Site 3555 scored 0.55” says nothing a reader can act on; “it sees 0.36 sr of usable sky across a 3.1 km gap with 780,000 g/cm² of rock behind it, and every criterion but the accepted solid angle is satisfied outright” says what the ground is actually like.

Parameters:

arrival_scandict

A site’s arrival_scan record.

statisticstr, optional

Which per-site statistic to read, "mean", "p50" or "p90".

thresholdfloat, optional

Score at or above which a component counts as satisfied. 0.75 rather than 1.0 because a band score falls off smoothly either side of its plateau, so insisting on exactly 1 would report nothing on most real sites.

Returns:

list of dict

One entry per satisfied component, strongest first, each with name, label, score, means and – where the record carries the observable behind it – evidence. Empty when the record has no components.

Examples

>>> from oroscope import explain
>>> rec = {"score_solid_angle_p50": 0.9, "score_depth_p50": 1.0,
...        "solid_angle_sr_p50": 1.08, "max_depth_gcm2_p50": 784440.0}
>>> [s["name"] for s in explain.site_strengths(rec)]
['depth', 'solid_angle']
>>> explain.site_strengths(rec)[1]["evidence"]
'1.08 sr'

oroscope.explain.constraint_overlap(params_a, params_b, bands=None)

Where two experiments’ screening bands agree, and by how little.

This is what decides whether two experiments can share ground, and it is decided before any arrival geometry is considered: a pixel has one slope, and both experiments must accept it. Measured at Colca, that is the whole story – GRAND’s 3-25 degree deployable band against a canyon’s ~40 degree walls leaves a 20-25 degree sliver, and the joint area follows from that rather than from anything about neutrinos.

Parameters:

params_a, params_bdict

The two runs’ recorded parameters blocks.

bandssequence, optional

Which bands to compare, as (label, low key, high key, unit). Defaults to the properties of the ground itself – slope, altitude, aspect – which are the only ones both experiments must agree on. Pass _VIEWING_BANDS to compare what each asks of the view instead, which need not agree at all.

Returns:

list of dict

One entry per band both runs recorded: label, a, b, overlap (a (low, high) pair or None), width, and share_of_narrower – the fraction of the tighter of the two bands that the overlap covers, which is the number that says how much room there is to share.

Examples

>>> from oroscope import explain
>>> a = {"min_slope_deg": 3.0, "max_slope_deg": 25.0}
>>> b = {"min_slope_deg": 20.0, "max_slope_deg": 60.0}
>>> band = explain.constraint_overlap(a, b)[0]
>>> band["overlap"], round(band["share_of_narrower"], 3)
((20.0, 25.0), 0.227)

oroscope.explain.closing_inflation(funnel, candidate_stride=1)

How much morphological closing grew this run’s mask, measured from its own funnel.

The 2.35x quoted from Colca is a property of that terrain and a 1 km element, not a constant. This run has the number in it: the stage before closing counts the accepted candidates, closing counts the pixels after, and the only correction needed between them is the stride – which samples one candidate in candidate_stride and was measured unbiased.

A ratio below 1 is not a bug and is worth reading carefully: it means the closing element was too small to bridge the gaps striding left, so the mask under-reports the accepted set rather than inflating it. That happens when the element is a few pixels across, as it is at a 100 m detector spacing.

Parameters:

funneldict

Ordered stage-name to survivor-count mapping.

candidate_strideint, optional

The run’s candidate_stride, to scale the accepted count back to full resolution.

Returns:

float or None

Closed pixels divided by estimated accepted pixels, or None when the funnel does not record both.

Examples

>>> from oroscope import explain
>>> f = {"directions accepted": 100, "after gap closing": 450}
>>> round(explain.closing_inflation(f, candidate_stride=5), 2)
0.9

oroscope.explain.selected_sites(results)

The sites actually in the result, separated from the ones that merely qualified.

results["results"]["sites"] lists every site that cleared the area and capacity thresholds, which with stop_at_target is more than were selected: selection walks the capacity-sorted list until the target is met and stops. total_sites, total_capacity and the exported mask all cover the selection only, so summing the list over-reports area and site count against every other number in the file.

Parameters:

resultsdict

A results dictionary.

Returns:

selectedlist of dict

Sites in the result.

rejectedlist of dict

Sites that qualified but were not selected. Usually empty.

Notes

Prefers each record’s selected flag. Files written before that flag existed fall back to the first total_sites entries, which is exact: the list is sorted by capacity and selection takes it in order.

Examples

>>> from oroscope import explain
>>> r = {"results": {"total_sites": 1, "sites": [
...     {"site_id": 2, "capacity_exact": 252, "selected": True},
...     {"site_id": 1, "capacity_exact": 36, "selected": False}]}}
>>> chosen, rest = explain.selected_sites(r)
>>> [s["site_id"] for s in chosen], [s["site_id"] for s in rest]
([2], [1])

Tools

Cuts a lat/lon window out of a DEM into a smaller GeoTIFF.

Two reasons this exists. A regional study wants one canyon or one plateau rather than the whole 10204x12603 tile the download gives you. And comparing experiments requires them to cover identical ground: the combiner overlays masks pixel for pixel, so GRAND and TAMBO have to be run against the same crop, not two crops that happen to overlap.

The window is snapped outward to whole pixels and the tiepoint recomputed for the crop’s own corner, so the output is georeferenced in its own right.

python crop_dem.py ../input/dem/arequipa_SRTMGL1.tif ../input/dem/colca.tif –north -15.30 –south -15.80 –west -72.40 –east -71.60

oroscope.crop_dem.crop(src: str, dst: str, north: float, south: float, west: float, east: float) → dict

Cuts a geographic window out of a DEM, writing a GeoTIFF that stands on its own.

The point of cropping is not disk space, it is sampling: a department at downsample_factor 4 and candidate_stride 5 costs area and fragments a thin mask, while a crop small enough to run at 1 and 1 does neither. Every unbiased number this project quotes comes from a crop made here.

The window is snapped outward to whole pixels, so the result contains the requested box rather than approximating it, and is clipped to the DEM. A window that misses the DEM entirely raises rather than writing an empty file.

The crop carries its own tiepoint, not the parent’s. That is what lets it be searched with no reference back to where it came from — and it is why a crop’s origin_lat/origin_lon differ from the parent’s, which is correct and not drift.

Parameters:

srcstr

The GeoTIFF to cut from.

dststr

Where to write the crop.

north, southfloat

Latitude bounds in degrees. north is the larger (less negative) value.

west, eastfloat

Longitude bounds in degrees.

Returns:

dict

path, the grid as rows and cols, the crop’s own origin_lat and origin_lon with the matching south and east edges, cell_size_deg, and the elevation range as z_min and z_max.

Raises:

SystemExit

If the requested window does not overlap the DEM, reporting what the DEM actually covers.

Examples

Cutting the Colca crop out of the Arequipa department DEM, which is how input/dem/colca.tif was made:

from oroscope import crop_dem

info = crop_dem.crop("input/dem/arequipa_SRTMGL1.tif", "input/dem/colca.tif",
                     north=-15.30, south=-15.85, west=-72.40, east=-71.55)
print(info["rows"], info["cols"], info["z_min"], info["z_max"])

oroscope.crop_dem.read_geo(path: str) → tuple[float, float, float, float, int, int]

Pixel size in degrees and the north-west corner, from the GeoTIFF tags.

Parameters:

pathstr

Path to a geographic GeoTIFF carrying ModelPixelScaleTag and ModelTiepointTag.

Returns:

tuple

(cell_x_deg, cell_y_deg, lon0, lat0, rows, cols), with lon0/lat0 the north-west corner.

Combines the results of two or more experiment searches over the same ground.

GRAND and TAMBO ask the same structural question – from this patch of ground, is there a target surface at the right range, in the right direction, with the right matter behind it? – and differ in their numbers rather than their structure. So each experiment is one run of the searcher with its own configuration, and combining them is an overlay of the masks those runs produce.

Three questions get different answers, and all three are worth reporting:

joint terrain that satisfies every experiment at once. This is the co-location

case: one site, one road, one power feed, two experiments.

union terrain that satisfies any of them. This is the coverage case: how much of

the region is useful to the programme as a whole.

each what each experiment gets on its own, and what it would lose by being

confined to the joint area.

The inputs must be pixel-aligned: same shape, same pixel size, same corner. That is not a detail to paper over – two runs on differently-cropped DEMs would silently overlay the wrong ground – so it is checked and refused rather than resampled.

python combine_experiments.py ../output/grand_colca_config ../output/tambo_colca_config

–labels GRAND TAMBO –out ../output/combined_colca

oroscope.combine_experiments.load_run(run_dir: str) → dict

Loads one search’s mask, georeferencing and results JSON.

The searcher writes the mask as a downsampled GeoTIFF beside a world file and a results JSON, all sharing a base name.

Parameters:

run_dirstr

A run’s output directory.

Returns:

dict

dir, tif, mask, world and results. results is None when no results JSON is present.

Raises:

SystemExit

If no mask GeoTIFF is found, or if the world file beside it is missing – without which alignment cannot be confirmed.

oroscope.combine_experiments.check_alignment(runs: list[dict]) → None

Refuses to overlay masks that do not describe the same ground.

Comparing shapes is not enough: two crops of the same size taken from different corners would overlay cleanly and mean nothing. The world file’s pixel size and upper-left corner are what actually pin the ground down.

Parameters:

runslist of dict

Loaded runs, from load_run(). The first is taken as the reference.

Raises:

SystemExit

If any run differs from the reference in shape, pixel size or corner. Refusing is deliberate: resampling would silently compare the wrong terrain.

oroscope.combine_experiments.read_world_file(tfw_path: str) → tuple[float, ...]

Reads the six affine terms of an ESRI world file.

Parameters:

tfw_pathstr

Path to the .tfw file.

Returns:

tuple of float

(pixel_size_x, rot_y, rot_x, pixel_size_y, upper_left_x, upper_left_y).

Raises:

ValueError

If the file does not hold exactly six terms.

oroscope.combine_experiments.pixel_area_km2(world: tuple[float, ...], reference_latitude_deg: float) → float

Ground area of one mask pixel, in km^2.

The world file is in degrees, so the east-west size shrinks with the cosine of the latitude. This uses the same convention as the searcher’s own grid geometry.

Parameters:

worldtuple of float

The six affine terms, from read_world_file().

reference_latitude_degfloat

Latitude at which to evaluate the east-west pixel size, in degrees. Normally the centre of the map.

Returns:

float

Area of one pixel, in km^2.

oroscope.combine_experiments.capacity_of(results)

The capacity a run reported, or None when it did not report one.

search_mode: single writes the string 'N/A' rather than a number, so int() raises ValueError – which was not caught, and took the whole combination down with it. A run that does not report a capacity is an ordinary case here, not an error: the overlay is about ground, not detectors.

Parameters:

resultsdict or None

A run’s parsed results JSON, or None when the directory held none.

Returns:

int or None

The reported capacity, or None when the run did not report one.

oroscope.combine_experiments.colocation_capacity(partner_mask, world, centre_lat, centre_lon, radii_km=(5.0, 10.0, 20.0, 40.0), spacing_km=1.0, grid_type='hex')

How much partner ground lies within reach of a site, and what it could hold.

The question this answers is the deployment one rather than the site-finding one: given a site chosen for one experiment, is there enough ground nearby for the other? It is deliberately not an intersection. The joint mask is the ground both experiments accept at the same pixel, and a partner array does not have to stand there – measured on the Cajatambo crop, the GRAND-viable ground inside the joint mask was 22,577 fragments of which exactly one was large enough for a single 1 km lattice cell, while 976 km2 of perfectly good GRAND ground lay within 20 km (roadmap 6.52). An optimiser pointed at the intersection reports the partner array impossible; the site is fine.

What couples two arrays is a shared line of sight to the same massif, not a shared footprint. GRAND’s own targets are 10-40 km away, so a partner antenna 20 km from a TAMBO strip is still watching the same wall.

Parameters:

partner_maskndarray

2-D boolean mask of the other experiment’s viable ground, on the mask grid.

worldtuple of float

The six affine terms of that mask’s world file, from read_world_file().

centre_lat, centre_lonfloat

The site to measure from, in degrees. A site’s center_lat/center_lon from a results JSON is the usual source.

radii_kmiterable of float, optional

Distances to report, in km.

spacing_kmfloat, optional

Detector spacing of the partner array, in km. GRAND’s is 1.0, TAMBO’s 0.1.

grid_type{‘hex’, ‘square’}, optional

Lattice the partner array uses.

Returns:

list of dict

One entry per radius, ascending, with radius_km, area_km2 and capacity. capacity counts positions on an anchored lattice, as count_grid_capacity() does, so it is an estimate for an arbitrarily placed array rather than the best achievable packing.

Notes

Distances are great-circle-free: a local flat approximation with the east-west degree shrunk by the cosine of the site’s latitude, which is what the rest of this project uses and is good to much better than a pixel over the tens of kilometres involved.

Examples

>>> import numpy as np
>>> from oroscope import combine_experiments as ce
>>> mask = np.ones((200, 200), dtype=bool)
>>> world = (1/3600, 0.0, 0.0, -1/3600, -77.0, -10.0)
>>> rows = ce.colocation_capacity(mask, world, -10.02, -76.98, radii_km=(2.0,))
>>> rows[0]["radius_km"]
2.0
>>> rows[0]["capacity"] > 0
True

oroscope.combine_experiments.smallest_radius_for(partner_mask, world, centre_lat, centre_lon, wanted, spacing_km=1.0, grid_type='hex', radii_km=(2.0, 5.0, 10.0, 20.0, 30.0, 40.0, 60.0, 80.0))

The nearest radius at which a partner array of wanted detectors fits.

A convenience over colocation_capacity() for the question actually asked of a candidate site — how far would the partner array have to spread? — rather than the table.

Parameters:

partner_maskndarray

As colocation_capacity().

worldtuple of float

As colocation_capacity().

centre_lat, centre_lonfloat

The site to measure from, in degrees.

wantedint

Detectors the partner array needs.

spacing_kmfloat, optional

Partner detector spacing, in km.

grid_type{‘hex’, ‘square’}, optional

Partner lattice.

radii_kmiterable of float, optional

Radii to try, ascending.

Returns:

float or None

The smallest radius tried at which the capacity reaches wanted, or None if none of them does. None means “not within the largest radius tried”, not “impossible”.

Examples

>>> import numpy as np
>>> from oroscope import combine_experiments as ce
>>> mask = np.ones((400, 400), dtype=bool)
>>> world = (1/3600, 0.0, 0.0, -1/3600, -77.0, -10.0)
>>> ce.smallest_radius_for(mask, world, -10.05, -76.95, 4, radii_km=(2.0, 5.0))
2.0

oroscope.combine_experiments.dem_for_run(run: dict) → str | None

The DEM a run searched, if it can still be found.

The combiner works from the masks the runs wrote, so it never needed the DEM – which is why the overview lost its terrain. The path is recorded twice, in the results parameters and in provenance.json, so it can be recovered when the file is still there and reported as absent when it is not.

Parameters:

rundict

A run, as load_run() returns.

Returns:

str or None

Absolute path to the DEM, or None when it is not recorded or not present.

oroscope.combine_experiments.relief_for_mask(dem_path: str, shape: tuple[int, int]) → ndarray | None

Hillshaded relief for the mask grid, as a greyscale array in [0, 1].

Shaded relief rather than a colour elevation ramp, deliberately: the overlay needs the colour for its categories, and a terrain colourmap underneath would compete with them for the reader’s attention. Relief carries the shape of the ground without spending any hue on it.

The mask is the DEM downsampled, so the DEM is read strided to match rather than resampled. The float32 .npy cache the searcher leaves beside the DEM is memory-mapped when present, which makes a strided read of a 129 Mpx DEM cost almost nothing; otherwise the GeoTIFF is read whole.

Parameters:

dem_pathstr

Path to the DEM the run searched.

shapetuple of int

(rows, cols) of the mask to match.

Returns:

ndarray or None

Relief in [0, 1] with the given shape, or None if the DEM could not be read or does not divide onto the mask grid.

oroscope.combine_experiments.elevation_for_mask(dem_path: str, shape: tuple[int, int]) → ndarray | None

The run’s DEM, read strided onto the mask grid.

The mask is the DEM downsampled, so it is read strided to match rather than resampled. The float32 .npy cache the searcher leaves beside the DEM is memory-mapped when present, which makes a strided read of a 129 Mpx DEM cost almost nothing; otherwise the GeoTIFF is read whole.

Parameters:

dem_pathstr

Path to the DEM the run searched.

shapetuple of int

(rows, cols) of the mask to match.

Returns:

ndarray or None

Elevation in metres with the given shape, or None if the DEM could not be read or is smaller than the mask.

oroscope.combine_experiments.geographic_extent(world: tuple[float, ...], shape: tuple[int, int]) → tuple

The (left, right, bottom, top) an overlay needs to sit in degrees.

Without it imshow puts a map in pixel coordinates, which is why the overview used to carry no axes at all: pixel indices are not worth labelling, so the ticks were switched off rather than made meaningful.

The world file’s terms are the centre of the top-left pixel, so half a pixel is added back at each edge to give the outer bounds of the raster.

Parameters:

worldtuple of float

The six affine terms, from read_world_file().

shapetuple of int

(rows, cols) of the raster.

Returns:

tuple of float

(left, right, bottom, top) in degrees, for imshow(extent=...).

Examples

>>> from oroscope import combine_experiments as ce
>>> world = (0.01, 0.0, 0.0, -0.01, -72.0, -15.0)     # 0.01 deg pixels
>>> left, right, bottom, top = ce.geographic_extent(world, (100, 200))
>>> round(left, 3), round(right, 3)
(-72.005, -70.005)
>>> round(bottom, 3), round(top, 3)
(-15.995, -14.995)

The extent spans exactly rows * cell degrees: 100 rows of 0.01 is 1.0, from -14.995 to -15.995. This example read -16.005 until it was first executed, which would have made the raster 1.01 degrees tall.

One-at-a-time sensitivity sweep over a search’s parameters.

A single site-search result is only as firm as the assumptions behind it, and several of the criteria are choices rather than measurements: what fraction of shower maximum still counts as a usable shower, how steep the far wall must be, which energy stands in for a spectrum in the tau-decay term, where the score cut sits. This runs the pipeline repeatedly, varying one parameter at a time about a baseline, and tabulates how much each one moves the answer.

One-at-a-time rather than a full grid on purpose: the question here is “which assumption is this result most sensitive to”, which OAT answers directly and cheaply. It does not capture interactions between parameters, and a result that hinges on an interaction will not show up – for that, sweep the pair explicitly.

python sensitivity.py ../config/tambo_colca_config.json

–sweep decay_energy_pev 3 10 55 100 1000 –sweep min_score 0.0 0.2 0.35 0.5

oroscope.sensitivity.run_once(config, out_dir, verbose=False, max_memory_gb=None, timeout=3600)

Runs the pipeline once, in a subprocess.

A subprocess rather than a function call, for two reasons that a sweep makes unavoidable. Memory is reclaimed completely between points: running ten searches in one process took 6.9 GB and was killed by the kernel, because matplotlib retains every figure it is not explicitly asked to close and the leak compounded. And one point that fails – an impossible parameter, or a genuine out-of-memory – reports a failed row rather than ending the sweep.

The cost is a few seconds of Numba compilation per point, which is the right trade for a sweep that otherwise cannot finish.

Parameters:

configdict

Parameters as they appear in a config file. Keys beginning with _ are treated as comments and dropped.

out_dirstr

Directory for this run’s outputs.

verbosebool, optional

Let the child’s output through rather than capturing it.

max_memory_gbfloat, optional

Address-space ceiling for the child.

timeoutfloat, optional

Seconds before the child is killed and the point reported as failed.

Returns:

dict or None

The parsed results JSON, or None when the run produced nothing. Read back from disk rather than returned directly: the pipeline hands its results to its own caller, and that caller is in another process here – which is the whole point of running each point in one.

oroscope.sensitivity.summarise(results: dict | None) → dict

Reduces a results JSON to the few numbers a sensitivity table compares.

Parameters:

resultsdict or None

A parsed results JSON, or None when the run produced nothing.

Returns:

dict

sites, capacity, area_km2, accepted, kept and acceptance. All zero for a run that produced nothing, so a sweep row still appears rather than vanishing.

Figures

The schematics used throughout this documentation, as importable functions so they can be restyled and reused.

Publication-quality schematics of the geometry the search computes.

Each routine returns a matplotlib.figure.Figure built from the same physics the code uses, so a diagram cannot drift away from the implementation it illustrates. They are here rather than in the documentation tree so they can be imported, restyled and reused — in a talk, a proposal, or a paper — without copying code out of an .rst file.

Styling is applied per-figure through a context manager rather than by mutating global rcParams, so importing this module does not change the appearance of anybody else’s plots.

Examples

>>> from oroscope import figures
>>> fig = figures.walk_mechanism()
>>> fig.savefig("walk.pdf", bbox_inches="tight")

oroscope.figures.walk_mechanism(earth_radius_km=6371.0, detector_elevation_km=1.05, bin_edges_deg=(-1.2, -0.4, 0.4, 1.2), figsize=(10.2, 4.15))

One profile walk fills every elevation bin at once.

The central algorithmic claim, drawn. Panel (a) is the terrain profile with the first terrain met for each elevation bin; panel (b) is the quantity that makes it work, the apparent elevation angle

\[\theta_{\rm terrain}(d) = \arctan\!\left( \frac{z(d) - d^2/2R - z_0}{d}\right),\]

together with its running maximum. Because that maximum only increases, each new value claims a contiguous band of bins, so one pass fills them all — which is why the elevation binning is nearly free and the azimuth count sets the cost.

Parameters:

earth_radius_kmfloat, optional

Earth radius used for the \(d^2/2R\) curvature drop. The true radius, not the inflated radio one: the particle trajectory is not refracted.

detector_elevation_kmfloat, optional

Elevation of the candidate pixel the walk starts from.

bin_edges_degsequence of float, optional

Elevation-bin edges to trace, in degrees. Drawn in the order given, coloured low to high.

figsizetuple of float, optional

Figure size in inches.

Returns:

matplotlib.figure.Figure

The finished two-panel figure.

Examples

>>> from oroscope import figures
>>> fig = figures.walk_mechanism()
>>> len(fig.axes)
2

oroscope.figures.canyon_geometry(depth_m=1500.0, floor_width_m=1000.0, wall_slope_deg=40.6, figsize=(7.6, 4.0))

The canyon-crossing geometry a particle array such as TAMBO selects on.

Drawn to scale, unlike the long-range figure: a canyon is a few kilometres across and one and a half deep, so no exaggeration is needed and none is applied. That is itself worth showing — the arrival directions really do span tens of degrees here, where GRAND’s span three.

Two criteria are marked because they are separate and are easy to conflate. The near wall is the ground the array stands on and must be deployable; the far wall is where the tau exits and must be steep. A single slope band cannot express both.

Parameters:

depth_mfloat, optional

Rim-to-floor depth. Colca is about 1500 m.

floor_width_mfloat, optional

Width of the flat valley floor.

wall_slope_degfloat, optional

Slope of both walls. Colca’s published depth and ~4.5 km rim separation imply about 40.6 degrees, which is far outside GRAND’s 3-25 degree deployable band.

figsizetuple of float, optional

Figure size in inches.

Returns:

matplotlib.figure.Figure

Examples

>>> from oroscope import figures
>>> fig = figures.canyon_geometry()
>>> round(float(fig.get_figwidth()), 1)
7.6

oroscope.figures.decay_and_shower(energies_pev=(3.0, 10.0, 55.0, 100.0, 1000.0), crossing_m=3000.0, figsize=(7.4, 3.9))

Why a single energy cannot stand in for a spectrum.

The tau must decay inside the gap for a shower to reach the detector, with probability \(1 - \exp(-d/L)\) for a boosted decay length \(L\). Across a canyon that probability runs from essentially one to a few per cent over a single experiment’s energy reach, which is why a capacity computed at one representative energy is an artefact of the energy chosen rather than a property of the terrain.

Parameters:

energies_pevsequence of float, optional

Energies to mark, in PeV.

crossing_mfloat, optional

Gap the tau must decay within, in metres.

figsizetuple of float, optional

Figure size in inches.

Returns:

matplotlib.figure.Figure

Examples

>>> from oroscope import figures
>>> fig = figures.decay_and_shower()
>>> len(fig.axes)
1

oroscope.figures.pipeline_stages(figsize=(9.2, 5.4))

How a DEM becomes a list of sites: the stages, and what each one removes.

The vocabulary this project uses — screening, striding, the arrival scan, scoring, closing, pruning — is introduced nowhere in one place, and the terms are not guessable. This is that place, drawn.

The widths are proportional to the survivors at each stage, on a logarithmic scale because the range is six orders of magnitude and a linear funnel would show one visible bar and six slivers. The numbers are a real run: TAMBO over the full Ancash DEM, 68.6 Mpx.

Read it as two halves. Everything down to the arrival scan removes candidates; everything below rebuilds a map from them, which is why the count rises again at closing. Confusing those two halves is the single commonest way to misread a funnel table.

The arrival scan and the scoring are two bars, which they could not be before. The pre-fix funnel recorded the post-cut count under both names — the defect fixed in run_arrival_scan() — so the two were drawn merged, as a single bar reading 1,022,530, and the arrival window carried the blame for a cut it had not made. Separated on a post-fix run, the scan keeps 82% of what striding hands it and min_score takes 8.4x. That is the reverse of what the merged bar implied, and it is the whole reason the stage is worth drawing on its own.

Parameters:

figsizetuple of float, optional

Figure size in inches.

Returns:

matplotlib.figure.Figure

Examples

>>> from oroscope import figures
>>> fig = figures.pipeline_stages()
>>> len(fig.axes)
1

oroscope.figures.striding_and_closing(stride=5, element_px=(3, 5), figsize=(10.4, 3.2))

Why the closing element has to outrun the gap that striding leaves.

Striding keeps one surviving pixel in stride, so the accepted set becomes a lattice of isolated marks. Morphological closing is what turns that back into a region — but only if its structuring element is larger than the gap. Below the gap the marks never touch and the mask stays a scatter; above it the region reappears almost intact.

The transition is at the gap and it is abrupt, not gradual. A second, far smaller step follows at twice the gap, where the element grows wide enough to bridge second-neighbour marks as well. That is the whole content of the figure, and it is the mechanism behind a real 1.51x under-report of TAMBO’s area at Colca — and a 23.0x one on the steeper ground of the Callejon de Huaylas, where the accepted strips are narrower still. Both were far larger at TAMBO’s old 100 m element, 4.75x and 291x, which is the cliff drawn here: 100 m is three pixels against a five-pixel gap and 150 m is five.

Parameters:

strideint, optional

Keeps every Nth surviving pixel; also the gap it leaves, in pixels.

element_pxtuple of int, optional

Closing element sizes to draw, in pixels. One below the gap, one at it, one above.

figsizetuple of float, optional

Figure size in inches.

Returns:

matplotlib.figure.Figure

Examples

>>> from oroscope import figures
>>> fig = figures.striding_and_closing()
>>> len(fig.axes)
5

oroscope.figures.score_composition(cut=0.35, figsize=(7.8, 3.4))

Why a threshold on a product of components sits on a cliff.

Each component scores a candidate in [0, 1] against one named criterion — depth, accepted solid angle, exit distance, and so on. They are combined by multiplication, so a candidate has to be good at everything, and the composed score of several components piles up near zero however good the terrain is.

The curves are synthetic. Six independent draws from a Beta(5, 2) stand in for six components; no search, terrain or stored result is involved. The distribution is deliberately generous — mean 0.714, mode 0.80 — so that the collapse cannot be blamed on poor terrain. This figure demonstrates the mechanism; it is not a measurement, and the page that publishes it prints a measured sentence directly beneath, which is exactly the confusion the note on the axes exists to prevent.

A cut placed in the middle of that pile is therefore not a mild preference: it is a cliff, and where it lands depends on how many components happen to be enabled. Adding a component moves every score down and so silently tightens the cut.

Parameters:

cutfloat, optional

Where min_score is placed, for illustration.

figsizetuple of float, optional

Figure size in inches.

Returns:

matplotlib.figure.Figure

Examples

>>> from oroscope import figures
>>> fig = figures.score_composition()
>>> len(fig.axes)
2

oroscope.figures.score_composition_measured(cut=0.35, figsize=(7.8, 3.4))

What the product does to real candidates, rather than to a model of them.

The companion to score_composition(), which demonstrates the mechanism with invented components. This one is measured: each curve is the running product over the same 360,939 geometrically accepted candidates of one Colca TAMBO run, with one more component multiplied in. No candidate is removed between curves — the population is fixed and only its scores move left. What survives the cut at the end is 64,152 of them.

The result does not resemble the mechanism figure, and that is the point. depth and distance are exactly 1.0 for every candidate and shower for 92.3% of them, so three of the six components do nothing whatever. Five of the six together still leave 96.3% above the shipped cut. The sixth, solid_angle, takes it to 17.8%.

So min_score is not really a threshold on a product of six criteria. It is a cut on solid_angle wearing a product as a disguise — the measured form of the caveat that solid_angle is the weakest component at every selected site in every region.

Components enter least-restrictive first, so the collapse is attributable to a named component rather than to the number of them.

Parameters:

cutfloat, optional

Where min_score is drawn. The stored percentages are those of the shipped 0.35 and are not recomputed, so another value moves the line without moving the curve beside it.

figsizetuple of float, optional

Figure size in inches.

Returns:

matplotlib.figure.Figure

Examples

>>> from oroscope import figures
>>> fig = figures.score_composition_measured()
>>> len(fig.axes)
2

oroscope.figures.grand_and_tambo_scales(xmax_km=50.0, figsize=(10.4, 5.8))

The same question, asked ten kilometres apart and four hundred metres apart.

canyon_geometry() draws what TAMBO asks of the ground. This draws what both experiments ask, on one shared horizontal axis, because the comparison is the point and a reader should not have to do it in their head.

The upper row is GRAND: antennas on ground inside its 3–25° band, watching a massif through a window 3° about the horizon. The lower row is the Colca cross-section at the same scale — the notch in the first 4.5 km — with a magnified detail tied back to its true position, so it cannot be mistaken for a second canyon further out.

The window stops where the highest ray first meets ground, which is the acceptance test itself rather than a decoration: shading past the massif would claim sight lines the massif blocks.

One asymmetry is deliberate. The lower edge of GRAND’s window is invisible because the ground beyond the detector is level, so a ray 3° below the horizon meets it at once. That is the honest picture of a detector on a plain.

Parameters:

xmax_kmfloat, optional

Ground distance spanned by both rows, in km. The comparison only works while both rows share it.

figsizetuple of float, optional

Figure size in inches.

Returns:

matplotlib.figure.Figure

Examples

>>> from oroscope import figures
>>> fig = figures.grand_and_tambo_scales()
>>> len(fig.axes)
2

Fetching data

The one-shot tools that bring a region’s inputs onto disk are setup rather than library, which is why import oroscope does not re-export them and why their command-line side is documented on The command line instead. Only the names other pages link into are listed here: oroscope.fetch_dem.REGIONS is cross-referenced from Getting the data, and a cross-reference with no target renders as plain text rather than failing the build, so it sat there unresolved and unnoticed.

Downloads the elevation models a search needs, and writes configurations for them.

Fetches the bundled regions from OpenTopography into input/dem/ and generates a ready-to-run JSON config for each in config/:

oroscope-fetch-dem --region arequipa --open_topography_api_key KEY

Four regions are defined. Three departments at 1 arc-second from the same dataset, so that runs over them are comparable – arequipa (129 Mpx), lima (105 Mpx) and ancash (69 Mpx), all SRTMGL1 – and peru, the whole country, at 3 arc-seconds (SRTMGL3) because that is the only resolution that fits either a desktop’s memory or the API’s own area limit. Omit --region to fetch all of them.

Getting a key. It is free and takes a minute. Register at https://portal.opentopography.org/myopentopo, sign in, then open myOpenTopo Authorizations and API Key from the account menu and copy the key. Pass it as --open_topography_api_key, or set OPENTOPOGRAPHY_API_KEY in the environment, which keeps it out of your shell history and out of any file that might be committed.

Requests are capped by area, and the cap is per dataset: 450,000 km² for every 30 m dataset, 4,050,000 km² for the 90 m ones. That is why peru is SRTMGL3 – its bounding box is about 2.86 million km², six times over the 30 m limit.

For any other region, download the tiles from the OpenTopography portal, merge them into one GeoTIFF if the area spans several, and cut the window you want with crop_dem.

This was setup.py, whose name made pip install run the downloader instead of building the package.

oroscope.fetch_dem.REGIONS = {'ancash': {'demtype': 'SRTMGL1', 'east': -76.7257441, 'filename': 'ancash_SRTMGL1.tif', 'north': -8.049709, 'preset': 'default', 'south': -10.7873076, 'west': -78.6584805}, 'arequipa': {'demtype': 'SRTMGL1', 'east': -70.0852632522583, 'filename': 'arequipa_SRTMGL1.tif', 'north': -14.555380967667489, 'preset': 'arequipa', 'south': -17.38995824658555, 'west': -73.58612537384033}, 'lima': {'demtype': 'SRTMGL1', 'east': -75.39955615997313, 'filename': 'lima_SRTMGL1.tif', 'north': -10.228479499469358, 'preset': 'lima', 'south': -13.252477566131276, 'west': -78.07665824890137}, 'peru': {'demtype': 'SRTMGL3', 'east': -68.6, 'filename': 'peru_SRTMGL3.tif', 'north': 0.0, 'preset': 'default', 'south': -18.4, 'west': -81.4}}

The regions oroscope-fetch-dem --region knows, each mapping to its dataset and its bounding box in degrees. A #: comment rather than a plain one so that autodoc picks it up: Getting the data cross-references this name, and module-level data with only an ordinary comment above it is documented nowhere, leaving the reference to render as unlinked text.

Bounds are given as (West, East, South, North).

oroscope.fetch_dem.download_dem(region_name, bounds, api_key, output_dir)

Fetches one region’s GeoTIFF from OpenTopography.

Parameters:

region_namestr

A key of REGIONS, naming the box and dataset to fetch.

boundsdict

That region’s entry: demtype and the south/north/west/east bounds in degrees, plus the filename to write.

api_keystr

An OpenTopography key. Free, from https://portal.opentopography.org/myopentopo.

output_dirstr

Directory to write the .tif into. Created if absent.

Returns:

str or None

The path written, or None if the request failed.

oroscope.fetch_dem.generate_and_patch_config(region_name, preset, dem_filepath, config_dir=None)

Writes a default config for the region and points it at the DEM just downloaded.

This used to shell out to site_searcher.py in the current directory, which was right when the modules were flat and is wrong now that they are a package: an installed copy has no such file anywhere, so the config step failed for every user who was not standing in src/. site_searcher.generate_config() is the same code path the --generate_config flag takes.

Parameters:

region_namestr

Key in REGIONS; names the config file.

presetstr

Config preset to seed from, one of site_searcher.CONFIG_PRESETS.

dem_filepathstr

The DEM to point dem_path at.

config_dirstr, optional

Where to write. Defaults to ../config for continuity with the old layout.