From 04253ac73a7f2a4f4de73722f541eca99d96e8fd Mon Sep 17 00:00:00 2001 From: Dev Null Date: Sun, 7 Jul 2024 17:19:15 +0800 Subject: [PATCH] fixes #27 fits sourcelist support (#33) * Add support for reading sky models in FITS files. * implement fits srclist writing * tests for parsing jack, gleam and lobes sourcelists * tests for shapelets, lists, multi-components * add fits source list documentation --------- Co-authored-by: Christopher H. Jordan --- CHANGELOG.md | 8 +- build.rs | 2 +- examples/read_fits_srclist.py | 6 + mdbook/src/defs/source_list_fits.md | 219 ++++++++ mdbook/src/defs/source_lists.md | 1 + src/cli/error.rs | 5 +- src/cli/srclist/verify.rs | 59 +- src/srclist/error.rs | 6 + src/srclist/fits/mod.rs | 15 + src/srclist/fits/read.rs | 778 ++++++++++++++++++++++++++ src/srclist/fits/write.rs | 330 +++++++++++ src/srclist/general_tests.rs | 267 ++++++++- src/srclist/mod.rs | 4 + src/srclist/read.rs | 49 +- src/srclist/types/flux_density/mod.rs | 2 +- src/srclist/write.rs | 3 +- test_files/gen_fits_srclists.py | 177 ++++++ test_files/gleam.fits | Bin 0 -> 8640 bytes test_files/gleam.yaml | 50 ++ test_files/jack.fits | 1 + test_files/jack.yaml | 119 ++++ test_files/lobes.fits | 1 + 22 files changed, 2050 insertions(+), 52 deletions(-) create mode 100644 examples/read_fits_srclist.py create mode 100644 mdbook/src/defs/source_list_fits.md create mode 100644 src/srclist/fits/mod.rs create mode 100644 src/srclist/fits/read.rs create mode 100644 src/srclist/fits/write.rs create mode 100755 test_files/gen_fits_srclists.py create mode 100644 test_files/gleam.fits create mode 100644 test_files/gleam.yaml create mode 100644 test_files/jack.fits create mode 100644 test_files/jack.yaml create mode 100644 test_files/lobes.fits diff --git a/CHANGELOG.md b/CHANGELOG.md index 50d477cc..09079d7d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,14 +7,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [0.4.0] - 2024-06-19 +### Added +- fits sourcelist support (including shapelets for Jack-style fits) +- hyperbeam@0.9.3 built@0.7 marlu@0.11.0 mwalib@1.3.3 birli@0.11.0 + ### Fixed - rocm6 support -- a bunch of really nasty segfaults that took a big toll of my sanity +- a bunch of really nasty segfaults that took a big toll on my sanity - Huge thanks to @robotopia for fixing https://github.com/MWATelescope/mwa_hyperbeam/issues/9 via hyperbeam 0.9.0 - performance optimizations in hyperbeam 0.9.3 -### Added -- hyperbeam@0.9.3 built@0.7 marlu@0.11.0 mwalib@1.3.3 birli@0.11.0 ## [0.3.0] - 2023-09-27 ### Added diff --git a/build.rs b/build.rs index 26ff35ad..e00b5313 100644 --- a/build.rs +++ b/build.rs @@ -184,7 +184,7 @@ mod gpu { "cargo:warning=HIP_FLAGS set from env {}", p.to_string_lossy() ); - hip_target.flag(&p.to_string_lossy()); + hip_target.flag(&*p.to_string_lossy()); } println!("cargo:rerun-if-env-changed=ROCM_VER"); diff --git a/examples/read_fits_srclist.py b/examples/read_fits_srclist.py new file mode 100644 index 00000000..77b15d92 --- /dev/null +++ b/examples/read_fits_srclist.py @@ -0,0 +1,6 @@ +from astropy.io import fits +import sys +from tabulate import tabulate + +for hdu in fits.open(sys.argv[-1])[1:]: + print(tabulate(hdu.data, headers=[c.name for c in hdu.columns], tablefmt="github")) \ No newline at end of file diff --git a/mdbook/src/defs/source_list_fits.md b/mdbook/src/defs/source_list_fits.md new file mode 100644 index 00000000..5f22d4fa --- /dev/null +++ b/mdbook/src/defs/source_list_fits.md @@ -0,0 +1,219 @@ +# FITS source list formats + +There are three supported fits file formats: +- LoBES: used in LoBES catalogue +- Jack: extended LoBES format for Jack Line's sourcelist repository, . +- Gleam: used in GLEAM-X pipeline + +These formats differ mostly in the names of columns, and component and flux types +supported. *LoBES* fits files support point, and Gaussian components with list, +power law and curved power law flux density models. *Jack* fits files extend the +LoBES format with an additional table for shapelet coefficients. *Gleam* fits are +similar to LoBES fits, but with different column names, and combine power law and +curved power law flux density models into a just two columns. + +More info from [woden docs](https://woden.readthedocs.io/en/latest/operating_principles/skymodel.html) + +## Source posititons + +Coordinates are right ascension (RA) and declination, both with units of degrees +in the J2000 epoch. All frequencies are in Hz and all flux densities are in Jy. + +Jack and LoBES fits formats use the columns `RA` and `DEC` for source positions, +while Gleam fits files use `RAJ2000` and `DEJ2000`. + +## Component types + +Jack and LoBES fits formats use the column `COMP_TYPE` for component types: +- `P` for point +- `G` for Gaussian +- `S` for shapelet (Jack only) + +Jack and LoBES fits formats use the columns `MAJOR_DC`, `MINOR_DC` and `PA_DC` +for Gaussian component sizes and position angles (in degrees), while Gleam +fits files use `a`, `b` (arcseconds) and `pa` (degrees). + +In an image space where RA increases from right to left (i.e. bigger RA +values are on the left), position angles rotate counter clockwise. A +position angle of 0 has the major axis aligned with the declination axis. + +## Flux density models + +Jack and LoBES fits formats use the column `MOD_TYPE` for flux density types: +- `pl` for power law +- `cpl` for curved power law +- `nan` for lists + +Jack and LoBES fits formats use the columns `NORM_COMP_PL` and `ALPHA_PL` for +power law flux density normalisation and spectral index; and `NORM_COMP_CPL`, +`ALPHA_CPL` and `CURVE_CPL` for curved power law flux density normalisation, +while Gleam fits files use `S_200`, `alpha` and `beta`. + +A reference frequency of 200MHz is assumed in all fits files. + +Jack and LoBES fits formats use the columns `INT_FLXnnn` for integrated flux +densities in Jy at frequencies `nnn` MHz, while Gleam fits files use only `s_200`. +These columns are used to construct flux lists if power law information is +missing, or `MOD_TYPE` is `nan`. + +Only Stokes I can be specified in fits sourcelists, Stokes Q, U and V are +assumed to have values of 0. + +## Examples + +Example Python code to display these files is in the [examples +directory](https://github.com/MWATelescope/mwa_hyperdrive/tree/main/examples). + +e.g. `python examples/read_fits_srclist.py test_files/jack.fits` + +| UNQ_SOURCE_ID | NAME | RA | DEC | INT_FLX100 | INT_FLX150 | INT_FLX200 | MAJOR_DC | MINOR_DC | PA_DC | MOD_TYPE | COMP_TYPE | NORM_COMP_PL | ALPHA_PL | NORM_COMP_CPL | ALPHA_CPL | CURVE_CPL | +|-----------------|---------------|------|-------|--------------|--------------|--------------|------------|------------|---------|------------|-------------|----------------|------------|-----------------|-------------|-------------| +| point-list | point-list_C0 | 0 | 1 | 3 | 2 | 1 | 0 | 0 | 0 | nan | P | 1 | 0 | 0 | 0 | 0 | +| point-pl | point-pl_C0 | 1 | 2 | 3.5 | 2.5 | 2 | 0 | 0 | 0 | pl | P | 2 | -0.8 | 0 | 0 | 0 | +| point-cpl | point-cpl_C0 | 3 | 4 | 5.6 | 3.8 | 3 | 0 | 0 | 0 | cpl | P | 0 | 0 | 3 | -0.9 | 0.2 | +| gauss-list | gauss-list_C0 | 0 | 1 | 3 | 2 | 1 | 20 | 10 | 75 | nan | G | 1 | 0 | 0 | 0 | 0 | +| gauss-pl | gauss-pl_C0 | 1 | 2 | 3.5 | 2.5 | 2 | 20 | 10 | 75 | pl | G | 2 | -0.8 | 0 | 0 | 0 | +| gauss-cpl | gauss-cpl_C0 | 3 | 4 | 5.6 | 3.8 | 3 | 20 | 10 | 75 | cpl | G | 0 | 0 | 3 | -0.9 | 0.2 | +| shape-pl | shape-pl_C0 | 1 | 2 | 3.5 | 2.5 | 2 | 20 | 10 | 75 | pl | S | 2 | -0.8 | 0 | 0 | 0 | +| shape-pl | shape-pl_C1 | 1 | 2 | 3.5 | 2.5 | 2 | 20 | 10 | 75 | pl | S | 2 | -0.8 | 0 | 0 | 0 | + +| NAME | N1 | N2 | COEFF | +|-------------|------|------|---------| +| shape-pl_C0 | 0 | 0 | 0.9 | +| shape-pl_C0 | 0 | 1 | 0.2 | +| shape-pl_C0 | 1 | 0 | -0.2 | +| shape-pl_C1 | 0 | 0 | 0.8 | + +e.g. `python examples/read_fits_srclist.py test_files/gleam.fits` + +| Name | RAJ2000 | DEJ2000 | S_200 | alpha | beta | a | b | pa | +|-----------|-----------|-----------|---------|---------|--------|-------|-------|------| +| point-pl | 1 | 2 | 2 | -0.8 | 0 | 0 | 0 | 0 | +| point-cpl | 3 | 4 | 3 | -0.9 | 0.2 | 0 | 0 | 0 | +| gauss-pl | 1 | 2 | 2 | -0.8 | 0 | 72000 | 36000 | 75 | +| gauss-cpl | 3 | 4 | 3 | -0.9 | 0.2 | 72000 | 36000 | 75 | + +these are both equivalent to the following YAML file (ignoring shapelets and +lists for the gleam example): + +```yaml +point-list: +- ra: 0.0 + dec: 1.0 + comp_type: point + flux_type: + list: + - freq: 100000000.0 + i: 3.0 + - freq: 150000000.0 + i: 2.0 + - freq: 200000000.0 + i: 1.0 +point-pl: +- ra: 1.0 + dec: 2.0 + comp_type: point + flux_type: + power_law: + si: -0.8 + fd: + freq: 200000000.0 + i: 2.0 +point-cpl: +- ra: 3.0000000000000004 + dec: 4.0 + comp_type: point + flux_type: + curved_power_law: + si: -0.9 + fd: + freq: 200000000.0 + i: 3.0 + q: 0.2 +gauss-list: +- ra: 0.0 + dec: 1.0 + comp_type: + gaussian: + maj: 72000.0 + min: 36000.0 + pa: 75.0 + flux_type: + list: + - freq: 100000000.0 + i: 3.0 + - freq: 150000000.0 + i: 2.0 + - freq: 200000000.0 + i: 1.0 +gauss-pl: +- ra: 1.0 + dec: 2.0 + comp_type: + gaussian: + maj: 72000.0 + min: 36000.0 + pa: 75.0 + flux_type: + power_law: + si: -0.8 + fd: + freq: 200000000.0 + i: 2.0 +gauss-cpl: +- ra: 3.0000000000000004 + dec: 4.0 + comp_type: + gaussian: + maj: 72000.0 + min: 36000.0 + pa: 75.0 + flux_type: + curved_power_law: + si: -0.9 + fd: + freq: 200000000.0 + i: 3.0 + q: 0.2 +shape-pl: +- ra: 1.0 + dec: 2.0 + comp_type: + shapelet: + maj: 72000.0 + min: 36000.0 + pa: 75.0 + coeffs: + - n1: 0 + n2: 0 + value: 0.9 + - n1: 0 + n2: 1 + value: 0.2 + - n1: 1 + n2: 0 + value: -0.2 + flux_type: + power_law: + si: -0.8 + fd: + freq: 200000000.0 + i: 2.0 +- ra: 1.0 + dec: 2.0 + comp_type: + shapelet: + maj: 72000.0 + min: 36000.0 + pa: 75.0 + coeffs: + - n1: 0 + n2: 0 + value: 0.8 + flux_type: + power_law: + si: -0.8 + fd: + freq: 200000000.0 + i: 2.0 +``` \ No newline at end of file diff --git a/mdbook/src/defs/source_lists.md b/mdbook/src/defs/source_lists.md index d6be26a3..1d471457 100644 --- a/mdbook/src/defs/source_lists.md +++ b/mdbook/src/defs/source_lists.md @@ -23,6 +23,7 @@ types. - [`hyperdrive` format](source_list_hyperdrive.md) - [André Offringa (`ao`) format](source_list_ao.md) - [`RTS` format](source_list_rts.md) +- [Jack, Gleam or LoBES style fits](source_list_fits.md) ~~~ ~~~admonish info title="Conversion" diff --git a/src/cli/error.rs b/src/cli/error.rs index fbdf3748..2faaff61 100644 --- a/src/cli/error.rs +++ b/src/cli/error.rs @@ -368,7 +368,10 @@ impl From for HyperdriveError { | WriteSourceListError::InvalidHyperdriveFormat(_) | WriteSourceListError::Sexagesimal(_) => Self::Srclist(s), WriteSourceListError::IO(e) => Self::from(e), - WriteSourceListError::Yaml(_) | WriteSourceListError::Json(_) => Self::Generic(s), + WriteSourceListError::Yaml(_) + | WriteSourceListError::Json(_) + | WriteSourceListError::Fitsio(_) + | WriteSourceListError::Fits(_) => Self::Generic(s), } } } diff --git a/src/cli/srclist/verify.rs b/src/cli/srclist/verify.rs index 83dc7d32..b129a5a6 100644 --- a/src/cli/srclist/verify.rs +++ b/src/cli/srclist/verify.rs @@ -16,8 +16,8 @@ use log::info; use crate::{ cli::common::{display_warnings, SOURCE_LIST_INPUT_TYPE_HELP}, srclist::{ - ao, hyperdrive, read::read_source_list_file, rts, woden, ComponentCounts, SourceListType, - SrclistError, + ao, fits, hyperdrive, read::read_source_list_file, rts, woden, ComponentCounts, + SourceListType, SrclistError, }, HyperdriveError, }; @@ -67,24 +67,45 @@ fn verify>( info!("{}:", source_list.as_ref().display()); let (sl, sl_type) = if let Some(input_type) = input_type { - let mut buf = std::io::BufReader::new(File::open(source_list)?); let result = match input_type { - SourceListType::Hyperdrive => crate::misc::expensive_op( - || hyperdrive::source_list_from_yaml(&mut buf), - "Still reading source list file", - ), - SourceListType::AO => crate::misc::expensive_op( - || ao::parse_source_list(&mut buf), - "Still reading source list file", - ), - SourceListType::Rts => crate::misc::expensive_op( - || rts::parse_source_list(&mut buf), - "Still reading source list file", - ), - SourceListType::Woden => crate::misc::expensive_op( - || woden::parse_source_list(&mut buf), - "Still reading source list file", - ), + SourceListType::Hyperdrive => { + let mut buf = std::io::BufReader::new(File::open(source_list)?); + crate::misc::expensive_op( + || hyperdrive::source_list_from_yaml(&mut buf), + "Still reading source list file", + ) + } + SourceListType::Fits => { + let source_list = source_list.as_ref(); + let sl = crate::misc::expensive_op( + || fits::parse_source_list(source_list), + "Still reading source list file", + ) + .unwrap(); + // TODO: Proper error handling + Ok(sl) + } + SourceListType::AO => { + let mut buf = std::io::BufReader::new(File::open(source_list)?); + crate::misc::expensive_op( + || ao::parse_source_list(&mut buf), + "Still reading source list file", + ) + } + SourceListType::Rts => { + let mut buf = std::io::BufReader::new(File::open(source_list)?); + crate::misc::expensive_op( + || rts::parse_source_list(&mut buf), + "Still reading source list file", + ) + } + SourceListType::Woden => { + let mut buf = std::io::BufReader::new(File::open(source_list)?); + crate::misc::expensive_op( + || woden::parse_source_list(&mut buf), + "Still reading source list file", + ) + } }; match result { Ok(sl) => (sl, input_type), diff --git a/src/srclist/error.rs b/src/srclist/error.rs index 5f54b24f..aac9b282 100644 --- a/src/srclist/error.rs +++ b/src/srclist/error.rs @@ -309,6 +309,12 @@ pub(crate) enum WriteSourceListError { #[error(transparent)] Sexagesimal(#[from] marlu::sexagesimal::SexagesimalError), + #[error(transparent)] + Fitsio(#[from] fitsio::errors::Error), + + #[error(transparent)] + Fits(#[from] crate::io::read::fits::FitsError), + /// An IO error. #[error(transparent)] IO(#[from] std::io::Error), diff --git a/src/srclist/fits/mod.rs b/src/srclist/fits/mod.rs new file mode 100644 index 00000000..113209a7 --- /dev/null +++ b/src/srclist/fits/mod.rs @@ -0,0 +1,15 @@ +// This Source Code Form is subject to the terms of the Mozilla Public +// 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/. + +//! Code to handle FITS source list files. + +// The reference frequency of the power laws. +const REF_FREQ_HZ: f64 = 200e6; + +mod read; +mod write; + +// Re-exports. +pub(crate) use read::parse_source_list; +pub(crate) use write::write_source_list_jack; diff --git a/src/srclist/fits/read.rs b/src/srclist/fits/read.rs new file mode 100644 index 00000000..57e69a4f --- /dev/null +++ b/src/srclist/fits/read.rs @@ -0,0 +1,778 @@ +// This Source Code Form is subject to the terms of the Mozilla Public +// 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 std::{collections::HashMap, path::Path}; + +use fitsio::{hdu::FitsHdu, FitsFile}; +use indexmap::IndexMap; +use itertools::izip; +use log::debug; +use marlu::RADec; +use vec1::Vec1; + +// The reference frequency of the power laws and curved power laws. +use super::REF_FREQ_HZ; +use crate::{ + io::read::fits::{fits_open, fits_open_hdu, FitsError}, + srclist::{ + ComponentType, FluxDensity, FluxDensityType, ShapeletCoeff, Source, SourceComponent, + SourceList, + }, +}; + +/// A supported "flux type", but without any data. +#[derive(Clone, Copy, Debug)] +enum FluxType { + PowerLaw, + CurvedPowerLaw, + List, +} + +/// Wrap a FITS call with our error. +macro_rules! fe { + ($file:expr, $result:expr) => {{ + $result.map_err(|e| FitsError::Fitsio { + fits_error: Box::new(e), + fits_filename: $file.to_path_buf().into_boxed_path(), + hdu_description: "1".to_string().into_boxed_str(), + source_file: file!(), + source_line: line!(), + source_column: column!(), + })? + }}; +} + +pub(crate) fn parse_source_list(file: &Path) -> Result { + let mut fptr = fits_open(file)?; + // We assume everything is on HDU 2. + let hdu = fits_open_hdu(&mut fptr, 1)?; + + // Get the fits column names. + let col_names: Vec = match &hdu.info { + fitsio::hdu::HduInfo::TableInfo { + column_descriptions, + num_rows: _, + } => column_descriptions + .iter() + .map(|cd| cd.name.clone()) + .collect(), + fitsio::hdu::HduInfo::ImageInfo { .. } => todo!(), + fitsio::hdu::HduInfo::AnyInfo => todo!(), + }; + + // Try to determine the format of the FITS file. + if col_names.iter().any(|name| *name == "UNQ_SOURCE_ID") { + debug!("Attempting to read the first NAME"); + let first_src_id = fe!(file, hdu.read_cell_value::(&mut fptr, "NAME", 0)); + if first_src_id.contains("_GID") { + debug!("I reckon this is a 'LoBES' source list"); + parse_lobes_source_list(file, fptr, hdu, col_names) + } else { + debug!("I reckon this is a 'Jack FITS' source list"); + parse_jack_source_list(file, fptr, hdu, col_names) + } + } else { + debug!("I reckon this is a 'GLEAM FITS' source list"); + parse_gleam_x_source_list(file, fptr, hdu, col_names) + } +} + +/// FITS file source lists are pretty similar. This struct +/// contains the common columns. +struct CommonCols { + unq_source_id: Vec, + names: Vec, + ra_degrees: Vec, + dec_degrees: Vec, + majors: Vec, + minors: Vec, + pas: Vec, + shapelet_sources: Option>>, + comp_types: Vec, + flux_types: Vec, + power_law_stokes_is: Vec, + power_law_alphas: Vec, + curved_power_law_stokes_is: Vec, + curved_power_law_alphas: Vec, + curved_power_law_qs: Vec, + list_flux_densities: Vec>, + list_flux_density_freqs: Vec, +} + +impl CommonCols { + fn new( + file: &Path, + fptr: &mut FitsFile, + hdu: &FitsHdu, + col_names: &[String], + ) -> Result { + macro_rules! read_optional_col { + ($possible_col_names: expr) => {{ + let mut maybe_col = None; + for possible_col_name in $possible_col_names { + if col_names + .iter() + .any(|col_name| col_name == possible_col_name) + { + maybe_col = Some(fe!(file, hdu.read_col(fptr, possible_col_name))); + } + } + if !maybe_col.is_some() { + debug!("None of {:?} were available columns!", $possible_col_names) + } + maybe_col + }}; + } + macro_rules! read_mandatory_col { + ($possible_col_names: expr) => {{ + read_optional_col!($possible_col_names).unwrap_or_else(|| { + panic!("None of {:?} were available columns!", $possible_col_names) + }) + }}; + } + + let unq_source_id = if col_names.iter().any(|col_name| col_name == "UNQ_SOURCE_ID") { + fe!(file, hdu.read_col(fptr, "UNQ_SOURCE_ID")) + } else { + vec![] + }; + + let names = read_mandatory_col!(["NAME", "Name"]); + let ra_degrees = read_mandatory_col!(["RA", "RAJ2000"]); + let dec_degrees = read_mandatory_col!(["DEC", "DEJ2000"]); + let majors = read_mandatory_col!(["MAJOR_DC", "a"]); + let minors = read_mandatory_col!(["MINOR_DC", "b"]); + let pas = read_mandatory_col!(["PA_DC", "pa"]); + + // Get any shapelet info ready. We assume that the info lives in HDU 3 + // (index 2 in sane languages), and if there's an error, we assume it's + // because the HDU isn't present. + let hdu_shapelets = fits_open_hdu(fptr, 2).ok(); + let shapelet_sources = match hdu_shapelets { + Some(hdu) => { + let mut map = HashMap::new(); + let shapelet_sources: Vec = fe!(file, hdu.read_col(fptr, "NAME")); + let shapelet_n1s: Vec = fe!(file, hdu.read_col(fptr, "N1")); + let shapelet_n2s: Vec = fe!(file, hdu.read_col(fptr, "N2")); + let shapelet_coeff_values: Vec = fe!(file, hdu.read_col(fptr, "COEFF")); + + // Pre-allocate vectors for the shapelet coeffs. + let mut count = 0; + let mut last_source = None; + for this_source in &shapelet_sources { + match last_source.as_mut() { + None => { + last_source = Some(this_source); + count += 1; + continue; + } + + Some(last_source) => { + if *last_source == this_source { + count += 1; + continue; + } else { + map.insert(last_source.clone(), Vec::with_capacity(count)); + count = 0; + *last_source = this_source; + } + } + } + } + // Don't forget the last source. + if let Some(last_source) = last_source { + map.insert(last_source.clone(), Vec::with_capacity(count)); + } + + // Now populate the shapelet coeffs associated with the sources. + for (source, n1, n2, value) in izip!( + shapelet_sources, + shapelet_n1s, + shapelet_n2s, + shapelet_coeff_values + ) { + let coeffs = map.get_mut(&source).expect("was populated previously"); + coeffs.push(ShapeletCoeff { + n1: n1.try_into().expect("n1 is not larger than u8::MAX"), + n2: n2.try_into().expect("n2 is not larger than u8::MAX"), + value, + }) + } + + Some(map) + } + + None => None, + }; + + let comp_types = if col_names.iter().any(|col_name| col_name == "COMP_TYPE") { + // This is pretty gross, but, so is handling bespoke formats. + let comp_types_as_strings: Vec = fe!(file, hdu.read_col(fptr, "COMP_TYPE")); + comp_types_as_strings + .into_iter() + .map(|mut s| s.pop().expect("COMP_TYPE strings aren't empty")) + .collect() + } else { + // We need to determine the component types here. + let mut comp_types = Vec::with_capacity(names.len()); + for (i_row, name) in names.iter().enumerate() { + if let Some(shapelet_sources) = shapelet_sources.as_ref() { + if shapelet_sources.contains_key(name) { + comp_types.push('S'); + continue; + } + } else { + let gauss_maj: f64 = majors[i_row]; + let gauss_min: f64 = minors[i_row]; + let gauss_pa: f64 = pas[i_row]; + if (gauss_maj.is_nan() && gauss_min.is_nan() && gauss_pa.is_nan()) + || (gauss_maj.abs() < f64::EPSILON + && gauss_min.abs() < f64::EPSILON + && gauss_pa.abs() < f64::EPSILON) + { + comp_types.push('P'); + } else { + comp_types.push('G'); + } + } + } + comp_types + }; + + let power_law_stokes_is: Vec = read_mandatory_col!(["NORM_COMP_PL", "S_200"]); + let power_law_alphas = read_mandatory_col!(["ALPHA_PL", "alpha"]); + + let (curved_power_law_stokes_is, curved_power_law_alphas, curved_power_law_qs): ( + Vec, + Vec, + Vec, + ) = if col_names.iter().any(|col_name| col_name == "NORM_COMP_CPL") { + ( + fe!(file, hdu.read_col(fptr, "NORM_COMP_CPL")), + fe!(file, hdu.read_col(fptr, "ALPHA_CPL")), + fe!(file, hdu.read_col(fptr, "CURVE_CPL")), + ) + } else if col_names.iter().any(|col_name| col_name == "beta") { + (vec![], vec![], fe!(file, hdu.read_col(fptr, "beta"))) + } else { + (vec![], vec![], vec![]) + }; + + // For "list" flux density types, first, find all the relevant + // columns. Then, pull out all their values. + let (list_flux_densities, list_flux_density_freqs): (Vec>, Vec) = if col_names + .iter() + .any(|col_name| col_name.starts_with("INT_FLX")) + { + { + let mut flux_densities = Vec::with_capacity(32); + let mut flux_density_freqs = Vec::with_capacity(32); + #[allow(non_snake_case)] + for (col_name, int_flx_freq_MHz) in col_names.iter().filter_map(|col_name| { + col_name + .strip_prefix("INT_FLX") + .map(|suffix| (col_name.as_str(), suffix)) + }) { + flux_densities.push(fe!(file, hdu.read_col(fptr, col_name))); + flux_density_freqs + .push(int_flx_freq_MHz.parse::().expect("is a number") * 1e6); + } + (flux_densities, flux_density_freqs) + } + } else { + (vec![], vec![]) + }; + + let flux_types = if col_names.iter().any(|col_name| col_name == "MOD_TYPE") { + let mod_types: Vec = fe!(file, hdu.read_col(fptr, "MOD_TYPE")); + mod_types.into_iter().enumerate().map(|(i_row, mod_type)| match mod_type.as_str() { + "nan" => FluxType::List, + "pl" => FluxType::PowerLaw, + "cpl" => FluxType::CurvedPowerLaw, + t => panic!("Got '{t}' in row {i_row} of the 'MOD_TYPE' column, which is none of 'nan', 'pl' or 'cpl'"), + }).collect() + } else { + // MOD_TYPE isn't here. We assume that there are always power laws, + // but check for curved power laws and lists. + if curved_power_law_stokes_is.is_empty() && list_flux_densities.is_empty() { + vec![FluxType::PowerLaw; power_law_stokes_is.len()] + } else { + // We have to iterate over everything... + let mut results = Vec::with_capacity(power_law_stokes_is.len()); + let mut indices_to_check = Vec::with_capacity(power_law_stokes_is.len()); + // We assume that if a CPL q is present, this should be a CPL. + if !curved_power_law_qs.is_empty() { + for (i, q) in curved_power_law_qs.iter().copied().enumerate() { + if !q.is_nan() && q != 0.0 { + results.push(FluxType::CurvedPowerLaw); + } else { + indices_to_check.push(i); + } + } + } else { + for i in 0..power_law_stokes_is.len() { + indices_to_check.push(i); + } + } + // Do the same thing for PL Stokes Is. + for stokes_i in indices_to_check.into_iter().map(|i| power_law_stokes_is[i]) { + if !stokes_i.is_nan() && stokes_i != 0.0 { + results.push(FluxType::PowerLaw); + } else { + results.push(FluxType::List); + } + } + + results + } + }; + + Ok(Self { + unq_source_id, + names, + ra_degrees, + dec_degrees, + majors, + minors, + pas, + shapelet_sources, + comp_types, + flux_types, + power_law_stokes_is, + power_law_alphas, + curved_power_law_stokes_is, + curved_power_law_alphas, + curved_power_law_qs, + list_flux_densities, + list_flux_density_freqs, + }) + } +} + +fn parse_lobes_source_list( + file: &Path, + mut fptr: FitsFile, + hdu: FitsHdu, + col_names: Vec, +) -> Result { + let CommonCols { + unq_source_id: _, + names, + ra_degrees, + dec_degrees, + majors, + minors, + pas, + flux_types, + power_law_stokes_is, + power_law_alphas, + curved_power_law_stokes_is, + curved_power_law_alphas, + curved_power_law_qs, + list_flux_densities, + list_flux_density_freqs, + .. + } = CommonCols::new(file, &mut fptr, &hdu, &col_names)?; + + // UNQ_SOURCE_ID comes in as strings of ints. Parse them. + + let mut map = IndexMap::with_capacity(names.len()); + // Get all of the source names. + for (i_row, name) in names.into_iter().enumerate() { + let prefix_len = name.rsplit_once("_GID").expect("contains '_GID'").0.len(); + // Here, we assume pure ASCII. + let name = String::from_utf8({ + let mut bytes = name.into_bytes(); + bytes.resize(prefix_len, 0); + bytes + }) + .expect("is valid"); + + map.entry(name) + .and_modify(|v: &mut Vec| v.push(i_row)) + .or_insert(vec![i_row]); + } + + // Find all of the components that belong with this source and populate it. + let mut source_list = IndexMap::with_capacity(map.len()); + for (name, rows) in map { + let src = source_list.entry(name).or_insert(vec![]); + + for i_row in rows { + let radec = RADec::from_degrees(ra_degrees[i_row], dec_degrees[i_row]); + let maj = majors[i_row]; + let min = minors[i_row]; + let pa = pas[i_row]; + let comp_type = if (maj.is_nan() && min.is_nan() && pa.is_nan()) + || (maj.abs() < f64::EPSILON && min.abs() < f64::EPSILON && pa.abs() < f64::EPSILON) + { + ComponentType::Point + } else { + ComponentType::Gaussian { + maj: maj.to_radians(), + min: min.to_radians(), + pa: pa.to_radians(), + } + }; + + let flux_type = match flux_types[i_row] { + FluxType::List => { + let mut list_fds = Vec::with_capacity(list_flux_densities.len()); + for (fds, &freq) in list_flux_densities + .iter() + .zip(list_flux_density_freqs.iter()) + { + let fd = fds[i_row]; + if !fd.is_nan() { + list_fds.push(FluxDensity { + freq, + i: fd, + q: 0.0, + u: 0.0, + v: 0.0, + }); + } + } + FluxDensityType::List(Vec1::try_from_vec(list_fds).unwrap_or_else(|_| { + panic!("No valid flux densities were provided on row {i_row}") + })) + } + + FluxType::PowerLaw => { + let si = power_law_alphas[i_row]; + let i = power_law_stokes_is[i_row]; + FluxDensityType::PowerLaw { + si, + fd: FluxDensity { + freq: REF_FREQ_HZ, + i, + q: 0.0, + u: 0.0, + v: 0.0, + }, + } + } + + FluxType::CurvedPowerLaw => { + let si = curved_power_law_alphas[i_row]; + let i = curved_power_law_stokes_is[i_row]; + let q = curved_power_law_qs[i_row]; + FluxDensityType::CurvedPowerLaw { + si, + fd: FluxDensity { + freq: REF_FREQ_HZ, + i, + q: 0.0, + u: 0.0, + v: 0.0, + }, + q, + } + } + }; + + src.push(SourceComponent { + radec, + comp_type, + flux_type, + }) + } + } + + // Box all of the vectors. + let source_list = source_list + .into_iter() + .map(|(name, comps)| { + ( + name, + Source { + components: comps.into_boxed_slice(), + }, + ) + }) + .collect::(); + + Ok(source_list) +} + +fn parse_jack_source_list( + file: &Path, + mut fptr: FitsFile, + hdu: FitsHdu, + col_names: Vec, +) -> Result { + let CommonCols { + unq_source_id: src_names, + names: comp_names, + ra_degrees, + dec_degrees, + majors, + minors, + pas, + mut shapelet_sources, + comp_types, + flux_types, + power_law_stokes_is, + power_law_alphas, + curved_power_law_stokes_is, + curved_power_law_alphas, + curved_power_law_qs, + list_flux_densities, + list_flux_density_freqs, + } = CommonCols::new(file, &mut fptr, &hdu, &col_names)?; + let mut map = IndexMap::with_capacity(src_names.len()); + // Get all of the source names. + for name in src_names { + map.entry(name).or_insert(vec![]); + } + + // Find all of the components that belong with this source and populate + // it. + for ( + i_row, + (ra_degrees, dec_degrees, flux_type, comp_name, comp_type, gauss_maj, gauss_min, gauss_pa), + ) in izip!( + ra_degrees.into_iter(), + dec_degrees.into_iter(), + flux_types.into_iter(), + comp_names.into_iter(), + comp_types.into_iter(), + majors.into_iter(), + minors.into_iter(), + pas.into_iter(), + ) + .enumerate() + { + let prefix = comp_name + .rsplit_once("_C") + .unwrap_or_else(|| panic!("{comp_name:?} does not contain '_C'")) + .0; + let src_comps = map.get_mut(prefix).unwrap_or_else(|| { + panic!("Component '{comp_name}' couldn't be matched against any of the UNQ_SOURCE_ID") + }); + + let radec = RADec::from_degrees(ra_degrees, dec_degrees); + let gaussian_params_arent_available = + (gauss_maj.is_nan() && gauss_min.is_nan() && gauss_pa.is_nan()) + || (gauss_maj.abs() < f64::EPSILON + && gauss_min.abs() < f64::EPSILON + && gauss_pa.abs() < f64::EPSILON); + let comp_type = match comp_type { + // The COMP_TYPE column indicates a point source + 'P' => { + // Check if Gaussian parameters were provided; if so, complain. + if gaussian_params_arent_available { + ComponentType::Point + } else { + panic!("Gaussian parameters were provided for COMP_TYPE 'P' on row {i_row}"); + } + } + + // The COMP_TYPE column indicates a Gaussian source + 'G' => { + if gaussian_params_arent_available { + panic!("Gaussian parameters weren't provided for COMP_TYPE 'G' on row {i_row}"); + } else { + ComponentType::Gaussian { + maj: gauss_maj.to_radians(), + min: gauss_min.to_radians(), + pa: gauss_pa.to_radians(), + } + } + } + + // The COMP_TYPE column indicates a shapelet source + 'S' => { + if gaussian_params_arent_available { + panic!("Gaussian parameters weren't provided for COMP_TYPE 'S' on row {i_row}"); + } else { + // We need to extract the shapelet coeffs associated with + // this source. + let coeffs = shapelet_sources + .as_mut() + .expect("shapelet coeffs available if a COMP_TYPE 'S' is present") + .remove(&comp_name) + .expect("shapelet coeffs available for COMP_TYPE 'S' source"); + + ComponentType::Shapelet { + maj: gauss_maj.to_radians(), + min: gauss_min.to_radians(), + pa: gauss_pa.to_radians(), + coeffs: coeffs.into_boxed_slice(), + } + } + } + + t => panic!("Got an unexpected COMP_TYPE '{t}' on row {i_row}"), + }; + + let flux_type = match flux_type { + FluxType::List => { + let mut list_fds = Vec::with_capacity(list_flux_densities.len()); + for (fds, &freq) in list_flux_densities + .iter() + .zip(list_flux_density_freqs.iter()) + { + let fd = fds[i_row]; + if !fd.is_nan() { + list_fds.push(FluxDensity { + freq, + i: fd, + q: 0.0, + u: 0.0, + v: 0.0, + }); + } + } + FluxDensityType::List(Vec1::try_from_vec(list_fds).unwrap_or_else(|_| { + panic!("No valid flux densities were provided on row {i_row}") + })) + } + + FluxType::PowerLaw => { + let si = power_law_alphas[i_row]; + let i = power_law_stokes_is[i_row]; + FluxDensityType::PowerLaw { + si, + fd: FluxDensity { + freq: REF_FREQ_HZ, + i, + q: 0.0, + u: 0.0, + v: 0.0, + }, + } + } + + FluxType::CurvedPowerLaw => { + let si = curved_power_law_alphas[i_row]; + let i = curved_power_law_stokes_is[i_row]; + let q = curved_power_law_qs[i_row]; + FluxDensityType::CurvedPowerLaw { + si, + fd: FluxDensity { + freq: REF_FREQ_HZ, + i, + q: 0.0, + u: 0.0, + v: 0.0, + }, + q, + } + } + }; + + src_comps.push(SourceComponent { + radec, + comp_type, + flux_type, + }); + } + + // Box all of the vectors. + let source_list = map + .into_iter() + .map(|(name, comps)| { + if comps.is_empty() { + panic!("No components are against source {name}; this is a programmer error"); + } + + ( + name, + Source { + components: comps.into_boxed_slice(), + }, + ) + }) + .collect::(); + + Ok(source_list) +} + +fn parse_gleam_x_source_list( + file: &Path, + mut fptr: FitsFile, + hdu: FitsHdu, + col_names: Vec, +) -> Result { + let CommonCols { + unq_source_id: _, + names: src_names, + ra_degrees, + dec_degrees, + majors, + minors, + pas, + power_law_stokes_is, + power_law_alphas, + curved_power_law_qs, + .. + } = CommonCols::new(file, &mut fptr, &hdu, &col_names)?; + + let mut source_list = IndexMap::with_capacity(src_names.len()); + + // It appears that each source has one component and everything is a power law. + for (src_name, ra_degrees, dec_degrees, maj, min, pa, stokes_i, si, q) in izip!( + src_names, + ra_degrees, + dec_degrees, + majors, + minors, + pas, + power_law_stokes_is, + power_law_alphas, + curved_power_law_qs + ) { + let radec = RADec::from_degrees(ra_degrees, dec_degrees); + let comp_type = if (maj.is_nan() && min.is_nan() && pa.is_nan()) + || (maj.abs() < f64::EPSILON && min.abs() < f64::EPSILON && pa.abs() < f64::EPSILON) + { + ComponentType::Point + } else { + ComponentType::Gaussian { + maj: (maj / 3600.0).to_radians(), + min: (min / 3600.0).to_radians(), + pa: pa.to_radians(), + } + }; + let flux_type = if q.is_nan() || q.abs() < f64::EPSILON { + FluxDensityType::PowerLaw { + si, + fd: FluxDensity { + freq: REF_FREQ_HZ, + i: stokes_i, + q: 0.0, + u: 0.0, + v: 0.0, + }, + } + } else { + FluxDensityType::CurvedPowerLaw { + si, + fd: FluxDensity { + freq: REF_FREQ_HZ, + i: stokes_i, + q: 0.0, + u: 0.0, + v: 0.0, + }, + q, + } + }; + + source_list.insert( + src_name, + Source { + components: vec![SourceComponent { + radec, + comp_type, + flux_type, + }] + .into_boxed_slice(), + }, + ); + } + + Ok(SourceList::from(source_list)) +} diff --git a/src/srclist/fits/write.rs b/src/srclist/fits/write.rs new file mode 100644 index 00000000..0c079fbc --- /dev/null +++ b/src/srclist/fits/write.rs @@ -0,0 +1,330 @@ +use std::{ffi::CString, path::Path}; + +use super::REF_FREQ_HZ; +use crate::{ + io::read::fits::fits_open_hdu, + srclist::{error::WriteSourceListError, ComponentType, FluxDensityType, SourceList}, +}; +use fitsio::{ + errors::check_status as fits_check_status, + tables::{ColumnDataType, ColumnDescription}, + FitsFile, +}; + +pub(crate) fn write_source_list_jack( + file: &Path, + sl: &SourceList, + num_sources: Option, +) -> Result<(), WriteSourceListError> { + if file.exists() { + std::fs::remove_file(file)?; + } + let mut fptr = FitsFile::create(file).open()?; + let hdu = fits_open_hdu(&mut fptr, 0)?; + let mut status = 0; + + // Signal that we're using long strings. + unsafe { + // ffplsw = fits_write_key_longwarn + fitsio_sys::ffplsw( + fptr.as_raw(), /* I - FITS file pointer */ + &mut status, /* IO - error status */ + ); + fits_check_status(status)?; + } + + // Write the documentation URL as a comment. + unsafe { + let comm = CString::new("The contents of this file are documented at:").unwrap(); + // ffpcom = fits_write_comment + fitsio_sys::ffpcom(fptr.as_raw(), comm.as_ptr(), &mut status); + fits_check_status(status)?; + let comm = CString::new( + "https://mwatelescope.github.io/mwa_hyperdrive/defs/source_list_fits_jack.html", + ) + .unwrap(); + fitsio_sys::ffpcom(fptr.as_raw(), comm.as_ptr(), &mut status); + fits_check_status(status)?; + } + + hdu.write_key( + &mut fptr, + "SOFTWARE", + format!( + "Created by {} v{}", + env!("CARGO_PKG_NAME"), + env!("CARGO_PKG_VERSION") + ), + )?; + + // Write out the current command-line call ("CMDLINE"). + unsafe { + // It's possible that the command-line call has invalid UTF-8. So use + // args_os and attempt to convert to UTF-8 strings. If there are + // problems on the way, don't bother trying to write the CMDLINE key. + match std::env::args_os() + .map(|a| a.into_string()) + .collect::, _>>() + .and_then(|v| CString::new(v.join(" ")).map_err(|_| std::ffi::OsString::from(""))) + { + // This represents failure to convert an argument to UTF-8. + Err(_) => (), + Ok(value) => { + let key_name = CString::new("CMDLINE").unwrap(); + let comment = CString::new("Command-line call").unwrap(); + let mut status = 0; + // ffpkls = fits_write_key_longstr + fitsio_sys::ffpkls( + fptr.as_raw(), /* I - FITS file pointer */ + key_name.as_ptr(), /* I - name of keyword to write */ + value.as_ptr(), /* I - keyword value */ + comment.as_ptr(), /* I - keyword comment */ + &mut status, /* IO - error status */ + ); + fits_check_status(status)?; + } + } + } + + // build the components and shapelets tables + let mut unq_source_id = Vec::new(); + let mut names = Vec::new(); + let mut ra_degrees = Vec::new(); + let mut dec_degrees = Vec::new(); + let mut majors = Vec::new(); + let mut minors = Vec::new(); + let mut pas = Vec::new(); + let mut comp_types = Vec::new(); + let mut mod_types = Vec::new(); + let mut norm_comp_pls = Vec::new(); + let mut alpha_pls = Vec::new(); + let mut norm_comp_cpls = Vec::new(); + let mut alpha_cpls = Vec::new(); + let mut curve_cpls = Vec::new(); + + let mut shapelet_sources = Vec::new(); + let mut shapelet_n1s: Vec = Vec::new(); + let mut shapelet_n2s: Vec = Vec::new(); + let mut shapelet_coeff_values = Vec::new(); + + // get unique flux list frequencies (to nearest MHz) + let mut flux_freqs = vec![]; + for comp in sl.iter().flat_map(|(_, src)| src.components.iter()) { + match &comp.flux_type { + FluxDensityType::List(l) => { + for f in l.iter() { + let freq_mhz = (f.freq / 1e6).round() as i32; + if !flux_freqs.contains(&freq_mhz) { + flux_freqs.push(freq_mhz); + } + } + } + _ => continue, + } + } + flux_freqs.sort(); + let mut flux_lists: Vec> = Vec::new(); + for _ in flux_freqs.iter() { + flux_lists.push(Vec::new()); + } + + let mut max_src_name = 0; + + for (sidx, (src_name, src)) in sl.iter().enumerate() { + if src_name.len() > max_src_name { + max_src_name = src_name.len(); + } + if let Some(num_sources) = num_sources { + if sidx >= num_sources { + break; + } + } + for (cidx, comp) in src.components.iter().enumerate() { + unq_source_id.push(src_name.to_string()); + let comp_name = format!("{src_name}_C{cidx:02}"); + names.push(comp_name.clone()); + ra_degrees.push(comp.radec.ra.to_degrees()); + dec_degrees.push(comp.radec.dec.to_degrees()); + match &comp.comp_type { + ComponentType::Point => { + majors.push(f64::NAN); + minors.push(f64::NAN); + pas.push(f64::NAN); + comp_types.push("P".to_string()); + } + + ComponentType::Gaussian { maj, min, pa } => { + majors.push(maj.to_degrees()); + minors.push(min.to_degrees()); + pas.push(pa.to_degrees()); + comp_types.push("G".to_string()); + } + + ComponentType::Shapelet { + maj, + min, + pa, + coeffs, + } => { + majors.push(maj.to_degrees()); + minors.push(min.to_degrees()); + pas.push(pa.to_degrees()); + comp_types.push("S".to_string()); + + for coeff in coeffs.iter() { + shapelet_sources.push(comp_name.clone()); + shapelet_n1s.push(coeff.n1 as i32); + shapelet_n2s.push(coeff.n2 as i32); + shapelet_coeff_values.push(coeff.value); + } + } + } + match &comp.flux_type { + FluxDensityType::List(l) => { + mod_types.push("nan".to_string()); + for (idx, freq) in flux_freqs.iter().enumerate() { + let flux = l.iter().find(|f| (f.freq / 1e6).round() as i32 == *freq); + match flux { + Some(f) => flux_lists[idx].push(f.i), + None => flux_lists[idx].push(f64::NAN), + } + } + norm_comp_pls.push(f64::NAN); + alpha_pls.push(f64::NAN); + norm_comp_cpls.push(f64::NAN); + alpha_cpls.push(f64::NAN); + curve_cpls.push(f64::NAN); + } + FluxDensityType::CurvedPowerLaw { si, fd, q } => { + mod_types.push("cpl".to_string()); + for (idx, _) in flux_freqs.iter().enumerate() { + flux_lists[idx].push(f64::NAN); + } + if fd.freq != REF_FREQ_HZ { + panic!("Curve Power-law flux densities must be at the reference frequency"); + } + norm_comp_pls.push(f64::NAN); + alpha_pls.push(f64::NAN); + norm_comp_cpls.push(fd.i); + alpha_cpls.push(*si); + curve_cpls.push(*q); + } + FluxDensityType::PowerLaw { si, fd } => { + mod_types.push("pl".to_string()); + for (idx, _) in flux_freqs.iter().enumerate() { + flux_lists[idx].push(f64::NAN); + } + if fd.freq != REF_FREQ_HZ { + panic!("Power-law flux densities must be at the reference frequency"); + } + norm_comp_pls.push(fd.i); + alpha_pls.push(*si); + norm_comp_cpls.push(f64::NAN); + alpha_cpls.push(f64::NAN); + curve_cpls.push(f64::NAN); + } + } + } + } + + // write the components HDU + let mut table_description = vec![ + ColumnDescription::new("UNQ_SOURCE_ID") + .with_type(ColumnDataType::String) + .that_repeats(max_src_name) + .create()?, + ColumnDescription::new("NAME") + .with_type(ColumnDataType::String) + .that_repeats(max_src_name + 4) + .create()?, + ColumnDescription::new("RA") + .with_type(ColumnDataType::Double) + .create()?, + ColumnDescription::new("DEC") + .with_type(ColumnDataType::Double) + .create()?, + ]; + for freq in flux_freqs.iter() { + table_description.push( + ColumnDescription::new(&format!("INT_FLX{freq}")) + .with_type(ColumnDataType::Double) + .create()?, + ); + } + let mut extra = vec![ + ColumnDescription::new("MAJOR_DC") + .with_type(ColumnDataType::Double) + .create()?, + ColumnDescription::new("MINOR_DC") + .with_type(ColumnDataType::Double) + .create()?, + ColumnDescription::new("PA_DC") + .with_type(ColumnDataType::Double) + .create()?, + ColumnDescription::new("MOD_TYPE") + .with_type(ColumnDataType::String) + .that_repeats(3) + .create()?, + ColumnDescription::new("COMP_TYPE") + .with_type(ColumnDataType::String) + .create()?, + ColumnDescription::new("NORM_COMP_PL") + .with_type(ColumnDataType::Double) + .create()?, + ColumnDescription::new("ALPHA_PL") + .with_type(ColumnDataType::Double) + .create()?, + ColumnDescription::new("NORM_COMP_CPL") + .with_type(ColumnDataType::Double) + .create()?, + ColumnDescription::new("ALPHA_CPL") + .with_type(ColumnDataType::Double) + .create()?, + ColumnDescription::new("CURVE_CPL") + .with_type(ColumnDataType::Double) + .create()?, + ]; + table_description.append(&mut extra); + let hdu = fptr.create_table("COMPONENTS", &table_description)?; + hdu.write_col(&mut fptr, "UNQ_SOURCE_ID", &unq_source_id)?; + hdu.write_col(&mut fptr, "NAME", &names)?; + hdu.write_col(&mut fptr, "RA", &ra_degrees)?; + hdu.write_col(&mut fptr, "DEC", &dec_degrees)?; + for (idx, freq) in flux_freqs.iter().enumerate() { + hdu.write_col(&mut fptr, &format!("INT_FLX{freq}"), &flux_lists[idx])?; + } + hdu.write_col(&mut fptr, "MAJOR_DC", &majors)?; + hdu.write_col(&mut fptr, "MINOR_DC", &minors)?; + hdu.write_col(&mut fptr, "PA_DC", &pas)?; + hdu.write_col(&mut fptr, "MOD_TYPE", &mod_types)?; + hdu.write_col(&mut fptr, "COMP_TYPE", &comp_types)?; + hdu.write_col(&mut fptr, "NORM_COMP_PL", &norm_comp_pls)?; + hdu.write_col(&mut fptr, "ALPHA_PL", &alpha_pls)?; + hdu.write_col(&mut fptr, "NORM_COMP_CPL", &norm_comp_cpls)?; + hdu.write_col(&mut fptr, "ALPHA_CPL", &alpha_cpls)?; + hdu.write_col(&mut fptr, "CURVE_CPL", &curve_cpls)?; + + let table_description = vec![ + ColumnDescription::new("NAME") + .with_type(ColumnDataType::String) + .that_repeats(max_src_name + 4) + .create()?, + ColumnDescription::new("N1") + .with_type(ColumnDataType::Int) + .create()?, + ColumnDescription::new("N2") + .with_type(ColumnDataType::Int) + .create()?, + ColumnDescription::new("COEFF") + .with_type(ColumnDataType::Double) + .create()?, + ]; + + let hdu = fptr.create_table("SHAPELETS", &table_description)?; + hdu.write_col(&mut fptr, "NAME", &shapelet_sources)?; + hdu.write_col(&mut fptr, "N1", &shapelet_n1s)?; + hdu.write_col(&mut fptr, "N2", &shapelet_n2s)?; + hdu.write_col(&mut fptr, "COEFF", &shapelet_coeff_values)?; + + Ok(()) +} diff --git a/src/srclist/general_tests.rs b/src/srclist/general_tests.rs index 6b830018..62c0c7b0 100644 --- a/src/srclist/general_tests.rs +++ b/src/srclist/general_tests.rs @@ -7,18 +7,26 @@ use std::{ fs::File, io::{BufReader, Cursor, Read}, + path::PathBuf, }; -use approx::assert_abs_diff_eq; +use approx::{abs_diff_eq, assert_abs_diff_eq}; +use fits::write_source_list_jack; use marlu::RADec; -use tempfile::NamedTempFile; +use tempfile::{NamedTempFile, TempDir}; use vec1::vec1; -use super::*; +use super::{fits::parse_source_list, *}; use crate::constants::DEFAULT_SPEC_INDEX; fn test_two_sources_lists_are_the_same(sl1: &SourceList, sl2: &SourceList) { - assert_eq!(sl1.len(), sl2.len()); + assert_eq!( + sl1.len(), + sl2.len(), + "length mismatch \n {:?}\n != {:?}", + sl1.iter().collect_vec(), + sl2.iter().collect_vec() + ); for ((sl1_name, s1), (sl2_name, s2)) in sl1.iter().zip(sl2.iter()) { assert_eq!(sl1_name, sl2_name); assert_eq!(s1.components.len(), s2.components.len()); @@ -107,10 +115,10 @@ fn test_two_sources_lists_are_the_same(sl1: &SourceList, sl2: &SourceList) { } FluxDensityType::PowerLaw { .. } => { - assert!(matches!( - s2_comp.flux_type, - FluxDensityType::PowerLaw { .. } - )); + assert!( + matches!(s2_comp.flux_type, FluxDensityType::PowerLaw { .. }), + "{sl1_name}: fdtype mismatch \n {s1_comp:?}\n != {s2_comp:?}" + ); match s2_comp.flux_type { FluxDensityType::PowerLaw { .. } => { // The parameters of the power law may not @@ -119,7 +127,10 @@ fn test_two_sources_lists_are_the_same(sl1: &SourceList, sl2: &SourceList) { let s1_fd = s1_comp.flux_type.estimate_at_freq(150e6); let s2_fd = s2_comp.flux_type.estimate_at_freq(150e6); assert_abs_diff_eq!(s1_fd.freq, s2_fd.freq, epsilon = 1e-10); - assert_abs_diff_eq!(s1_fd.i, s2_fd.i, epsilon = 1e-10); + assert!( + abs_diff_eq!(s1_fd.i, s2_fd.i, epsilon = 1e-10), + "{sl1_name}: i flux mismatch\n {s1_comp:?}\n != {s2_comp:?}\n {s1_fd:?}\n != {s2_fd:?}" + ); assert_abs_diff_eq!(s1_fd.q, s2_fd.q, epsilon = 1e-10); assert_abs_diff_eq!(s1_fd.u, s2_fd.u, epsilon = 1e-10); assert_abs_diff_eq!(s1_fd.v, s2_fd.v, epsilon = 1e-10); @@ -136,10 +147,10 @@ fn test_two_sources_lists_are_the_same(sl1: &SourceList, sl2: &SourceList) { } FluxDensityType::CurvedPowerLaw { .. } => { - assert!(matches!( - s2_comp.flux_type, - FluxDensityType::PowerLaw { .. } - )); + assert!( + matches!(s2_comp.flux_type, FluxDensityType::CurvedPowerLaw { .. }), + "{sl1_name}: fdtype mismatch {s1_comp:?} {s2_comp:?}" + ); match s2_comp.flux_type { FluxDensityType::CurvedPowerLaw { .. } => { // The parameters of the curved power law may @@ -873,3 +884,233 @@ fn read_invalid_json_file() { )); } } + +fn get_fits_expected_srclist( + ref_freq: f64, + incl_list: bool, + incl_cpl: bool, + incl_shape: bool, +) -> SourceList { + let mut expected_srclist = SourceList::new(); + let cmp_type_gaussian = ComponentType::Gaussian { + maj: 20.0_f64.to_radians(), + min: 10.0_f64.to_radians(), + pa: 75.0_f64.to_radians(), + }; + #[rustfmt::skip] + let flux_type_list = FluxDensityType::List(vec1![ + FluxDensity { freq: 100e6, i: 3.0, q: 0.0, u: 0.0, v: 0.0,}, + FluxDensity { freq: 150e6, i: 2.0, q: 0.0, u: 0.0, v: 0.0,}, + FluxDensity { freq: 200e6, i: 1.0, q: 0.0, u: 0.0, v: 0.0,}, + ]); + #[rustfmt::skip] + let flux_type_pl = FluxDensityType::PowerLaw { + si: -0.8, + fd: FluxDensity { + freq: ref_freq, i: 2.0, q: 0.0, u: 0.0, v: 0.0, + }, + }; + #[rustfmt::skip] + let flux_type_cpl = FluxDensityType::CurvedPowerLaw { + si: -0.9, + fd: FluxDensity { + freq: ref_freq, i: 3.0, q: 0.0, u: 0.0, v: 0.0, + }, + q: 0.2, + }; + if incl_list { + expected_srclist.insert( + "point-list".into(), + Source { + components: vec![SourceComponent { + radec: RADec::from_degrees(0.0, 1.0), + comp_type: ComponentType::Point, + flux_type: flux_type_list.clone(), + }] + .into(), + }, + ); + } + expected_srclist.insert( + "point-pl".into(), + Source { + components: vec![SourceComponent { + radec: RADec::from_degrees(1.0, 2.0), + comp_type: ComponentType::Point, + flux_type: flux_type_pl.clone(), + }] + .into(), + }, + ); + if incl_cpl { + expected_srclist.insert( + "point-cpl".into(), + Source { + components: vec![SourceComponent { + radec: RADec::from_degrees(3.0, 4.0), + comp_type: ComponentType::Point, + flux_type: flux_type_cpl.clone(), + }] + .into(), + }, + ); + } + if incl_list { + expected_srclist.insert( + "gauss-list".into(), + Source { + components: vec![SourceComponent { + radec: RADec::from_degrees(0.0, 1.0), + comp_type: cmp_type_gaussian.clone(), + flux_type: flux_type_list, + }] + .into(), + }, + ); + } + expected_srclist.insert( + "gauss-pl".into(), + Source { + components: vec![SourceComponent { + radec: RADec::from_degrees(1.0, 2.0), + comp_type: cmp_type_gaussian.clone(), + flux_type: flux_type_pl.clone(), + }] + .into(), + }, + ); + if incl_cpl { + expected_srclist.insert( + "gauss-cpl".into(), + Source { + components: vec![SourceComponent { + radec: RADec::from_degrees(3.0, 4.0), + comp_type: cmp_type_gaussian, + flux_type: flux_type_cpl, + }] + .into(), + }, + ); + } + + if incl_shape { + #[rustfmt::skip] + expected_srclist.insert( + "shape-pl".into(), + Source { + components: vec![ + SourceComponent { + radec: RADec::from_degrees(1.0, 2.0), + comp_type: ComponentType::Shapelet { + maj: 20.0_f64.to_radians(), + min: 10.0_f64.to_radians(), + pa: 75.0_f64.to_radians(), + coeffs: vec![ + ShapeletCoeff { n1: 0, n2: 0, value: 0.9, }, + ShapeletCoeff { n1: 0, n2: 1, value: 0.2, }, + ShapeletCoeff { n1: 1, n2: 0, value: -0.2, }, + ] + .into(), + }, + flux_type: flux_type_pl.clone(), + }, + SourceComponent { + radec: RADec::from_degrees(1.0, 2.0), + comp_type: ComponentType::Shapelet { + maj: 20.0_f64.to_radians(), + min: 10.0_f64.to_radians(), + pa: 75.0_f64.to_radians(), + coeffs: vec![ + ShapeletCoeff { n1: 0, n2: 0, value: 0.8, }, + ] + .into(), + }, + flux_type: flux_type_pl.clone(), + }, + ] + .into(), + }, + ); + } + + expected_srclist +} + +#[test] +fn test_parse_gleam_fits() { + // TODO(Dev): CPL? + + // python -c 'from astropy.io import fits; import sys; from tabulate import tabulate; print(tabulate((i:=fits.open(sys.argv[-1])[1]).data, headers=[c.name for c in i.columns]))' /home/dev/src/hyperdrive_main/test_files/gleam.fits + // Name RAJ2000 DEJ2000 S_200 alpha beta a b pa + // ---------- --------- --------- ------- ------- ------ --- --- ---- + // point-pl 1 2 1 -0.8 0 0 0 0 + // gauss-pl 1 2 1 -0.8 0 20 10 75 + + let res_srclist = parse_source_list(&PathBuf::from("test_files/gleam.fits")).unwrap(); + let expected_srclist = get_fits_expected_srclist(200e6, false, true, false); + // dbg!(&res_srclist, &expected_srclist); + test_two_sources_lists_are_the_same(&res_srclist, &expected_srclist); +} + +#[test] +fn test_parse_jack_fits() { + // python -c 'from astropy.io import fits; import sys; from tabulate import tabulate; [print(tabulate(hdu.data, headers=[c.name for c in hdu.columns], tablefmt="plain")) for hdu in fits.open(sys.argv[-1])[1:]]' test_files/jack.fits + // UNQ_SOURCE_ID NAME RA DEC INT_FLX100 INT_FLX150 INT_FLX200 MAJOR_DC MINOR_DC PA_DC MOD_TYPE COMP_TYPE NORM_COMP_PL ALPHA_PL NORM_COMP_CPL ALPHA_CPL CURVE_CPL + // point-list point-list_C0 0 1 3 2 1 0 0 0 nan P 1 0 0 0 0 + // point-pl point-pl_C0 1 2 3.5 2.5 2 0 0 0 pl P 2 -0.8 0 0 0 + // point-cpl point-cpl_C0 3 4 5.6 3.8 3 0 0 0 cpl P 0 0 3 -0.9 0.2 + // gauss-list gauss-list_C0 0 1 3 2 1 20 10 75 nan G 1 0 0 0 0 + // gauss-pl gauss-pl_C0 1 2 3.5 2.5 2 20 10 75 pl G 2 -0.8 0 0 0 + // gauss-cpl gauss-cpl_C0 3 4 5.6 3.8 3 20 10 75 cpl G 0 0 3 -0.9 0.2 + // shape-pl shape-pl_C0 1 2 3.5 2.5 2 20 10 75 pl S 2 -0.8 0 0 0 + // shape-pl shape-pl_C1 1 2 3.5 2.5 2 20 10 75 pl S 2 -0.8 0 0 0 + // NAME N1 N2 COEFF + // shape-pl_C0 0 0 0.9 + // shape-pl_C0 0 1 0.2 + // shape-pl_C0 1 0 -0.2 + // shape-pl_C1 0 0 0.8 + + // setup logging + // use crate::cli::setup_logging; + // setup_logging(3).expect("Failed to setup logging"); + + let res_srclist = parse_source_list(&PathBuf::from("test_files/jack.fits")).unwrap(); + let expected_srclist = get_fits_expected_srclist(200e6, true, true, true); + // dbg!(&res_srclist, &expected_srclist); + test_two_sources_lists_are_the_same(&res_srclist, &expected_srclist); +} + +#[test] +fn test_parse_lobes_fits() { + // python -c 'from astropy.io import fits; import sys; from tabulate import tabulate; [print(tabulate(hdu.data, headers=[c.name for c in hdu.columns], tablefmt="plain")) for hdu in fits.open(sys.argv[-1])[1:]]' test_files/lobes.fits + // UNQ_SOURCE_ID NAME RA DEC INT_FLX100 INT_FLX150 INT_FLX200 MAJOR_DC MINOR_DC PA_DC MOD_TYPE COMP_TYPE NORM_COMP_PL ALPHA_PL NORM_COMP_CPL ALPHA_CPL CURVE_CPL + // point-list point-list_GID0 0 1 3 2 1 0 0 0 nan P 1 0 0 0 0 + // point-pl point-pl_GID0 1 2 3.5 2.5 2 0 0 0 pl P 2 -0.8 0 0 0 + // point-cpl point-cpl_GID0 3 4 5.6 3.8 3 0 0 0 cpl P 0 0 3 -0.9 0.2 + // gauss-list gauss-list_GID0 0 1 3 2 1 20 10 75 nan G 1 0 0 0 0 + // gauss-pl gauss-pl_GID0 1 2 3.5 2.5 2 20 10 75 pl G 2 -0.8 0 0 0 + // gauss-cpl gauss-cpl_GID0 3 4 5.6 3.8 3 20 10 75 cpl G 0 0 3 -0.9 0.2 + + // setup logging + // use crate::cli::setup_logging; + // setup_logging(3).expect("Failed to setup logging"); + + let res_srclist = parse_source_list(&PathBuf::from("test_files/lobes.fits")).unwrap(); + let expected_srclist = get_fits_expected_srclist(200e6, true, true, false); + // dbg!(&res_srclist, &expected_srclist); + test_two_sources_lists_are_the_same(&res_srclist, &expected_srclist); +} + +#[test] +fn test_convert_yaml_to_fits() { + let yaml_path = PathBuf::from("test_files/jack.yaml"); + let mut f = BufReader::new(File::open(yaml_path).unwrap()); + let expected_srclist = hyperdrive::source_list_from_yaml(&mut f).unwrap(); + + let tmp_dir = TempDir::new().unwrap(); + let fits_path = tmp_dir.path().join("test.fits"); + write_source_list_jack(&fits_path, &expected_srclist, None).unwrap(); + + let res_srclist = parse_source_list(&fits_path).unwrap(); + test_two_sources_lists_are_the_same(&res_srclist, &expected_srclist); +} diff --git a/src/srclist/mod.rs b/src/srclist/mod.rs index 1cb95503..311162f8 100644 --- a/src/srclist/mod.rs +++ b/src/srclist/mod.rs @@ -6,6 +6,7 @@ //! pub(crate) mod ao; +pub(crate) mod fits; pub(crate) mod hyperdrive; pub(crate) mod read; pub(crate) mod rts; @@ -41,6 +42,9 @@ pub(crate) enum SourceListType { #[strum(serialize = "hyperdrive")] Hyperdrive, + #[strum(serialize = "fits")] + Fits, + #[strum(serialize = "rts")] Rts, diff --git a/src/srclist/read.rs b/src/srclist/read.rs index 479fb2fa..869c422e 100644 --- a/src/srclist/read.rs +++ b/src/srclist/read.rs @@ -9,7 +9,7 @@ use std::{fs::File, path::Path}; use log::{debug, trace}; use super::{error::ReadSourceListError, SourceList, SourceListType}; -use crate::srclist::{ao, hyperdrive, rts, woden}; +use crate::srclist::{ao, fits, hyperdrive, rts, woden}; /// Given the path to a sky-model source list file (and optionally its type, /// e.g. "RTS style"), return a [SourceList] object. The [SourceListType] is @@ -23,7 +23,6 @@ pub(crate) fn read_source_list_file>( sl_type: Option, ) -> Result<(SourceList, SourceListType), ReadSourceListError> { debug!("Attempting to read source list"); - let mut f = std::io::BufReader::new(File::open(path)?); // If the file extension corresponds to YAML or JSON, we know what to // target. @@ -34,11 +33,13 @@ pub(crate) fn read_source_list_file>( match ext.as_deref() { Some("yaml" | "yml") => { debug!("Read as hyperdrive yaml"); + let mut f = std::io::BufReader::new(File::open(path)?); return hyperdrive::source_list_from_yaml(&mut f) .map(|r| (r, SourceListType::Hyperdrive)); } Some("json") => { debug!("Read as hyperdrive json"); + let mut f = std::io::BufReader::new(File::open(path)?); return hyperdrive::source_list_from_json(&mut f) .map(|r| (r, SourceListType::Hyperdrive)); } @@ -49,11 +50,13 @@ pub(crate) fn read_source_list_file>( match sl_type { Some(SourceListType::Hyperdrive) => { // Can be either yaml or json. + let mut f = std::io::BufReader::new(File::open(path)?); let json_err = match hyperdrive::source_list_from_json(&mut f) { Ok(sl) => return Ok((sl, SourceListType::Hyperdrive)), Err(e @ ReadSourceListError::Json(_)) => e.to_string(), Err(e) => return Err(e), }; + let mut f = std::io::BufReader::new(File::open(path)?); let yaml_err = match hyperdrive::source_list_from_yaml(&mut f) { Ok(sl) => return Ok((sl, SourceListType::Hyperdrive)), // If there was an error in parsing, then pass the parse @@ -64,23 +67,43 @@ pub(crate) fn read_source_list_file>( Err(ReadSourceListError::FailedToDeserialise { yaml_err, json_err }) } - Some(SourceListType::Rts) => match rts::parse_source_list(&mut f) { - Ok(sl) => Ok((sl, SourceListType::Rts)), - Err(e) => Err(e), + Some(SourceListType::Fits) => match fits::parse_source_list(path) { + Ok(sl) => Ok((sl, SourceListType::Fits)), + Err(e) => panic!("{:?}", e), }, - Some(SourceListType::AO) => match ao::parse_source_list(&mut f) { - Ok(sl) => Ok((sl, SourceListType::AO)), - Err(e) => Err(e), - }, + Some(SourceListType::Rts) => { + let mut f = std::io::BufReader::new(File::open(path)?); + match rts::parse_source_list(&mut f) { + Ok(sl) => Ok((sl, SourceListType::Rts)), + Err(e) => Err(e), + } + } - Some(SourceListType::Woden) => match woden::parse_source_list(&mut f) { - Ok(sl) => Ok((sl, SourceListType::Woden)), - Err(e) => Err(e), - }, + Some(SourceListType::AO) => { + let mut f = std::io::BufReader::new(File::open(path)?); + match ao::parse_source_list(&mut f) { + Ok(sl) => Ok((sl, SourceListType::AO)), + Err(e) => Err(e), + } + } + + Some(SourceListType::Woden) => { + let mut f = std::io::BufReader::new(File::open(path)?); + match woden::parse_source_list(&mut f) { + Ok(sl) => Ok((sl, SourceListType::Woden)), + Err(e) => Err(e), + } + } None => { // Try all kinds. + match fits::parse_source_list(path) { + Ok(sl) => return Ok((sl, SourceListType::Fits)), + Err(_) => trace!("Failed to read source list as a fits file"), + } + + let mut f = std::io::BufReader::new(File::open(path)?); match rts::parse_source_list(&mut f) { Ok(sl) => return Ok((sl, SourceListType::Rts)), Err(_) => { diff --git a/src/srclist/types/flux_density/mod.rs b/src/srclist/types/flux_density/mod.rs index 3ba05adf..74f4e73f 100644 --- a/src/srclist/types/flux_density/mod.rs +++ b/src/srclist/types/flux_density/mod.rs @@ -133,7 +133,7 @@ pub enum FluxDensityType { /// Similar to a power law. See Callingham et al. 2017, section 4.1. /// - /// S_\nu = a \nu^{\alpha} e^{q(\ln{v})^2} + /// S_\nu = a \nu^{\alpha} e^{q(\ln{\nu})^2} CurvedPowerLaw { /// Spectral index (alpha) si: f64, diff --git a/src/srclist/write.rs b/src/srclist/write.rs index fb0cca63..e3c5914a 100644 --- a/src/srclist/write.rs +++ b/src/srclist/write.rs @@ -12,7 +12,7 @@ use std::{ use log::{info, trace}; use super::{ - ao, hyperdrive, rts, woden, HyperdriveFileType, SourceList, SourceListType, + ao, fits, hyperdrive, rts, woden, HyperdriveFileType, SourceList, SourceListType, WriteSourceListError, }; @@ -43,6 +43,7 @@ pub(crate) fn write_source_list( output_ext.unwrap_or("").to_string(), )) } + (SourceListType::Fits, _) => fits::write_source_list_jack(path, sl, num_sources)?, (SourceListType::Rts, _) => { rts::write_source_list(&mut f, sl, num_sources)?; info!("Wrote rts-style source list to {}", path.display()); diff --git a/test_files/gen_fits_srclists.py b/test_files/gen_fits_srclists.py new file mode 100755 index 00000000..756f7234 --- /dev/null +++ b/test_files/gen_fits_srclists.py @@ -0,0 +1,177 @@ +#!/usr/bin/env python + +import pandas as pd +from astropy.coordinates import Angle +import astropy.units as u +from astropy.io import fits +from astropy.table import Table +import numpy as np +import re +import io + +arrays_jack = { + "UNQ_SOURCE_ID": [], + "NAME": [], + "RA": [], + "DEC": [], + "INT_FLX100": [], + "INT_FLX150": [], + "INT_FLX200": [], + "MAJOR_DC": [], + "MINOR_DC": [], + "PA_DC": [], + "MOD_TYPE": [], + "COMP_TYPE": [], + "NORM_COMP_PL": [], + "ALPHA_PL": [], + "NORM_COMP_CPL": [], + "ALPHA_CPL": [], + "CURVE_CPL": [], +} +arrays_jack_2 = { + "NAME": [], + "N1": [], + "N2": [], + "COEFF": [], +} +arrays_lobes = { + "UNQ_SOURCE_ID": [], + "NAME": [], + "RA": [], + "DEC": [], + "INT_FLX100": [], + "INT_FLX150": [], + "INT_FLX200": [], + "MAJOR_DC": [], + "MINOR_DC": [], + "PA_DC": [], + "MOD_TYPE": [], + "COMP_TYPE": [], + "NORM_COMP_PL": [], + "ALPHA_PL": [], + "NORM_COMP_CPL": [], + "ALPHA_CPL": [], + "CURVE_CPL": [], +} +arrays_gleam = { + "Name": [], + "RAJ2000": [], + "DEJ2000": [], + "S_200": [], + "alpha": [], + "beta": [], + "a": [], + "b": [], + "pa": [], +} +for (i,(name, cmp, ra, de, fd100, fd150, fd200, maj, min, pa, mtyp, ctyp, alpha, q, )) in enumerate([ + ["point-list", 0, 0., 1., 3.0, 2.0, 1., 0., 0., 0., "nan", "P", 0.0, 0.0 ], + ["point-pl", 0, 1., 2., 3.5, 2.5, 2., 0., 0., 0., "pl", "P", -0.8, 0.0 ], + ["point-cpl", 0, 3., 4., 5.6, 3.8, 3., 0., 0., 0., "cpl", "P", -0.9, 0.2 ], + ["gauss-list", 0, 0., 1., 3., 2., 1., 20., 10., 75., "nan", "G", 0.0, 0.0 ], + ["gauss-pl", 0, 1., 2., 3.5, 2.5, 2., 20., 10., 75., "pl", "G", -0.8, 0.0 ], + ["gauss-cpl", 0, 3., 4., 5.6, 3.8, 3., 20., 10., 75., "cpl", "G", -0.9, 0.2 ], + ["shape-pl", 0, 1., 2., 3.5, 2.5, 2., 20., 10., 75., "pl", "S", -0.8, 0.0 ], + ["shape-pl", 1, 1., 2., 3.5, 2.5, 2., 20., 10., 75., "pl", "S", -0.8, 0.0 ], +]): + + # i_ = f"{i:04d}" + arrays_jack["UNQ_SOURCE_ID"].append(f"{name}") + arrays_jack["NAME"].append(f"{name}_C{cmp}") + arrays_jack["RA"].append(ra) + arrays_jack["DEC"].append(de) + arrays_jack["INT_FLX100"].append(fd100) + arrays_jack["INT_FLX150"].append(fd150) + arrays_jack["INT_FLX200"].append(fd200) + arrays_jack["MAJOR_DC"].append(maj) + arrays_jack["MINOR_DC"].append(min) + arrays_jack["PA_DC"].append(pa) + arrays_jack["MOD_TYPE"].append(mtyp) + arrays_jack["COMP_TYPE"].append(ctyp) + + if mtyp == "cpl": + arrays_jack["NORM_COMP_PL"].append(0.0) + arrays_jack["ALPHA_PL"].append(0.0) + arrays_jack["NORM_COMP_CPL"].append(fd200) + arrays_jack["ALPHA_CPL"].append(alpha) + arrays_jack["CURVE_CPL"].append(q) + else: + arrays_jack["NORM_COMP_PL"].append(fd200) + arrays_jack["ALPHA_PL"].append(alpha) + arrays_jack["NORM_COMP_CPL"].append(0.0) + arrays_jack["ALPHA_CPL"].append(0.0) + arrays_jack["CURVE_CPL"].append(0.0) + + if ctyp == "S": # no shapelet support in gleam or lobes + continue + + arrays_lobes["UNQ_SOURCE_ID"].append(f"{name}") + arrays_lobes["NAME"].append(f"{name}_GID{cmp}") + arrays_lobes["RA"].append(ra) + arrays_lobes["DEC"].append(de) + arrays_lobes["INT_FLX100"].append(fd100) + arrays_lobes["INT_FLX150"].append(fd150) + arrays_lobes["INT_FLX200"].append(fd200) + arrays_lobes["MAJOR_DC"].append(maj) + arrays_lobes["MINOR_DC"].append(min) + arrays_lobes["PA_DC"].append(pa) + arrays_lobes["MOD_TYPE"].append(mtyp) + arrays_lobes["COMP_TYPE"].append(ctyp) + + if mtyp == "cpl": + arrays_lobes["NORM_COMP_PL"].append(0.0) + arrays_lobes["ALPHA_PL"].append(0.0) + arrays_lobes["NORM_COMP_CPL"].append(fd200) + arrays_lobes["ALPHA_CPL"].append(alpha) + arrays_lobes["CURVE_CPL"].append(q) + else: + arrays_lobes["NORM_COMP_PL"].append(fd200) + arrays_lobes["ALPHA_PL"].append(alpha) + arrays_lobes["NORM_COMP_CPL"].append(0.0) + arrays_lobes["ALPHA_CPL"].append(0.0) + arrays_lobes["CURVE_CPL"].append(0.0) + + if mtyp == "nan": # no list support in gleam + continue + + arrays_gleam['Name'].append(f"{name}") + arrays_gleam['RAJ2000'].append(ra) + arrays_gleam['DEJ2000'].append(de) + arrays_gleam['S_200'].append(fd200) + arrays_gleam['alpha'].append(alpha) + arrays_gleam['beta'].append(q) + arrays_gleam['a'].append(maj * 3600) + arrays_gleam['b'].append(min * 3600) + arrays_gleam['pa'].append(pa) + +for (i, (name, cmp, n1, n2, coeff)) in enumerate([ + ["shape-pl", 0, 0, 0, 0.9], + ["shape-pl", 0, 0, 1, 0.2], + ["shape-pl", 0, 1, 0, -0.2], + ["shape-pl", 1, 0, 0, 0.8], +]): + arrays_jack_2["NAME"].append(f"{name}_C{cmp}") + arrays_jack_2["N1"].append(n1) + arrays_jack_2["N2"].append(n2) + arrays_jack_2["COEFF"].append(coeff) + + +table = Table(arrays_jack) +table2 = Table(arrays_jack_2) +# convert table to fits hdus +hdus = fits.HDUList([fits.PrimaryHDU(), fits.table_to_hdu(table), fits.table_to_hdu(table2)]) +hdus.writeto('test_files/jack.fits', overwrite=True) +df = table.to_pandas() +print("\njack\n", df[[df.columns[0], *df.columns[2:]]].to_string(index=False)) +df = table2.to_pandas() +print("\njack2\n", df[[df.columns[0], *df.columns[1:]]].to_string(index=False)) + +table = Table(arrays_lobes) +table.write('test_files/lobes.fits', overwrite=True) +df = table.to_pandas() +print("\nlobes\n", df.to_string(index=False)) + +table = Table(arrays_gleam, ) +table.write('test_files/gleam.fits', overwrite=True) +df = table.to_pandas() +print("\ngleam\n", df.to_string(index=False)) \ No newline at end of file diff --git a/test_files/gleam.fits b/test_files/gleam.fits new file mode 100644 index 0000000000000000000000000000000000000000..ececb6a2e2ba304b4dfd3bcec6b0ae94bb39ec7e GIT binary patch literal 8640 zcmeI0&uiN-6vutsaesntK8>bnw>F1TCyO&YQYXY3WT&CjN@{RyV=J)yAG@soOm`h@ ze?*Uk9k%~qvYta6kV`_uVDR1{`oYUbUwx&+lLPN&Jah@UBATG&gxsddJk6H{5h)pX zG$4gwNyze0lP{-w9vyqr2HvIQu{>vM5;DPvSZ7>M+0OdisTathR$oiv0B1)5D3#s+oDRupR1CQ1c_v|Xzp1GKY4xj_*06KsUpabXtI`IEHP?nQ{H}bED z-Sd3f?G0VKwyP5|ixQTviD0vs6aHC1aav2+?f%pHxqW#3tm^Soef@lKCV4SW7UHg~ z$HPq}PGo*lC|+C7Bj*h~nco!Sx1$fftmF51FD&wOm643)Y{>=B_3wpu&+$%o{8KzV zWviABJa^a+^!hqGJh(FRh{s_wmC}!6SDlCL%YBp!)$iYS;XQkFFq+%|uhWJ99LFp5 zsq?ib-S_QQt7SkR@Z>8_#qs*i%KKvXT~qoho-uvrN?-6v+9#RE?YBH%;FsI+7CSE>$9Y+`)%fAp zk1~E