Merge pull request #10 from AlecThomson/two_factor

Factor of two option

README.md

@@ -12,7 +12,7 @@ ASKAP MSs are produced in a way that breaks compatibility with most other imager
 
 1. `fix_ms_dir` : ASKAP MeasurementSets are phased towards the centre of field, but not the centre of its given beam. This utility reads the appropriate offsets to the beam centre from the `BEAM_OFFSET` and updates the `FIELD` table, as well as the phase and delay reference columns.
 
-2. `fix_ms_corrs` : ASKAP MeasurementSets, as calibrated by the obervatory, provide correlations in the instrument frame. ASKAP has a unique 'roll' axis which means, in principle, the instrument frame can be at any arbitrary rotation on the sky. This utility applies the appropriate rotation matrix to the visibilities such the 'X' is aligned North-South and 'Y' is aligned East-West (IAU convention). Further, ASKAPsoft defines Stokes I as $I=XX+YY$, whereas most other telescopes use $I=\frac{1}{2}(XX+YY)$ (note this also applies to all other Stokes paramters). This factor is also corrected for here at the same time as the rotation.
+2. `fix_ms_corrs` : ASKAP MeasurementSets, as calibrated by the obervatory, provide correlations in the instrument frame. ASKAP has a unique 'roll' axis which means, in principle, the instrument frame can be at any arbitrary rotation on the sky. This utility applies the appropriate rotation matrix to the visibilities such the 'X' is aligned North-South and 'Y' is aligned East-West (IAU convention). Further, ASKAPsoft defines Stokes I as $I=XX+YY$, whereas most other telescopes use $I=\frac{1}{2}(XX+YY)$ (note this also applies to all other Stokes paramters). This factor is also corrected for here at the same time as the rotation. If you have calibrated with non-ASKAPsoft tools, you may need to use the `--no-fix-stokes-factor` option, which will disable the factor of two correction and just do a rotation.
 
 For convenience, we also provide `fix_ms` which does both of the above!
 

fixms/fix_ms.py

@@ -39,7 +39,12 @@ def get_parser() -> ArgumentParser:
         default="CORRECTED_DATA",
         help="The column to write the corrected data to",
     )
-
+    parser.add_argument(
+        "--no-fix-stokes-factor",
+        dest="no_fix_stokes_factor",
+        action="store_true",
+        help="Don't fix the Stokes factor. Use this if you have *not* used ASKAPsoft. If you have used ASKAPsoft, you should leave this option alone.",
+    )
     return parser
 
 
@@ -55,6 +60,7 @@ def cli() -> None:
         chunksize=args.chunksize,
         data_column=args.data_column,
         corrected_data_column=args.corrected_data_column,
+        fix_stokes_factor=not args.no_fix_stokes_factor,
     )
 
 

fixms/fix_ms_corrs.py

@@ -58,13 +58,16 @@ def get_pol_axis(ms: Path, feed_idx: Optional[int] = None) -> u.Quantity:
     return pol_ang
 
 
-def convert_correlations(correlations: np.ndarray, pol_axis: u.Quantity) -> np.ndarray:
+def convert_correlations(
+    correlations: np.ndarray, pol_axis: u.Quantity, fix_stokes_factor: bool = True
+) -> np.ndarray:
     """
     Convert ASKAP standard correlations to the 'standard' correlations
 
     Args:
         correlations (np.ndarray): The correlations from the MS. Has shape (NCOR, NCHAN, 4)
         pol_axis (astropy.units.Quantity): The polarization axis angle of the MS
+        fix_stokes_factor (bool, optional): Whether to fix the Stokes factor. Defaults to True.
 
     Returns:
         np.ndarray: The correlations in the 'standard' format
@@ -201,7 +204,12 @@ def convert_correlations(correlations: np.ndarray, pol_axis: u.Quantity) -> np.n
         ]
     )
     # This is a matrix multiplication broadcasted along the time and channel axes
-    return np.einsum("ij,...j->...i", correction_matrix, correlations)
+    rotated_correlations = np.einsum("ij,...j->...i", correction_matrix, correlations)
+    if not fix_stokes_factor:
+        logger.warning("Not fixing Stokes factor! Multiplying by 0.5...")
+        rotated_correlations *= 0.5
+
+    return rotated_correlations
 
 
 def get_data_chunk(
@@ -270,6 +278,7 @@ def fix_ms_corrs(
     chunksize: int = 10_000,
     data_column: str = "DATA",
     corrected_data_column: str = "CORRECTED_DATA",
+    fix_stokes_factor: bool = True,
 ) -> None:
     """Apply corrections to the ASKAP visibilities to bring them inline with
     what is expectede from other imagers, including CASA and WSClean. The
@@ -285,6 +294,7 @@ def fix_ms_corrs(
         chunksize (int, optional): Size of chunked data to correct. Defaults to 10_000.
         data_column (str, optional): The name of the data column to correct. Defaults to "DATA".
         corrected_data_column (str, optional): The name of the corrected data column. Defaults to "CORRECTED_DATA".
+        fix_stokes_factor (bool, optional): Whether to fix the Stokes factor. Defaults to True.
     """
     logger.info(f"Correcting {data_column} of {str(ms)}.")
 
@@ -358,6 +368,7 @@ def fix_ms_corrs(
                 data_chunk_cor = convert_correlations(
                     data_chunk,
                     pol_axis,
+                    fix_stokes_factor=fix_stokes_factor,
                 )
                 tab.putcol(
                     corrected_data_column,
@@ -396,12 +407,19 @@ def cli():
         default="CORRECTED_DATA",
         help="The column to write the corrected data to",
     )
+    parser.add_argument(
+        "--no-fix-stokes-factor",
+        dest="no_fix_stokes_factor",
+        action="store_true",
+        help="Don't fix the Stokes factor. Use this if you have *not* used ASKAPsoft. If you have used ASKAPsoft, you should leave this option alone.",
+    )
     args = parser.parse_args()
     fix_ms_corrs(
         Path(args.ms),
         chunksize=args.chunksize,
         data_column=args.data_column,
         corrected_data_column=args.corrected_data_column,
+        fix_stokes_factor=not args.no_fix_stokes_factor,
     )
 
 

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "fixms"
-version = "0.1.4"
+version = "0.1.5"
 description = ""
 authors = ["Alec Thomson (S&A, Kensington WA) <[email protected]>"]
 readme = "README.md"