add m3 diat to polarimetry model footprint

.gitignore

@@ -166,4 +166,6 @@ docs/**/data
 docs/.jupyter_cache
 docs/jupyter_execute
 
-.DS_store
+.DS_store
+
+tests/data

pyproject.toml

@@ -25,10 +25,11 @@ dependencies = [
     "photutils",
     "pydantic",
     "reproject",
-    "synphot",
+    "result",
     "scikit_image",
     "scipy",
     "sep",
+    "synphot",
     "tomli_w",
     "tomli",
     "tqdm",

src/vampires_dpp/__init__.py

@@ -1 +1 @@
-__version__ = "0.12.9"
+__version__ = "0.12.10"

src/vampires_dpp/image_processing.py

@@ -164,7 +164,7 @@ def weighted_collapse(data: ArrayLike, angles: ArrayLike, **kwargs) -> NDArray:
     # if the variance is zero, return the mean
     if np.allclose(variance_frame, 0):
         derotated = derotate_cube(data, angles, **kwargs)
-        return np.nanmean(derotated, 0)
+        return bn.nanmean(derotated, 0)
 
     # expand the variance frame into a cube
     variance_cube = np.repeat(variance_frame, data.shape[0], axis=0)
@@ -173,8 +173,8 @@ def weighted_collapse(data: ArrayLike, angles: ArrayLike, **kwargs) -> NDArray:
     derotated_variance = derotate_cube(variance_cube, angles, **kwargs)
     derotated_variance[derotated_variance == 0] = np.inf
     # calculate weighted sum
-    numer = np.nansum(derotated_data / derotated_variance, axis=0)
-    denom = np.nansum(1 / derotated_variance, axis=0)
+    numer = bn.nansum(derotated_data / derotated_variance, axis=0)
+    denom = bn.nansum(1 / derotated_variance, axis=0)
     weighted_frame = numer / denom
     return weighted_frame
 
@@ -217,10 +217,14 @@ def collapse_cube(
                 weights = 1 / bn.nanvar(cube, axis=(1, 2), keepdims=True)
                 frame = bn.nansum(cube * weights, axis=0) / bn.nansum(weights)
         case "biweight":
-            frame = biweight_location(cube, axis=0, c=6, skip_nans=True)
+            # suppress all-nan axis warnings
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore")
+                frame = biweight_location(cube, axis=0, c=6, ignore_nan=True)
 
     if header is not None:
-        header["COL_METH"] = method, "DPP cube collapse method"
+        header["hierarch DPP COL METH"] = method, "DPP cube collapse method"
+        header["hierarch DPP NCOADD"] = cube.shape[0], "Num. of coadded frames in cube"
 
     return frame, header
 
@@ -505,3 +509,9 @@ def adaptive_sigma_clip_mask(data, sigma=10, boxsize=8):
             output_mask[inds] = np.abs(cutout - med) > sigma * std
 
     return output_mask
+
+
+def create_footprint(cube, angles):
+    mask = np.isfinite(cube)
+    derot = derotate_cube(mask.astype(float), angles)
+    return bn.nanmean(derot, axis=0)

src/vampires_dpp/indexing.py

@@ -80,7 +80,7 @@ def frame_angles_image(frame, center):
 def frame_angles_astro(frame, center):
     Ys, Xs = np.ogrid[0 : frame.shape[-2], 0 : frame.shape[-1]]
     # degrees East of North: phi = arctan(-x, y)
-    thetas = np.arctan2(center[-2] - Xs, Ys - center[-1])
+    thetas = np.arctan2(center[-1] - Xs, Ys - center[-2])
     return thetas
 
 

src/vampires_dpp/pdi/models.py

@@ -35,7 +35,8 @@ def load_calibration_file(header):
 
 class VAMPIRESMuellerMatrix(BaseModel):
     name: str = "ideal"
-    pbs_ratio: float = 1  # cam1 / cam2 ratio
+    m3_diat: float = 0
+    m3_offset: float = 0  # deg
     hwp_offset: float = 0  # deg
     hwp_phi: float = 0.5  # wave
     imr_offset: float = 0  # deg
@@ -49,7 +50,8 @@ class VAMPIRESMuellerMatrix(BaseModel):
     def common_path_mm(self, pa, alt, az, hwp, imr, hwp_adi_sync=True):
         # telescope
         pa_theta = np.deg2rad(pa)
-        m3 = mm.generic(epsilon=0, delta=np.pi)
+        m3_theta = np.deg2rad(self.m3_offset)
+        m3 = mm.generic(epsilon=self.m3_diat, theta=m3_theta, delta=np.pi)
         alt_theta = np.deg2rad(alt)
         tel_mm = mm.rotator(-alt_theta) @ m3 @ mm.rotator(pa_theta)
 
@@ -82,11 +84,9 @@ def __call__(self, flc_state: str, camera: int, *args, **kwargs) -> NDArray:
         flc_theta = np.deg2rad(self.flc_theta[flc_state])
         flc_mm = mm.waveplate(flc_theta, self.flc_phi * 2 * np.pi)
 
-        # beamsplitter
+        # beamsplitter - horiztonal to camera 1
         is_ordinary = camera == 1
         pbs_mm = mm.wollaston(is_ordinary)
-        if is_ordinary:
-            pbs_mm *= self.pbs_ratio
 
         M = pbs_mm @ flc_mm @ cp_mm
         return M.astype("f4")
@@ -111,7 +111,8 @@ def load_calib_file(cls, header):
         flc_theta = {k: t + table["flc_theta"] for k, t in zip(("A", "B"), (0, 45), strict=True)}
         return cls(
             name=table.name,
-            pbs_ratio=table["emgain"],
+            m3_diat=table["m3_diat"],
+            m3_offset=table["m3_theta"],
             hwp_offset=table["hwp_delta"],
             hwp_phi=table["hwp_phi"],
             imr_offset=table["imr_delta"],
@@ -171,6 +172,8 @@ def load_calib_file(cls, header):
         table = load_calibration_file(header)
         return cls(
             name=table.name,
+            m3_diat=table["m3_diat"],
+            m3_offset=table["m3_theta"],
             hwp_offset=table["hwp_delta"],
             hwp_phi=table["hwp_phi"],
             imr_offset=table["imr_delta"],

src/vampires_dpp/pdi/processing.py

@@ -343,15 +343,14 @@ def get_triplediff_set(table, row) -> dict | None:
             (table["RET-ANG1"] == key[0])
             & (table["U_FLC"] == key[1])
             & (table["U_CAMERA"] == key[2])
-            & (delta_angs < 2)  # 1 l/D rotation @ 45 l/d sep
+            & (delta_angs < 4)  # 2 l/D rotation @ 45 l/d sep
         )
         if np.any(mask):
             idx = delta_time[mask].argmin()
         else:
             return None
 
         output_set[key] = table.loc[mask, "path"].iloc[idx]
-
     return output_set
 
 
@@ -507,6 +506,8 @@ def polarization_ip_correct(stokes_data, phot_rad, method, header=None):
     if header is not None:
         header["IP_PQ"] = pQ, "I -> Q IP correction value"
         header["IP_PU"] = pU, "I -> U IP correction value"
+        header["IP_POL"] = np.hypot(pU, pQ) * 100, "[%] Residual IP"
+        header["IP_ANG"] = 0.5 * np.rad2deg(np.arctan2(pU, pQ)), "[deg] Residual IP angle"
         header["IP_METH"] = method, "IP measurement method"
     return stokes_data, header
 

src/vampires_dpp/pipeline/pipeline.py

@@ -460,7 +460,6 @@ def polarimetry_difference(self, table, method, num_proc=None, force=False):
                     full_paths.append(tuple(sorted(stokes_set.values())))
 
         full_path_set = list(set(full_paths))
-
         stokes_files = []
         # stokes_sets_path = self.paths.pdi_dir / f"{self.config.name}_stokes_sets.txt"
         # with stokes_sets_path.open("w") as fh:
@@ -523,7 +522,12 @@ def polarimetry_difference(self, table, method, num_proc=None, force=False):
         ## Collapse outputs
         stokes_data = np.array(stokes_data)
         stokes_err = np.array(stokes_err)
-        coll_frame, _ = collapse_frames(stokes_data)
+        coll_frame, _ = collapse_frames(np.nan_to_num(stokes_data))
+        footprint = np.mean(np.isfinite(stokes_data).astype("f4"), axis=0)
+        fits.writeto(
+            self.paths.pdi_dir / f"{self.config.name}_footprint.fits", footprint, overwrite=True
+        )
+        # coll_frame *= footprint
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
             coll_err = np.sqrt(np.nansum(stokes_err**2, axis=0)) / stokes_err.shape[0]

src/vampires_dpp/wcs.py

@@ -47,10 +47,15 @@ def get_coord_header(header, time=None):
 GAIA_CATALOGS = {"dr1": "I/337/gaia", "dr2": "I/345/gaia2", "dr3": "I/355/gaiadr3"}
 
 
-def get_gaia_astrometry(target, catalog="dr3", radius=1):
+def get_gaia_astrometry(target: str, catalog="dr3", radius=1):
+    """
+    Get coordinate from GAIA catalogue with proper motions and parallax
+
+    Radius is in arcminute.
+    """
     # get precise RA and DEC
     gaia_catalog_list = Vizier.query_object(
-        target, radius=radius * u.arcsec, catalog=GAIA_CATALOGS[catalog.lower()]
+        target, radius=radius * u.arcminute, catalog=GAIA_CATALOGS[catalog.lower()]
     )
     if len(gaia_catalog_list) == 0:
         return None
@@ -68,6 +73,7 @@ def get_gaia_astrometry(target, catalog="dr3", radius=1):
     return coord
 
 
-def get_precise_coord(coord, time, scale="utc"):
+def get_precise_coord(coord: SkyCoord, time: str, scale="utc"):
+    """Use astropy to get proper-motion corrected coordinate"""
     _time = Time(time, scale=scale, location=SUBARU_LOC)
     return coord.apply_space_motion(_time)