diff --git a/src/cli.rs b/src/cli.rs index 4855011..125338a 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -724,6 +724,8 @@ impl<'a> BirliContext<'a> { .help_heading("FLAGGING"), // corrections + arg!(--"van-vleck" "Apply Van Vleck corrections") + .help_heading("CORRECTION"), arg!(--"no-cable-delay" "Do not perform cable length corrections") .help_heading("CORRECTION"), arg!(--"no-geometric-delay" "Do not perform geometric corrections") @@ -1291,6 +1293,20 @@ impl<'a> BirliContext<'a> { (_, true) => RADec::from_mwalib_tile_pointing(&corr_ctx.metafits_context), _ => RADec::from_mwalib_phase_or_pointing(&corr_ctx.metafits_context), }; + prep_ctx.correct_van_vleck = match ( + matches.is_present("van-vleck"), + corr_ctx.mwa_version + ) { + (true, MWAVersion::CorrLegacy) => true, + (true, _) => { + return Err(BirliError::CLIError(InvalidCommandLineArgument { + option: "--van-vleck".into(), + expected: "only available for legacy".into(), + received: "van-vleck".into(), + })) + } + _ => false, + }; prep_ctx.correct_cable_lengths = { let cable_delays_disabled = matches.is_present("no-cable-delay"); let cable_delays_applied = corr_ctx.metafits_context.cable_delays_applied; diff --git a/src/error.rs b/src/error.rs index e718813..a545f69 100644 --- a/src/error.rs +++ b/src/error.rs @@ -3,7 +3,10 @@ use marlu::{io::error::BadArrayShape, mwalib}; use thiserror::Error; -use crate::corrections::{DigitalGainCorrection, PassbandCorrection}; +use crate::{ + corrections::{DigitalGainCorrection, PassbandCorrection}, + van_vleck::VanVleckCorrection, +}; /// Errors relating to CI #[cfg(feature = "cli")] @@ -66,6 +69,10 @@ pub enum BirliError { /// Error derived from [`crate::corrections::DigitalGainCorrection`] DigitalGainCorrection(#[from] DigitalGainCorrection), + #[error(transparent)] + /// Error derived from [`crate::van_vleck::VanVleckCorrection`] + VanVleckCorrection(#[from] VanVleckCorrection), + #[error("You selected dry run")] /// enum variant for when a dry run is selected DryRun {}, diff --git a/src/lib.rs b/src/lib.rs index a2cc932..da32ee0 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -193,6 +193,7 @@ pub mod preprocessing; pub use preprocessing::PreprocessContext; pub mod van_vleck; +pub use van_vleck::correct_van_vleck; cfg_if! { if #[cfg(feature = "cli")] { diff --git a/src/preprocessing.rs b/src/preprocessing.rs index afcdad2..db4763f 100644 --- a/src/preprocessing.rs +++ b/src/preprocessing.rs @@ -4,11 +4,12 @@ use crate::{ correct_cable_lengths, correct_geometry, corrections::{correct_coarse_passband_gains, correct_digital_gains, ScrunchType}, marlu::{mwalib::CorrelatorContext, ndarray::prelude::*, Jones, LatLngHeight, RADec}, + van_vleck::correct_van_vleck, with_increment_duration, BirliError, VisSelection, }; use cfg_if::cfg_if; use derive_builder::Builder; -use log::trace; +use log::{trace, warn}; use std::{ fmt::{Debug, Display}, time::Duration, @@ -31,6 +32,9 @@ pub struct PreprocessContext<'a> { /// The phase centre used for geometric corrections pub phase_centre: RADec, + /// Whether Van Vleck corrections are enabled + #[builder(default = "false")] + pub correct_van_vleck: bool, /// Whether cable length corrections are enabled #[builder(default = "true")] pub correct_cable_lengths: bool, @@ -57,6 +61,15 @@ pub struct PreprocessContext<'a> { impl Display for PreprocessContext<'_> { fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + writeln!( + f, + "{} correct Van Vleck.", + if self.correct_van_vleck { + "Will" + } else { + "Will not" + } + )?; writeln!( f, "{} correct cable lengths.", @@ -110,6 +123,11 @@ impl<'a> PreprocessContext<'a> { /// A one line description of the tasks preprocessing will do. pub fn as_comment(&self) -> String { [ + if self.correct_van_vleck { + Some("Van Vleck corrections".to_string()) + } else { + None + }, if self.correct_cable_lengths { Some("cable length corrections".to_string()) } else { @@ -162,6 +180,17 @@ impl<'a> PreprocessContext<'a> { mut flag_array: ArrayViewMut3, vis_sel: &VisSelection, ) -> Result<(), BirliError> { + let sel_ant_pairs = vis_sel.get_ant_pairs(&corr_ctx.metafits_context); + + if self.correct_van_vleck { + warn!("Van Vleck correction is a work in progress!"); + trace!("correcting van vleck"); + with_increment_duration!( + "correct_van_vleck", + correct_van_vleck(corr_ctx, jones_array.view_mut(), &sel_ant_pairs)? + ); + } + if self.correct_cable_lengths { trace!("correcting cable lengths"); with_increment_duration!( @@ -176,8 +205,6 @@ impl<'a> PreprocessContext<'a> { ); } - let sel_ant_pairs = vis_sel.get_ant_pairs(&corr_ctx.metafits_context); - if self.correct_digital_gains { trace!("correcting digital gains"); with_increment_duration!( diff --git a/src/van_vleck.rs b/src/van_vleck.rs index eb69a5e..f2db675 100644 --- a/src/van_vleck.rs +++ b/src/van_vleck.rs @@ -65,7 +65,7 @@ pub enum VanVleckCorrection { /// # Examples /// /// ```rust -/// use birli::{correct_cable_lengths, mwalib::CorrelatorContext, VisSelection, io::read_mwalib}; +/// use birli::{correct_van_vleck, mwalib::CorrelatorContext, VisSelection, io::read_mwalib}; /// /// // define our input files /// let metafits_path = "tests/data/1297526432_mwax/1297526432.metafits"; @@ -90,7 +90,7 @@ pub enum VanVleckCorrection { /// /// let ant_pairs = vis_sel.get_ant_pairs(&corr_ctx.metafits_context); /// -/// van_vleck_correction(&corr_ctx, jones_array.view_mut(), &ant_pairs); +/// correct_van_vleck(&corr_ctx, jones_array.view_mut(), &ant_pairs); /// ``` /// /// # Errors @@ -222,7 +222,7 @@ pub fn correct_van_vleck( ] }) .collect_vec(); - dbg!(&sighat_xxr_yyr); + // dbg!(&sighat_xxr_yyr); // partitioning of auto jones matrix into vectors of [time * freq * n_autos]: // - `sighat_xxr`, `sighat_yyr` <- sqrt of xx real and yy real, for correction with with van_vleck_autos @@ -292,16 +292,16 @@ pub fn correct_van_vleck( // map sigma_xxr_yyr pol index 0,1 to jones index 0,3 (xx, yy) by multiplying by 3 let jp_idx = p_idx * 3; let re = (nsamples as f64) * s.powi(2); - println!( - "j@({:?},{:?}).tf[({:?},{:?})][{:?}]={:?}<-{:?}", - ant1, - ant2, - t_idx, - f_idx, - jp_idx, - j_tf[(t_idx, f_idx)][jp_idx].re, - re - ); + // println!( + // "j@({:?},{:?}).tf[({:?},{:?})][{:?}]={:?}<-{:?}", + // ant1, + // ant2, + // t_idx, + // f_idx, + // jp_idx, + // j_tf[(t_idx, f_idx)][jp_idx].re, + // re + // ); j_tf[(t_idx, f_idx)][jp_idx].re = re as f32; // j_tf[(t_idx, f_idx)][p_idx * 3].im = 0.0; ? }); @@ -360,7 +360,7 @@ pub fn van_vleck_autos(hat: &[f64]) -> Vec { hat.par_iter() .map(|&sighat| { van_vleck_auto(sighat).map_or(sighat, |sigma| { - println!("sigma={sigma} <- sighat={sighat}"); + // println!("sigma={sigma} <- sighat={sighat}"); sigma }) })