Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 10 additions & 1 deletion src/cli/common/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ use crate::{
AverageFactorError,
},
beam::Beam,
constants::{DEFAULT_VETO_THRESHOLD, MWA_HEIGHT_M, MWA_LAT_DEG, MWA_LONG_DEG},
constants::{DEFAULT_ELEVATION_LIMIT, DEFAULT_VETO_THRESHOLD, MWA_HEIGHT_M, MWA_LAT_DEG, MWA_LONG_DEG},
io::{
get_single_match_from_glob,
write::{can_write_to_file, VisOutputType, VIS_OUTPUT_EXTENSIONS},
Expand Down Expand Up @@ -340,6 +340,12 @@ pub(super) struct SkyModelWithVetoArgs {

#[clap(long, help = VETO_THRESHOLD_HELP.as_str(), help_heading = "SKY MODEL")]
pub(super) veto_threshold: Option<f64>,

/// Minimum elevation for a source to be included in the sky model [degrees].
/// Sources with any component below this elevation are discarded. Default: 0.
#[clap(long, help_heading = "SKY MODEL")]
pub(super) elevation_limit: Option<f64>,

/// Optional source names to include (or exclude if --invert).
/// Skips vetoing if provided, unless --invert is enabled.
#[clap(long, multiple_values(true), help_heading = "SKY MODEL SOURCES")]
Expand All @@ -364,6 +370,7 @@ impl SkyModelWithVetoArgs {
num_sources: self.num_sources.or(other.num_sources),
source_dist_cutoff: self.source_dist_cutoff.or(other.source_dist_cutoff),
veto_threshold: self.veto_threshold.or(other.veto_threshold),
elevation_limit: self.elevation_limit.or(other.elevation_limit),
named_sources: self.named_sources.or(other.named_sources),
invert: self.invert.or(other.invert),
}
Expand All @@ -387,6 +394,7 @@ impl SkyModelWithVetoArgs {
num_sources,
source_dist_cutoff,
veto_threshold,
elevation_limit,
named_sources,
invert,
} = self;
Expand Down Expand Up @@ -506,6 +514,7 @@ impl SkyModelWithVetoArgs {
num_sources,
source_dist_cutoff.unwrap_or(f64::MAX),
veto_threshold.unwrap_or(DEFAULT_VETO_THRESHOLD),
elevation_limit.unwrap_or(DEFAULT_ELEVATION_LIMIT),
)?;
if sl.is_empty() {
return Err(ReadSourceListError::NoSourcesAfterVeto);
Expand Down
10 changes: 9 additions & 1 deletion src/cli/srclist/by_beam/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ use crate::{
display_warnings, BeamArgs, Warn, ARRAY_POSITION_HELP, SOURCE_LIST_INPUT_TYPE_HELP,
SOURCE_LIST_OUTPUT_TYPE_HELP, VETO_THRESHOLD_HELP,
},
constants::DEFAULT_VETO_THRESHOLD,
constants::{DEFAULT_ELEVATION_LIMIT, DEFAULT_VETO_THRESHOLD},
metafits::get_dipole_delays,
srclist::{
read::read_source_list_file, veto_sources, write_source_list, ReadSourceListError,
Expand Down Expand Up @@ -121,6 +121,11 @@ pub struct SrclistByBeamArgs {
#[clap(long, help = VETO_THRESHOLD_HELP.as_str(), help_heading = "SOURCE FILTERING")]
veto_threshold: Option<f64>,

/// Minimum elevation for a source to be included in the sky model [degrees].
/// Sources with any component below this elevation are discarded. Default: 0.
#[clap(long, help_heading = "SOURCE FILTERING")]
elevation_limit: Option<f64>,

/// Don't include point components from the input sky model.
#[clap(long, help_heading = "SOURCE FILTERING")]
filter_points: bool,
Expand Down Expand Up @@ -171,6 +176,7 @@ impl SrclistByBeamArgs {
self.freqs_hz.as_deref(),
self.source_dist_cutoff,
self.veto_threshold,
self.elevation_limit,
self.filter_points,
self.filter_gaussians,
self.filter_shapelets,
Expand Down Expand Up @@ -204,6 +210,7 @@ fn by_beam(
freqs_hz: Option<&[f64]>,
source_dist_cutoff: Option<f64>,
veto_threshold: Option<f64>,
elevation_limit: Option<f64>,
filter_points: bool,
filter_gaussians: bool,
filter_shapelets: bool,
Expand Down Expand Up @@ -359,6 +366,7 @@ fn by_beam(
None,
source_dist_cutoff.unwrap_or(f64::MAX),
veto_threshold.unwrap_or(DEFAULT_VETO_THRESHOLD),
elevation_limit.unwrap_or(DEFAULT_ELEVATION_LIMIT),
)?;
// Were any sources left after vetoing?
if sl.is_empty() {
Expand Down
1 change: 1 addition & 0 deletions src/cli/srclist/by_beam/tests.rs
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ fn test_srclist_by_beam() {
number: n,
source_dist_cutoff: None,
veto_threshold: None,
elevation_limit: None,
filter_points: false,
filter_gaussians: false,
filter_shapelets: false,
Expand Down
3 changes: 2 additions & 1 deletion src/cli/vis_utils/simulate/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ use vec1::Vec1;
use crate::{
averaging::{parse_freq_average_factor, parse_time_average_factor, timesteps_to_timeblocks},
beam::{create_fee_beam_object, create_no_beam_object, Beam, Delays},
constants::{DEFAULT_CUTOFF_DISTANCE, DEFAULT_VETO_THRESHOLD},
constants::{DEFAULT_CUTOFF_DISTANCE, DEFAULT_ELEVATION_LIMIT, DEFAULT_VETO_THRESHOLD},
context::Polarisations,
help_texts::{
ARRAY_POSITION_HELP, DIPOLE_DELAYS_HELP, SOURCE_DIST_CUTOFF_HELP, VETO_THRESHOLD_HELP,
Expand Down Expand Up @@ -651,6 +651,7 @@ impl VisSimParams {
*num_sources,
source_dist_cutoff.unwrap_or(DEFAULT_CUTOFF_DISTANCE),
veto_threshold.unwrap_or(DEFAULT_VETO_THRESHOLD),
DEFAULT_ELEVATION_LIMIT,
)?;
if source_list.is_empty() {
return Err(VisSimulateError::NoSourcesAfterVeto);
Expand Down
2 changes: 1 addition & 1 deletion src/constants.rs
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ pub(crate) const DEFAULT_VETO_THRESHOLD: f64 = 0.01;

/// Sources with elevations less than this value are discarded from sky-model
/// source lists \[degrees\].
pub(crate) const ELEVATION_LIMIT: f64 = 0.0;
pub(crate) const DEFAULT_ELEVATION_LIMIT: f64 = 0.0;

// sqrt(pi^2 / (2 ln(2)))
pub(crate) const SQRT_FRAC_PI_SQ_2_LN_2: f64 = 2.6682231283184983;
Expand Down
127 changes: 124 additions & 3 deletions src/srclist/veto.rs
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@ use rayon::{iter::Either, prelude::*};

use crate::{
beam::Beam,
constants::*,
srclist::{FluxDensity, ReadSourceListError, SourceList},
};

Expand Down Expand Up @@ -51,6 +50,7 @@ pub(crate) fn veto_sources(
num_sources: Option<usize>,
source_dist_cutoff_deg: f64,
veto_threshold: f64,
min_elevation_deg: f64,
) -> Result<(), ReadSourceListError> {
let dist_cutoff = source_dist_cutoff_deg.to_radians();

Expand All @@ -70,9 +70,9 @@ pub(crate) fn veto_sources(
let mut azels = vec![];
for comp in source.components.iter() {
let azel = comp.radec.to_hadec(lst_rad).to_azel(array_latitude_rad);
if azel.el.to_degrees() < ELEVATION_LIMIT {
if azel.el.to_degrees() < min_elevation_deg {
if log_enabled!(Trace) {
trace!("A component's elevation ({}°, source {source_name}) was below the limit ({ELEVATION_LIMIT}°)", azel.el.to_degrees());
trace!("A component's elevation ({}°, source {source_name}) was below the limit ({min_elevation_deg}°)", azel.el.to_degrees());
}
return Either::Left(Ok(source_name));
}
Expand Down Expand Up @@ -216,6 +216,7 @@ mod tests {
use super::*;
use crate::{
beam::{Delays, FEEBeam, NoBeam},
constants::DEFAULT_ELEVATION_LIMIT,
srclist::{
read::read_source_list_file, ComponentType, FluxDensityType, Source, SourceComponent,
},
Expand Down Expand Up @@ -362,6 +363,7 @@ mod tests {
None,
180.0,
0.1,
DEFAULT_ELEVATION_LIMIT,
);
assert!(result.is_ok());
result.unwrap();
Expand Down Expand Up @@ -407,6 +409,7 @@ mod tests {
Some(3),
180.0,
0.1,
DEFAULT_ELEVATION_LIMIT,
);
assert!(result.is_ok(), "{:?}", result.unwrap_err());
result.unwrap();
Expand Down Expand Up @@ -434,6 +437,7 @@ mod tests {
None,
180.0,
0.1,
DEFAULT_ELEVATION_LIMIT,
);
assert!(result.is_ok(), "{:?}", result.unwrap_err());
result.unwrap();
Expand All @@ -442,4 +446,121 @@ mod tests {
assert_eq!(source_list.get_index(0).unwrap().0, "J004616-420739");
assert_eq!(source_list.get_index(99).unwrap().0, "J000217-253912");
}

/// Sources whose components are below `elevation_limit` must be rejected,
/// and a higher elevation limit must reject more sources than a lower one.
#[test]
fn elevation_limit_rejects_sources() {
let beam = NoBeam { num_tiles: 1 };

// Build a tiny source list with three point sources:
// "above" – always well above the horizon (zenith)
// "horizon" – sitting right at the horizon (el ≈ 0°)
// "below" – well below the horizon
//
// With LST = 0 and the MWA latitude, a source at RA = LST and
// Dec = latitude is at the zenith, while Dec = -(90° - lat) puts it
// on the celestial equator and Dec < -(90° + lat) takes it below the
// horizon.
let lat_deg = MWA_LAT_RAD.to_degrees(); // ≈ -26.7°

let make_point = |ra_deg: f64, dec_deg: f64| Source {
components: vec![SourceComponent {
radec: RADec::from_degrees(ra_deg, dec_deg),
comp_type: ComponentType::Point,
flux_type: FluxDensityType::PowerLaw {
si: -0.8,
fd: FluxDensity {
freq: 180e6,
i: 100.0,
q: 0.0,
u: 0.0,
v: 0.0,
},
},
}]
.into_boxed_slice(),
};

// At LST=0:
// zenith source → dec = lat_deg, RA = 0 → el = 90°
// below-horizon → dec = -(90.0 - lat_deg).abs() - 10.0 → el < 0°
let zenith_dec = lat_deg;
let below_dec = -(90.0 + lat_deg.abs()) - 10.0; // deeply below the horizon

let mut sl = SourceList::new();
sl.insert("zenith".into(), make_point(0.0, zenith_dec));
sl.insert("below_horizon".into(), make_point(0.0, below_dec));

// With DEFAULT_ELEVATION_LIMIT (0°), only the below-horizon source is vetoed.
let mut sl_default = sl.clone();
veto_sources(
&mut sl_default,
RADec::from_degrees(0.0, zenith_dec),
0.0,
MWA_LAT_RAD,
&[180e6],
&beam,
None,
f64::MAX,
0.0,
DEFAULT_ELEVATION_LIMIT,
)
.unwrap();
assert!(
sl_default.contains_key("zenith"),
"zenith source should survive the default elevation limit"
);
assert!(
!sl_default.contains_key("below_horizon"),
"below-horizon source should be rejected by the default elevation limit"
);

// With a high elevation limit (45°), even the zenith source survives
// only if its elevation exceeds 45°. The below-horizon source must
// still be rejected.
let mut sl_high = sl.clone();
veto_sources(
&mut sl_high,
RADec::from_degrees(0.0, zenith_dec),
0.0,
MWA_LAT_RAD,
&[180e6],
&beam,
None,
f64::MAX,
0.0,
45.0,
)
.unwrap();
assert!(
sl_high.contains_key("zenith"),
"zenith source (el=90°) should survive a 45° elevation limit"
);
assert!(
!sl_high.contains_key("below_horizon"),
"below-horizon source should be rejected by a 45° elevation limit"
);

// With an elevation limit above 90°, every source is rejected and
// veto_sources should return an error (no sources remain).
let mut sl_all = sl;
let result = veto_sources(
&mut sl_all,
RADec::from_degrees(0.0, zenith_dec),
0.0,
MWA_LAT_RAD,
&[180e6],
&beam,
None,
f64::MAX,
0.0,
91.0,
);
assert!(result.is_ok(), "veto_sources itself should not error");
assert!(
sl_all.is_empty(),
"all sources should be rejected when elevation limit > 90°"
);
}
}
Loading