add contrast scaling option and remove optionality of specphot config

Author: mileslucas

src/vampires_dpp/__init__.py

@@ -1 +1 @@
-__version__ = "0.13.3"
+__version__ = "0.13.4"

src/vampires_dpp/cli/new.py

@@ -122,7 +122,7 @@ def get_base_settings(template: PipelineConfig) -> PipelineConfig:
 def get_calib_settings(template: PipelineConfig) -> PipelineConfig:
     click.secho("Frame Calibration", bold=True)
     readline.set_completer(pathCompleter)
-    calib_dir = click.prompt("Enter path to calibration files", default="")
+    calib_dir = click.prompt("Enter path to calibration files (or press enter to skip)", default="")
     readline.set_completer()
     if calib_dir != "":
         template.calibrate.calib_directory = Path(calib_dir)
@@ -305,16 +305,26 @@ def get_diff_image_config(template: PipelineConfig) -> PipelineConfig:
 def get_specphot_settings(template: PipelineConfig) -> PipelineConfig:
     click.secho("Spectrophotometric Calibration", bold=True)
     ## Specphot Cal
-    if click.confirm(
-        "Would you like to do flux calibration?", default=template.specphot is not None
-    ):
+    unit_choices = ["e-/s", "contrast", "Jy", "Jy/arcsec^2"]
+    readline.set_completer(createListCompleter(unit_choices))
+    unit = click.prompt(
+        "Choose output units",
+        type=click.Choice(unit_choices, case_sensitive=False),
+        default=template.specphot.unit,
+    )
+    readline.set_completer()
+    # default values
+    source = sptype = mag = mag_band = None
+    metric = "photometry"
+    if "Jy" in unit:
         readline.set_completer(pathCompleter)
         source = click.prompt(
-            ' - Enter source type ("pickles" or path to spectrum)', default="pickles"
+            ' - Enter spectrum source ("pickles" or path to spectrum)', default="pickles"
         )
         readline.set_completer()
         if source == "pickles":
             if template.target is not None:
+                click.echo("...Attempting to look up stellar flux from UCAC4/SIMBAD")
                 simbad_table = get_simbad_table(template.target.name)
                 sptype = re.match(r"\w\d[IV]{0,3}", simbad_table["SP_TYPE"][0]).group()
                 if len(sptype) == 2:
@@ -348,19 +358,15 @@ def get_specphot_settings(template: PipelineConfig) -> PipelineConfig:
                 default=mag_band,
                 type=click.Choice(list(FILTERS.keys()), case_sensitive=False),
             )
-        else:
-            sptype = mag = mag_band = None
-
+    if unit != "e-/s":
         metric = click.prompt(
             " - Select which metric to use for flux",
             default="photometry",
             type=click.Choice(["photometry", "sum"]),
         )
-        template.specphot = SpecphotConfig(
-            source=source, sptype=sptype, mag=mag, mag_band=mag_band, flux_metric=metric
-        )
-    else:
-        template.specphot = None
+    template.specphot = SpecphotConfig(
+        unit=unit, source=source, sptype=sptype, mag=mag, mag_band=mag_band, flux_metric=metric
+    )
 
     return template
 

src/vampires_dpp/constants.py

@@ -11,11 +11,6 @@
 SUBARU_LOC: Final[EarthLocation] = EarthLocation(lat=19.825504 * u.deg, lon=-155.4760187 * u.deg)
 
 
-SATSPOT_REF_NACT = 46.7
-# SATSPOT_REF_ANGLE = -83.2
-SATSPOT_REF_ANGLE = -88.2
-
-
 class InstrumentInfo(BaseModel):
     @property
     def pa_offset(self):

src/vampires_dpp/pipeline/config.py

@@ -128,13 +128,29 @@ class SpecphotConfig(BaseModel):
         Which frame analysis statistic to use for determining flux. "photometry" uses an aperture sum, while "sum" uses the sum in the analysis cutout window.
     """
 
-    source: Literal["pickles"] | Path = "pickles"
+    unit: Literal["e-/s", "contrast", "Jy", "Jy/arcsec^2"] = "e-/s"
+    source: Literal["pickles"] | Path | None = None
     sptype: str | None = None
     mag: float | None = None
-    mag_band: Literal["U", "B", "V", "r", "i", "J", "H", "K"] | None = "V"
-    unit: Literal["e-/s", "Jy", "Jy/arcsec^2"] = "Jy"
+    mag_band: Literal["U", "B", "V", "r", "i", "J", "H", "K"] | None = None
     flux_metric: Literal["photometry", "sum"] = "photometry"
 
+    def model_post_init(self, __context: Any) -> None:
+        if "Jy" in self.unit:
+            if self.source is None:
+                msg = "Must provide a spectrum or specify stellar model if you want to calibrate to Jy"
+                raise ValueError(msg)
+            if (
+                self.source == "pickles"
+                and self.sptype is None
+                or self.mag is None
+                or self.mag_band is None
+            ):
+                msg = "Must specify target magnitude (and filter) as well as spectral type to use 'pickles' stellar model"
+                raise ValueError(msg)
+
+        return super().model_post_init(__context)
+
 
 class CalibrateConfig(BaseModel):
     """Config for general image calibration.
@@ -444,7 +460,7 @@ class PipelineConfig(BaseModel):
     frame_select: FrameSelectConfig = FrameSelectConfig()
     align: AlignmentConfig = AlignmentConfig()
     coadd: CoaddConfig = CoaddConfig()
-    specphot: SpecphotConfig | None = None
+    specphot: SpecphotConfig = SpecphotConfig()
     diff_images: DiffImageConfig = DiffImageConfig()
     polarimetry: PolarimetryConfig | None = None
 

src/vampires_dpp/pipeline/pipeline.py

@@ -272,10 +272,9 @@ def process_group(self, group, group_key: str, output_path: Path):
         )
         logger.debug(f"Finished frame alignment for group {group_key}")
         ## Step 5: Spectrophotometric calibration
-        if self.config.specphot is not None:
-            logger.debug(f"Starting specphot cal for group {group_key}")
-            hdul = specphot_cal_hdul(hdul, metrics=metrics, config=self.config.specphot)
-            logger.debug(f"Finished specphot cal for group {group_key}")
+        logger.debug(f"Starting specphot cal for group {group_key}")
+        hdul = specphot_cal_hdul(hdul, metrics=metrics, config=self.config)
+        logger.debug(f"Finished specphot cal for group {group_key}")
         # Awkward: save registered data AFTER specphot calibration
         if self.config.align.save_intermediate and self.config.coadd.coadd:
             _, outpath = get_paths(output_path, output_directory=self.paths.aligned)

src/vampires_dpp/specphot/specphot.py

@@ -14,31 +14,43 @@
 SpecphotUnits: TypeAlias = Literal["e-/s", "Jy", "Jy/arcsec^2"]
 FluxMetric: TypeAlias = Literal["photometry", "sum"]
 
+
 # vampires_dpp/data
 SCEXAO_AREA: Final = 40.64 * u.m**2
 VEGASPEC: Final[SourceSpectrum] = SourceSpectrum.from_vega()
+SATSPOT_CONTRAST_POLY: Final = {  # Lucas et al. 2024 table 9
+    10.3: np.polynomial.Polynomial((0, -0.111, 14.217)),
+    15.5: np.polynomial.Polynomial((0, -0.043, 7.634)),
+    31.0: np.polynomial.Polynomial((0, -0.002, 0.561)),
+}
 
 
 def specphot_cal_hdul(hdul: fits.HDUList, metrics, config: SpecphotConfig):
     # determine any conversion factors
-    if config.unit in ("Jy", "Jy/arcsec^2"):
+    if config.specphot.unit in ("Jy", "Jy/arcsec^2"):
         assert metrics, "Must provide metrics to calculate photometry"
 
-    hdul, inst_flux = measure_inst_flux(hdul, metrics, config.flux_metric)
+    hdul, inst_flux = measure_inst_flux(
+        hdul, metrics, config.specphot.flux_metric, satspots=config.coronagraphic
+    )
 
-    match config.unit:
+    match config.specphot.unit:
         case "e-/s":
             conv_factor = 1
         case "Jy":
-            conv_factor = determine_jy_factor(hdul, inst_flux, config)
+            conv_factor = determine_jy_factor(hdul, inst_flux, config.specphot)
         case "Jy/arcsec^2":
-            conv_factor = determine_jy_factor(hdul, inst_flux, config) / hdul[0].header["PXAREA"]
+            conv_factor = (
+                determine_jy_factor(hdul, inst_flux, config.specphot) / hdul[0].header["PXAREA"]
+            )
+        case "contrast":
+            conv_factor = determine_contrast_factor(hdul, inst_flux)
 
     hdul[0].data *= conv_factor
     hdul["ERR"].data *= conv_factor
 
     info = fits.Header()
-    info["BUNIT"] = config.unit
+    info["BUNIT"] = config.specphot.unit
 
     for hdu in hdul:
         hdu.header.update(info)
@@ -50,7 +62,7 @@ def _format(number, sigfigs=4):
     return float(f"%.{sigfigs-1}g" % number)
 
 
-def measure_inst_flux(hdul, metrics, flux_metric: FluxMetric):
+def measure_inst_flux(hdul, metrics, flux_metric: FluxMetric, satspots: bool = False):
     info = fits.Header()
 
     match flux_metric:
@@ -62,13 +74,19 @@ def measure_inst_flux(hdul, metrics, flux_metric: FluxMetric):
     # collapse all but wavelength axis
     inst_flux = np.nanmedian(np.where(flux <= 0, np.nan, flux), axis=(1, 2))
     for flux, hdu in zip(inst_flux, hdul[2:], strict=True):
+        _, obs_filt = update_header_with_filt_info(hdu.header)
         field = hdu.header["FIELD"]
         info[f"hierarch DPP SPECPHOT INSTFLUX {field}"] = _format(flux), "[e-/s] Instrumental flux"
         info[f"hierarch DPP SPECPHOT INSTMAG {field}"] = (
             np.round(-2.5 * np.log10(flux), 3),
             "[mag] Instrumental magnitude",
         )
-
+        # add contrast
+        contrast = satellite_spot_contrast(hdu.header) if satspots else 1
+        info[f"hierarch DPP SPECPHOT CONTRAST {field}"] = (
+            contrast,
+            "Stellar flux to satspot flux ratio",
+        )
     for hdu in hdul:
         hdu.header.update(info)
     return hdul, inst_flux
@@ -176,3 +194,36 @@ def get_flux_from_metrics(metrics, config: SpecphotConfig) -> float:
             fluxes = metrics["sum"]
             weights = 1 / metrics["var"]
     return np.nansum(fluxes * weights) / np.nansum(weights)
+
+
+def satellite_spot_contrast(header: fits.Header) -> float:
+    # handle required headers
+    for key in ("X_GRDAMP", "X_GRDSEP", "WAVEAVE"):
+        if key not in header:
+            msg = f"Cannot calculate satspot flux ratio\n'{key}' was not found in header."
+            raise ValueError(msg)
+
+    amp = header["X_GRDAMP"]  # in um
+    sep = header["X_GRDSEP"]  # in lam/d
+    wave = header["WAVEAVE"]  # in nm
+    if sep not in SATSPOT_CONTRAST_POLY:
+        msg = f"No calibration data for astrogrid separation {sep} lam/D"
+        raise ValueError(msg)
+    poly = SATSPOT_CONTRAST_POLY[sep]
+    opd = amp * 1e3 / wave
+    contrast = poly(opd)  # satspot flux to stellar flux ratio
+    return contrast
+
+
+def determine_contrast_factor(hdul: fits.HDUList, inst_flux):
+    header = hdul[0].header
+    factors = []
+    for hdu in hdul[2:]:
+        field = hdu.header["FIELD"]
+        contrast = header[f"hierarch DPP SPECPHOT CONTRAST {field}"]
+        spotflux = header[f"hierarch DPP SPECPHOT INSTFLUX {field}"]
+
+        starflux = spotflux / contrast
+        factors.append(1 / starflux)
+
+    return np.array(factors)[None, :, None, None]