Allow Aman's dipole amps to be ingested.

Aman Chokshi is working with adjusting MWA dipole gains/amps, and this commit
uses these gains if available from a metafits file ("DipAmps" column in the
second HDU's table).

Author: cjordan

CHANGELOG.md

@@ -8,6 +8,8 @@ Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [0.3.0] - Unreleased
 ### Added
+- Support for the "DipAmps" column in a metafits file. This allows users to
+  control dipole gains in beam code.
 - Support for averaging incoming visibilities in time and frequency *before*
   doing any work on them.
 - When writing out visibilities, it is now possible to write out the smallest

mdbook/src/defs/mwa/dead_dipoles.md

@@ -10,8 +10,8 @@ account for dead dipoles if possible.
 Beam responses are generated with
 [`hyperbeam`](https://github.com/MWATelescope/mwa_hyperbeam) and dead dipole
 information is encoded as a "dipole gain" of 1 ("alive") or 0 ("dead"). It is
-possible to supply other values for dipole gains, although at the time of
-writing `hyperdrive` only uses ones or zeros.
+possible to supply other values for dipole gains with a "DipAmps" column; see
+the [metafits page](metafits.md).
 
 For the relevant functions, dead dipole information can be ignored by supplying
 a flag `--unity-dipole-gains`. This sets all dipole gains to 1.

mdbook/src/defs/mwa/metafits.md

@@ -29,3 +29,12 @@ that is required to be present by
 [`mwalib`](https://github.com/MWATelescope/mwalib), which is used by
 `hyperdrive` when reading metafits files.
 ~~~
+
+## Controlling dipole gains
+
+If the "TILEDATA" HDU of a metafits contains a "DipAmps" column, each row
+containing 16 double-precision values for bowties in the M&C order, these are
+used as the dipole gains in beam calculations. If the "DipAmps" column isn't
+available, the default behaviour is to use gains of 1.0 for all dipoles, except
+those that have delays of 32 in the "Delays" column (they will have a gain of
+0.0, and are considered [dead](dead_dipoles.md)).

src/cli/common/beam/mod.rs

@@ -177,18 +177,31 @@ impl BeamArgs {
                     dipole_gains
                 };
                 if let Some(dipole_gains) = dipole_gains.as_ref() {
-                    let num_tiles_with_dead_dipoles = dipole_gains
-                        .outer_iter()
-                        .filter(|tile_dipole_gains| {
-                            tile_dipole_gains.iter().any(|g| g.abs() < f64::EPSILON)
-                        })
-                        .count();
-                    printer.push_line(
-                        format!(
-                            "Using dead dipole information ({num_tiles_with_dead_dipoles} tiles affected)"
-                        )
-                        .into(),
-                    );
+                    trace!("Attempting to use dipole gains:");
+                    for row in dipole_gains.outer_iter() {
+                        trace!("{row}");
+                    }
+
+                    // Currently, the only way to have dipole gains other than
+                    // zero or one is by using Aman's "DipAmps" metafits column.
+                    if dipole_gains.iter().any(|&g| g != 0.0 && g != 1.0) {
+                        printer.push_line(
+                            "Using Aman's 'DipAmps' dipole gains from the metafits".into(),
+                        );
+                    } else {
+                        let num_tiles_with_dead_dipoles = dipole_gains
+                            .outer_iter()
+                            .filter(|tile_dipole_gains| {
+                                tile_dipole_gains.iter().any(|g| g.abs() < f64::EPSILON)
+                            })
+                            .count();
+                        printer.push_line(
+                            format!(
+                                "Using dead dipole information ({num_tiles_with_dead_dipoles} tiles affected)"
+                            )
+                            .into(),
+                        );
+                    }
                 } else {
                     // If we don't have dipole gains, we must assume all dipoles
                     // are "alive". But, if any dipole delays are 32, then the

src/cli/common/beam/tests.rs

@@ -2,6 +2,7 @@
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
+use approx::assert_abs_diff_eq;
 use ndarray::array;
 
 use super::BeamArgs;
@@ -91,3 +92,23 @@ fn test_unity_dipole_gains() {
     // Verify that there are no dead dipoles in the delays.
     assert!(beam.get_dipole_delays().unwrap().iter().all(|d| *d != 32));
 }
+
+#[test]
+fn test_aman_dipole_gains() {
+    let f = |metafits| {
+        let metafits = mwalib::MetafitsContext::new(metafits, None).unwrap();
+        let delays = crate::metafits::get_dipole_delays(&metafits);
+        let gains = crate::metafits::get_dipole_gains(&metafits);
+        (delays, gains)
+    };
+
+    let (vanilla_delays, vanilla_gains) = f("test_files/1120082744/1120082744.metafits");
+    let (dipamps_delays, dipamps_gains) = f("test_files/1120082744/1120082744_DipAmps.metafits");
+    assert_eq!(vanilla_delays, dipamps_delays);
+    assert_ne!(vanilla_gains, dipamps_gains);
+
+    // First X dipole for Tile011
+    assert_abs_diff_eq!(dipamps_gains[(0, 0)] as f32, 0.89985347);
+    // First Y dipole for Tile011
+    assert_abs_diff_eq!(dipamps_gains[(0, 16)] as f32, 0.8930142);
+}

src/io/read/fits/error.rs

@@ -75,4 +75,11 @@ pub enum FitsError {
         source_line: u32,
         source_column: u32,
     },
+
+    #[error("When attempting to read column {col_name} from HDU {hdu_num}, cfitsio gave an error: {err}")]
+    ReadCellArray {
+        col_name: &'static str,
+        hdu_num: usize,
+        err: Box<fitsio::errors::Error>,
+    },
 }

src/io/read/fits/mod.rs

@@ -11,6 +11,7 @@ pub(crate) use error::FitsError;
 use std::{
     ffi::{CStr, CString},
     fmt::Display,
+    os::raw::c_char,
     ptr,
 };
 
@@ -329,3 +330,57 @@ pub(crate) fn _fits_get_float_image_into_buffer(
 
     Ok(())
 }
+
+/// Pull out fits array-in-a-cell values; useful for e.g. STABXYZ. This function
+/// must be pointed at data containing f64s (i.e. TDOUBLE).
+#[track_caller]
+pub(crate) fn fits_read_cell_f64_array<const NUM_ELEM: usize>(
+    fits_ptr: &mut fitsio::FitsFile,
+    hdu: &fitsio::hdu::FitsHdu,
+    col_name: &'static str,
+    row: i64,
+) -> Result<[f64; NUM_ELEM], FitsError> {
+    unsafe {
+        // With the column name, get the column number.
+        let mut status = 0;
+        let mut col_num = -1;
+        let keyword = std::ffi::CString::new(col_name).expect("CString::new failed");
+        // ffgcno = fits_get_colnum
+        fitsio_sys::ffgcno(
+            fits_ptr.as_raw(),
+            0,
+            keyword.as_ptr() as *mut c_char,
+            &mut col_num,
+            &mut status,
+        );
+        // Check the status.
+        fitsio::errors::check_status(status).map_err(|err| FitsError::ReadCellArray {
+            col_name,
+            hdu_num: hdu.number + 1,
+            err: Box::new(err),
+        })?;
+
+        // Now get the specified row from that column.
+        let mut array = [0.0; NUM_ELEM];
+        // ffgcv = fits_read_col
+        fitsio_sys::ffgcv(
+            fits_ptr.as_raw(),
+            82, // TDOUBLE (fitsio.h)
+            col_num,
+            row + 1,
+            1,
+            NUM_ELEM.try_into().expect("not larger than i64::MAX"),
+            std::ptr::null_mut(),
+            array.as_mut_ptr().cast(),
+            &mut 0,
+            &mut status,
+        );
+        fitsio::errors::check_status(status).map_err(|err| FitsError::ReadCellArray {
+            col_name,
+            hdu_num: hdu.number + 1,
+            err: Box::new(err),
+        })?;
+
+        Ok(array)
+    }
+}

src/io/read/uvfits/error.rs

@@ -95,13 +95,6 @@ pub enum UvfitsReadError {
     #[error("When attempting to read uvfits baseline metadata, cfitsio gave an error: {0}")]
     Metadata(fitsio::errors::Error),
 
-    #[error("When attempting to read uvfits column {col_name} from HDU {hdu_num}, cfitsio gave an error: {err}")]
-    ReadCellArray {
-        col_name: &'static str,
-        hdu_num: usize,
-        err: fitsio::errors::Error,
-    },
-
     #[error("When attempting to read uvfits row {row_num}, cfitsio gave an error: {err}")]
     ReadVis {
         row_num: usize,

src/io/read/uvfits/mod.rs

@@ -17,7 +17,6 @@ use std::{
     borrow::Cow,
     collections::{HashMap, HashSet},
     num::NonZeroU16,
-    os::raw::c_char,
     path::{Path, PathBuf},
 };
 
@@ -40,6 +39,7 @@ use crate::{
     io::read::{
         fits::{
             fits_get_col, fits_get_optional_key, fits_get_required_key, fits_open, fits_open_hdu,
+            fits_read_cell_f64_array,
         },
         VisRead, VisReadError,
     },
@@ -149,7 +149,7 @@ impl UvfitsReader {
             let mut tile_xyzs: Vec<XyzGeodetic> = Vec::with_capacity(total_num_tiles);
             let mut average_xyz = XyzGeocentric::default();
             for i in 0..total_num_tiles {
-                let fits_xyz = read_cell_array::<3>(
+                let fits_xyz = fits_read_cell_f64_array::<3>(
                     &mut uvfits_fptr,
                     &antenna_table_hdu,
                     "STABXYZ",
@@ -1845,57 +1845,3 @@ impl Indices {
         })
     }
 }
-
-/// Pull out fits array-in-a-cell values; useful for e.g. STABXYZ. This function
-/// assumes that the output datatype is f64, and that the fits datatype is
-/// TDOUBLE, so it is not to be used generally!
-fn read_cell_array<const NUM_ELEM: usize>(
-    fits_ptr: &mut fitsio::FitsFile,
-    hdu: &fitsio::hdu::FitsHdu,
-    col_name: &'static str,
-    row: i64,
-) -> Result<[f64; NUM_ELEM], UvfitsReadError> {
-    unsafe {
-        // With the column name, get the column number.
-        let mut status = 0;
-        let mut col_num = -1;
-        let keyword = std::ffi::CString::new(col_name).expect("CString::new failed");
-        // ffgcno = fits_get_colnum
-        fitsio_sys::ffgcno(
-            fits_ptr.as_raw(),
-            0,
-            keyword.as_ptr() as *mut c_char,
-            &mut col_num,
-            &mut status,
-        );
-        // Check the status.
-        fits_check_status(status).map_err(|err| UvfitsReadError::ReadCellArray {
-            col_name,
-            hdu_num: hdu.number + 1,
-            err,
-        })?;
-
-        // Now get the specified row from that column.
-        let mut array = [0.0; NUM_ELEM];
-        // ffgcv = fits_read_col
-        fitsio_sys::ffgcv(
-            fits_ptr.as_raw(),
-            82, // TDOUBLE (fitsio.h)
-            col_num,
-            row + 1,
-            1,
-            NUM_ELEM.try_into().expect("not larger than i64::MAX"),
-            std::ptr::null_mut(),
-            array.as_mut_ptr().cast(),
-            &mut 0,
-            &mut status,
-        );
-        fits_check_status(status).map_err(|err| UvfitsReadError::ReadCellArray {
-            col_name,
-            hdu_num: hdu.number + 1,
-            err,
-        })?;
-
-        Ok(array)
-    }
-}