diff --git a/src/vampires_dpp/__init__.py b/src/vampires_dpp/__init__.py index 26c36ca..f2b93b1 100644 --- a/src/vampires_dpp/__init__.py +++ b/src/vampires_dpp/__init__.py @@ -1 +1 @@ -__version__ = "0.13.3" +__version__ = "0.13.4" diff --git a/src/vampires_dpp/cli/new.py b/src/vampires_dpp/cli/new.py index b8f4b9c..044460a 100644 --- a/src/vampires_dpp/cli/new.py +++ b/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 diff --git a/src/vampires_dpp/constants.py b/src/vampires_dpp/constants.py index 3d3b359..4f6602d 100644 --- a/src/vampires_dpp/constants.py +++ b/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): diff --git a/src/vampires_dpp/pipeline/config.py b/src/vampires_dpp/pipeline/config.py index 071809f..2a811ea 100644 --- a/src/vampires_dpp/pipeline/config.py +++ b/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 diff --git a/src/vampires_dpp/pipeline/pipeline.py b/src/vampires_dpp/pipeline/pipeline.py index 9d9c196..3c70db5 100644 --- a/src/vampires_dpp/pipeline/pipeline.py +++ b/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) diff --git a/src/vampires_dpp/specphot/specphot.py b/src/vampires_dpp/specphot/specphot.py index 7b41af0..fd3fdee 100644 --- a/src/vampires_dpp/specphot/specphot.py +++ b/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]