calibration

Calibrate non-redundant-mask interferometric data with calibrator observations.

average_calib_files(list_nrm, sig_thres=2, display=False)

Average data extracted from multiple calibrator files.

Parameters:

Name	Type	Description	Default
list_nrm	list of object	Extracted NRM data for calibrator files.	required
sig_thres	float	Sigma-clipping threshold around the median.	2
display	bool	Whether to display diagnostic plots.	False

Returns:

Type	Description
Munch	Averaged calibrator observables and uncertainties.

Source code in src/amical/calibration.py

def average_calib_files(list_nrm, sig_thres=2, display=False):
    """Average data extracted from multiple calibrator files.

    Parameters
    ----------
    list_nrm : list of object
        Extracted NRM data for calibrator files.
    sig_thres : float, default=2
        Sigma-clipping threshold around the median.
    display : bool, default=False
        Whether to display diagnostic plots.

    Returns
    -------
    munch.Munch
        Averaged calibrator observables and uncertainties.
    """
    nfiles = len(list_nrm)
    l_pa = np.zeros(nfiles)
    cp_vs_file, e_cp_vs_file = [], []
    vis2_vs_file, e_vis2_vs_file = [], []

    # Fill array containing each vis2 and cp across files.
    for n in range(nfiles):
        nrm = list_nrm[n]
        # Use the header stored at extraction: the original FITS file may not
        # be available where the calibration runs.
        hdr = nrm.infos.get("hdr") or {}
        # todo: Check parallactic angle param of a real NIRISS header.
        l_pa[n] = hdr.get("PARANG", 0)

        cp = nrm.cp
        e_cp = nrm.e_cp
        vis2 = nrm.vis2
        e_vis2 = nrm.e_vis2

        cp_vs_file.append(cp)
        e_cp_vs_file.append(e_cp)
        vis2_vs_file.append(vis2)
        e_vis2_vs_file.append(e_vis2)

    bl = list_nrm[0].bl

    cp_vs_file = np.array(cp_vs_file)
    e_cp_vs_file = np.array(e_cp_vs_file)
    vis2_vs_file = np.array(vis2_vs_file)
    e_vis2_vs_file = np.array(e_vis2_vs_file)

    zero_uncer = e_vis2_vs_file == 0
    e_vis2_vs_file[zero_uncer] = np.max(e_vis2_vs_file)

    cmn_vis2, std_vis2 = wtmn(vis2_vs_file, e_vis2_vs_file)
    cmn_cp, std_cp = wtmn(cp_vs_file, e_cp_vs_file)

    # Apply sigma clipping on the averages
    cmn_vis2_clip, std_vis2_clip = _apply_sig_clip(
        vis2_vs_file, e_vis2_vs_file, sig_thres=sig_thres, var="V2", display=display
    )
    cmn_cp_clip, std_cp_clip = _apply_sig_clip(
        cp_vs_file,
        e_cp_vs_file,
        sig_thres=sig_thres,
        ymax=10,
        var="CP",
        display=display,
    )

    res = {
        "f_v2_clip": np.array(cmn_vis2_clip),
        "f_v2": np.array(cmn_vis2),
        "std_vis2_clip": np.array(std_vis2_clip),
        "std_vis2": np.array(std_vis2),
        "f_cp_clip": np.array(cmn_cp_clip),
        "f_cp": np.array(cmn_cp),
        "std_cp_clip": np.array(std_cp_clip),
        "std_cp": np.array(std_cp),
        "bl": bl,
        "pa": l_pa,
    }

    return dict2class(res)

calibrate(res_t, res_c, clip=False, sig_thres=2, apply_phscorr=False, apply_atmcorr=False, normalize_err_indep=False, display=False)

Calibrate squared visibilities and closure phases.

Parameters:

Name	Type	Description	Default
res_t	object	Extracted NRM data for the science target.	required
res_c	object or list of object	Extracted NRM data for one or more calibrators.	required
clip	bool	Whether to sigma-clip calibrator observables.	False
sig_thres	float	Sigma-clipping threshold around the median.	2
apply_phscorr	bool	Whether to apply the piston phasor correction. Not recommended: it is known to bias the calibrated V2 (see GitHub issue #347).	False
apply_atmcorr	bool	Whether to apply the seeing and windshake correction to squared visibilities.	False
normalize_err_indep	bool	Whether to normalize closure-phase uncertainties for independent closure phases.	False
display	bool	Whether to display figures.	False

Returns:

Type	Description
Munch	Calibrated data including squared visibilities, closure phases, visibility amplitudes and phases, coordinates, wavelength, and raw target and calibrator inputs.

Source code in src/amical/calibration.py

def calibrate(
    res_t,
    res_c,
    clip=False,
    sig_thres=2,
    apply_phscorr=False,
    apply_atmcorr=False,
    normalize_err_indep=False,
    display=False,
):
    """Calibrate squared visibilities and closure phases.

    Parameters
    ----------
    res_t : object
        Extracted NRM data for the science target.
    res_c : object or list of object
        Extracted NRM data for one or more calibrators.
    clip : bool, default=False
        Whether to sigma-clip calibrator observables.
    sig_thres : float, default=2
        Sigma-clipping threshold around the median.
    apply_phscorr : bool, default=False
        Whether to apply the piston phasor correction. Not recommended: it is
        known to bias the calibrated V2 (see GitHub issue #347).
    apply_atmcorr : bool, default=False
        Whether to apply the seeing and windshake correction to squared
        visibilities.
    normalize_err_indep : bool, default=False
        Whether to normalize closure-phase uncertainties for independent
        closure phases.
    display : bool, default=False
        Whether to display figures.

    Returns
    -------
    munch.Munch
        Calibrated data including squared visibilities, closure phases,
        visibility amplitudes and phases, coordinates, wavelength, and raw
        target and calibrator inputs.
    """

    if not isinstance(res_c, list):
        res_c = [res_c]

    calib_tab = average_calib_files(res_c, sig_thres=sig_thres, display=display)

    if clip:
        cmn_v2_c, cmn_cp_c, std_v2_c, std_cp_c = (
            calib_tab.f_v2_clip,
            calib_tab.f_cp_clip,
            calib_tab.std_vis2_clip,
            calib_tab.std_cp_clip,
        )
    else:
        cmn_v2_c, cmn_cp_c, std_v2_c, std_cp_c = (
            calib_tab.f_v2,
            calib_tab.f_cp,
            calib_tab.std_vis2,
            calib_tab.std_cp,
        )

    v2_corr_t = 1
    if apply_atmcorr:
        v2_corr_t = _calc_correction_atm_vis2(res_t)

    if apply_phscorr:
        warnings.warn(
            "apply_phscorr is known to over-correct V2: on simulated data the "
            "calibrated V2 is worse than with no correction. See "
            "https://github.com/SAIL-Labs/AMICAL/issues/347",
            UserWarning,
            stacklevel=2,
        )
        v2_corr_t *= res_t.matrix.phs_v2corr

    # Raw V2 target (corrected from atm correction and phasors.)
    v2_t = res_t.vis2 / v2_corr_t
    e_v2_t = res_t.e_vis2 / v2_corr_t
    #  e_v2_c = res_c[0].e_vis2/v2_corr_t

    # Raw CP target
    cp_t = res_t.cp
    e_cp_t = res_t.e_cp

    # Calibration by the weighted averages and taking into accound the std of the calibrators.
    # ---------------------------------------------
    vis2_calib = v2_t / cmn_v2_c
    cp_calib = cp_t - cmn_cp_c

    #
    n_holes = res_t.mask.n_holes
    if normalize_err_indep:
        err_scale = np.sqrt(n_holes / 3.0)
    else:
        err_scale = 1

    # Quadratic added error due to calibrator dispersion (the average is weightened (see wtmn from amical.tools)).
    weightened_error = True
    if weightened_error:
        e_vis2_calib = np.sqrt(
            e_v2_t**2 / cmn_v2_c**2 + std_v2_c**2 * v2_t**2 / cmn_v2_c**4
        )
    else:
        e_vis2_calib = np.sqrt(e_v2_t**2 + std_v2_c**2)

    e_cp_calib = np.sqrt(e_cp_t**2 + std_cp_c**2) * err_scale

    u1 = res_t.u[res_t.mask.bs2bl_ix[0, :]]
    v1 = res_t.v[res_t.mask.bs2bl_ix[0, :]]
    u2 = res_t.u[res_t.mask.bs2bl_ix[1, :]]
    v2 = res_t.v[res_t.mask.bs2bl_ix[1, :]]

    cal = {
        "vis2": vis2_calib,
        "e_vis2": e_vis2_calib,
        "cp": cp_calib,
        "e_cp": e_cp_calib,
        "u": res_t.u,
        "v": res_t.v,
        "wl": res_t.wl,
        "u1": u1,
        "v1": v1,
        "u2": u2,
        "v2": v2,
        "raw_t": res_t,
        "raw_c": res_c,
    }
    return dict2class(cal)