tools

General utilities for aperture-masking data processing and analysis.

apply_windowing(img, window=80.0, m=3.0)

Function for applying a super-Gaussian window to an image.

Parameters:

Name	Type	Description	Default
img	ndarray	2D array with the input image.	required
window	float	Half width at half maximum (HWHM) of the window function (default: 80.0).	80.0
m	float	Exponent used for the super-Gaussian function (default: 3.0).	3.0

Returns:

Type	Description
ndarray	2D array with the windowed input image.

Source code in src/amical/tools.py

def apply_windowing(
    img: np.ndarray, window: float = 80.0, m: float = 3.0
) -> np.ndarray:
    """
    Function for applying a super-Gaussian window to an image.

    Parameters
    ----------
    img : np.ndarray
        2D array with the input image.
    window : float
        Half width at half maximum (HWHM) of the window function
        (default: 80.0).
    m : float
        Exponent used for the super-Gaussian function (default: 3.0).

    Returns
    -------
    np.ndarray
        2D array with the windowed input image.
    """
    super_gauss = _window_map(len(img), float(window), float(m))

    # Apply the windowing
    return img * super_gauss

check_seeing_cond(list_nrm)

Extract seeing conditions and observables from NRM results.

Parameters:

Name	Type	Description	Default
list_nrm	list of object	Results from the extraction pipeline.	required

Returns:

Type	Description
Munch	Seeing, parallactic-angle, visibility, and closure-phase information.

Source code in src/amical/tools.py

def check_seeing_cond(list_nrm):  # pragma: no cover
    """Extract seeing conditions and observables from NRM results.

    Parameters
    ----------
    list_nrm : list of object
        Results from the extraction pipeline.

    Returns
    -------
    munch.Munch
        Seeing, parallactic-angle, visibility, and closure-phase information.
    """
    from astropy.io import fits

    l_seeing, l_vis2, l_cp, l_pa, l_mjd = [], [], [], [], []

    with fits.open(list_nrm[0].infos.filename) as fd:
        hdr = fd[0].header
    for nrm in list_nrm:
        with fits.open(nrm.infos.filename) as fd:
            hdr = fd[0].header
        pa = np.mean(sphere_parang(hdr))
        seeing = nrm.infos.seeing
        mjd = hdr["MJD-OBS"]
        l_vis2.append(np.mean(nrm.vis2))
        l_cp.append(np.mean(nrm.cp))
        l_seeing.append(seeing)
        l_pa.append(pa)
        l_mjd.append(mjd)

    res = {
        "pa": l_pa,
        "seeing": l_seeing,
        "vis2": l_vis2,
        "cp": l_cp,
        "mjd": l_mjd,
        "target": hdr["OBJECT"],
    }

    return dict2class(sanitize_array(res))

cov2cor(cov)

Convert a covariance matrix to a correlation matrix.

Parameters:

Name	Type	Description	Default
cov	ndarray	Square covariance matrix.	required

Returns:

Name	Type	Description
cor	ndarray	Correlation matrix.
sigma	ndarray	Standard deviations.

Source code in src/amical/tools.py

def cov2cor(cov):
    """Convert a covariance matrix to a correlation matrix.

    Parameters
    ----------
    cov : numpy.ndarray
        Square covariance matrix.

    Returns
    -------
    cor : numpy.ndarray
        Correlation matrix.
    sigma : numpy.ndarray
        Standard deviations.
    """
    diag = np.diag(cov)
    negative = np.flatnonzero(diag < 0.0)
    if negative.size:
        ix = negative[0]
        raise ValueError(f"diagonal cov[{ix},{ix}]={diag[ix]:e} is not positive")

    sigma = np.sqrt(diag)
    # cor[i, j] = cov[i, j] / sqrt(cov[i, i] * cov[j, j])
    cor = cov / np.sqrt(np.outer(diag, diag))

    return cor, sigma

crop_max(img, dim, offx=0, offy=0, filtmed=True, f=3)

Crop an image around its brightest pixel.

Parameters:

Name	Type	Description	Default
img	ndarray	Input image.	required
dim	int	Output dimension.	required
offx	int	Crop-center offsets.	0
offy	int	Crop-center offsets.	0
filtmed	bool	Whether to median-filter before locating the maximum.	True
f	int	Median-filter kernel size.	3

Returns:

Name	Type	Description
image	ndarray	Cropped image.
position	tuple of int	Brightest-pixel coordinates.

Source code in src/amical/tools.py

def crop_max(img, dim, offx=0, offy=0, filtmed=True, f=3):
    """Crop an image around its brightest pixel.

    Parameters
    ----------
    img : numpy.ndarray
        Input image.
    dim : int
        Output dimension.
    offx, offy : int, default=0
        Crop-center offsets.
    filtmed : bool, default=True
        Whether to median-filter before locating the maximum.
    f : int, default=3
        Median-filter kernel size.

    Returns
    -------
    image : numpy.ndarray
        Cropped image.
    position : tuple of int
        Brightest-pixel coordinates.
    """
    from astropy.nddata import Cutout2D

    xmax, ymax = find_max(img, filtmed=filtmed, f=f)

    X = xmax + offx
    Y = ymax + offy
    isz_max = 2 * np.min([X, img.shape[1] - X - 1, Y, img.shape[0] - Y - 1]) + 1
    if isz_max < dim:
        size_msg = (
            f"The specified cropped image size, {dim}, is greater than the distance to"
            " the PSF center in at least one dimension. The max size for this image is"
            f" {isz_max}"
        )
        raise ValueError(size_msg)
    cutout = Cutout2D(img, (X, Y), dim)
    return cutout.data, (X, Y)

eris_parang(hdr, n_dit=None)

Reads the header and creates an array giving the paralactic angle for each frame, taking into account the inital derotator position. The parallactic angles are linearly interpolated between the start and end of the template. If some frames are removed from the cube beforehand but the header is kept the same, then this approximation becomes even worse. TO-DO: Ideally, the indices of the frames that are removed from the cube should be passed as argument, and the true parallactic angles between START and END should be calculated instead of linearly interpolated.

Source code in src/amical/tools.py

def eris_parang(hdr, n_dit=None):
    """
    Reads the header and creates an array giving the paralactic angle for each frame,
    taking into account the inital derotator position.
    The parallactic angles are linearly interpolated between the start and end of the template.
    If some frames are removed from the cube beforehand but the header is kept the same,
    then this approximation becomes even worse.
    TO-DO:
    Ideally, the indices of the frames that are removed from the cube should be passed as argument,
    and the true parallactic angles between START and END should be calculated instead of linearly interpolated.
    """

    if n_dit is None:
        n_frames = hdr["NAXIS3"]
    else:
        n_frames = n_dit

    pupil_pos = hdr["HIERARCH ESO ADA PUPILPOS"]
    par_ang_start = hdr["HIERARCH ESO TEL PARANG START"]
    par_ang_end = hdr["HIERARCH ESO TEL PARANG END"]
    par_angles = np.linspace(par_ang_start, par_ang_end, n_frames)

    pos_angles = pupil_pos - par_angles + 2

    return pos_angles

find_max(img, filtmed=True, f=3)

Find the brightest pixel in an image.

Parameters:

Name	Type	Description	Default
img	ndarray	Input image.	required
filtmed	bool	Whether to median-filter before locating the maximum.	True
f	int	Median-filter kernel size.	3

Returns:

Type	Description
x, y : int	Brightest-pixel coordinates.

Source code in src/amical/tools.py

def find_max(img, filtmed=True, f=3):
    """Find the brightest pixel in an image.

    Parameters
    ----------
    img : numpy.ndarray
        Input image.
    filtmed : bool, default=True
        Whether to median-filter before locating the maximum.
    f : int, default=3
        Median-filter kernel size.

    Returns
    -------
    x, y : int
        Brightest-pixel coordinates.
    """
    from scipy.signal import medfilt2d

    if filtmed:
        fast = _find_max_median_fast(img, f)
        if fast is not None:
            return fast
        try:
            im_med = medfilt2d(img, f)
        except ValueError:
            img = img.astype(float)
            im_med = medfilt2d(img, f)
    else:
        im_med = img.copy()

    pos_max = np.where(im_med == im_med.max())

    X = pos_max[1][0]
    Y = pos_max[0][0]

    return X, Y

gauss_2d_asym(X, param)

Create an oriented two-dimensional Gaussian on an asymmetric grid.

Parameters:

Name	Type	Description	Default
X	sequence	X coordinates, y coordinates, and pixel scale in milliarcseconds.	required
param	dict	Gaussian amplitude, offsets, widths, and orientation.	required

Returns:

Type	Description
ndarray	Gaussian image.

Source code in src/amical/tools.py

def gauss_2d_asym(X, param):
    """Create an oriented two-dimensional Gaussian on an asymmetric grid.

    Parameters
    ----------
    X : sequence
        X coordinates, y coordinates, and pixel scale in milliarcseconds.
    param : dict
        Gaussian amplitude, offsets, widths, and orientation.

    Returns
    -------
    numpy.ndarray
        Gaussian image.
    """

    x_1d = X[0]
    y_1d = X[1]
    pixel_scale = X[2]

    dim = len(x_1d)

    x, y = np.meshgrid(x_1d, y_1d)

    fwhmx = param["fwhm_x"] / pixel_scale
    fwhmy = param["fwhm_y"] / pixel_scale

    sigma_x = fwhmx / np.sqrt(8 * np.log(2))
    sigma_y = fwhmy / np.sqrt(8 * np.log(2))

    amplitude = param["A"]
    x0 = dim // 2 + param["x0"] / pixel_scale
    y0 = dim // 2 + param["y0"] / pixel_scale
    theta = np.deg2rad(param["theta"])
    size_x = len(x)
    size_y = len(y)
    im = np.zeros([size_y, size_x])
    x0 = float(x0)
    y0 = float(y0)
    a = (np.cos(theta) ** 2) / (2 * sigma_x**2) + (np.sin(theta) ** 2) / (
        2 * sigma_y**2
    )
    b = -(np.sin(2 * theta)) / (4 * sigma_x**2) + (np.sin(2 * theta)) / (4 * sigma_y**2)
    c = (np.sin(theta) ** 2) / (2 * sigma_x**2) + (np.cos(theta) ** 2) / (
        2 * sigma_y**2
    )
    im = amplitude * np.exp(
        -(a * ((x - x0) ** 2) + 2 * b * (x - x0) * (y - y0) + c * ((y - y0) ** 2))
    )
    return im

jd2lst(lng, jd)

Convert a Julian date to local sidereal time.

Parameters:

Name	Type	Description	Default
lng	float	Longitude.	required
jd	float	Julian date.	required

Returns:

Type	Description
float	Local sidereal time.

Source code in src/amical/tools.py

def jd2lst(lng, jd):
    """Convert a Julian date to local sidereal time.

    Parameters
    ----------
    lng : float
        Longitude.
    jd : float
        Julian date.

    Returns
    -------
    float
        Local sidereal time.
    """
    c = [280.46061837, 360.98564736629, 0.000387933, 38710000.0]
    jd2000 = 2451545.0
    t0 = jd - jd2000
    t = t0 / 36525.0

    # Compute GST in seconds.
    theta = c[0] + (c[1] * t0) + t**2 * (c[2] - t / c[3])

    # Compute LST in hours.
    lst = (theta + lng) / 15.0
    cond_neg = lst < 0.0
    n = lst[cond_neg].size
    if n > 0:
        lst[cond_neg] = 24.0 + (lst[cond_neg] % 24)
    lst = lst % 24
    return lst

linear(x, param)

Evaluate the linear model used by the fitting utilities.

Parameters:

Name	Type	Description	Default
x	array_like	Independent-variable values.	required
param	dict[str, float]	Linear-model coefficients with "a" and "b" keys.	required

Returns:

Type	Description
array_like	Values of a * x + b.

Source code in src/amical/tools.py

def linear(x, param):
    """Evaluate the linear model used by the fitting utilities.

    Parameters
    ----------
    x : array_like
        Independent-variable values.
    param : dict[str, float]
        Linear-model coefficients with ``"a"`` and ``"b"`` keys.

    Returns
    -------
    array_like
        Values of ``a * x + b``.
    """
    a = param["a"]
    b = param["b"]
    y = a * x + b
    return y

load_bs_hdf5(filename)

Load extraction results from an HDF5 file.

Parameters:

Name	Type	Description	Default
filename	str or path - like	HDF5 filename.	required

Returns:

Type	Description
Munch	Loaded extraction result.

Source code in src/amical/tools.py

def load_bs_hdf5(filename):
    """Load extraction results from an HDF5 file.

    Parameters
    ----------
    filename : str or path-like
        HDF5 filename.

    Returns
    -------
    munch.Munch
        Loaded extraction result.
    """
    import h5py

    dict_bs = {"matrix": {}, "infos": {"hdr": {}}, "mask": {}}
    with h5py.File(filename, "r") as hf2:
        obs = hf2["obs"]
        for o in obs:
            dict_bs[o] = np.array(obs.get(o))

        matrix = hf2["matrix"]
        for key in matrix:
            dict_bs["matrix"][key] = np.array(matrix.get(key))

        if len(dict_bs["matrix"]["cp_cov"]) == 1:
            dict_bs["matrix"]["cp_cov"] = None

        mask = hf2["mask"]
        for key in mask:
            if key not in ["u1", "u2", "v1", "v2"]:
                dict_bs["mask"][key] = np.array(mask.get(key))

        t3_coord = {
            "u1": np.array(mask.get("u1")),
            "u2": np.array(mask.get("u2")),
            "v1": np.array(mask.get("v1")),
            "v2": np.array(mask.get("v2")),
        }

        dict_bs["mask"]["t3_coord"] = t3_coord

        infos = hf2["infos"]
        for key in hf2["infos"].attrs.keys():
            dict_bs["infos"][key] = infos.attrs[key]

        hdr = hf2["hdr"]
        for key in hdr.attrs.keys():
            dict_bs["infos"]["hdr"][key] = hdr.attrs[key]

        bs_save = dict2class(dict_bs)
    return bs_save

mas2rad(mas)

Convert milliarcseconds to radians.

Parameters:

Name	Type	Description	Default
mas	float or array_like	Angle in milliarcseconds.	required

Returns:

Type	Description
float or array_like	Angle in radians.

Source code in src/amical/tools.py

def mas2rad(mas):
    """Convert milliarcseconds to radians.

    Parameters
    ----------
    mas : float or array_like
        Angle in milliarcseconds.

    Returns
    -------
    float or array_like
        Angle in radians.
    """
    rad = mas * (10 ** (-3)) / (3600 * 180 / np.pi)
    return rad

norm_max(tab)

Normalize values by their maximum.

Parameters:

Name	Type	Description	Default
tab	ndarray or list	Input values.	required

Returns:

Type	Description
ndarray or list	Normalized values.

Source code in src/amical/tools.py

def norm_max(tab):
    """Normalize values by their maximum.

    Parameters
    ----------
    tab : numpy.ndarray or list
        Input values.

    Returns
    -------
    numpy.ndarray or list
        Normalized values.
    """
    tab_norm = tab / np.max(tab)
    return tab_norm

plot_circle(d, x, y, hole_radius, sz=1, display=True)

Create an image containing a circular disk.

Parameters:

Name	Type	Description	Default
d	int	Image dimension.	required
x	float	Disk-center coordinates.	required
y	float	Disk-center coordinates.	required
hole_radius	float	Disk radius.	required
sz	float	Disk value.	1
display	bool	Whether to display the image.	True

Returns:

Type	Description
ndarray	Disk image.

Source code in src/amical/tools.py

def plot_circle(d, x, y, hole_radius, sz=1, display=True):
    """Create an image containing a circular disk.

    Parameters
    ----------
    d : int
        Image dimension.
    x, y : float
        Disk-center coordinates.
    hole_radius : float
        Disk radius.
    sz : float, default=1
        Disk value.
    display : bool, default=True
        Whether to display the image.

    Returns
    -------
    numpy.ndarray
        Disk image.
    """
    chipsz = np.shape(d)[0]

    im = np.zeros([chipsz, chipsz])
    info = [len(im.shape), im.shape[0], im.shape[1], 3, len(im.ravel())]

    if isinstance(x, float | int):
        n_circ = 1
        xx = [x]
        yy = [y]
    else:
        n_circ = len(x)
        xx = x
        yy = y

    r = hole_radius
    for c in range(n_circ):
        ind1 = int(max([0.0, xx[c] - r - 1]))
        ind2 = int(min([info[1] - 1, xx[c] + r + 1]))
        for i in np.arange(ind1, ind2, 1):
            ind3 = int(max([0.0, yy[c] - r - 1]))
            ind4 = int(min([info[2] - 1, yy[c] + r + 1]))
            for j in np.arange(ind3, ind4, 1):
                r_d = np.sqrt(
                    (float(i) - float(xx[c])) ** 2 + (float(j) - float(yy[c])) ** 2
                )
                if r_d <= r:
                    im[i, j] = sz

    if display:
        import matplotlib.pyplot as plt

        plt.figure()
        plt.imshow(im)

    return im

plot_seeing_cond(cond, lim_seeing=None)

Plot calibrator and target seeing conditions.

Parameters:

Name	Type	Description	Default
cond	object	Seeing-condition data.	required
lim_seeing	float or None	Upper seeing limit.	None

Returns:

Type	Description
Figure	Seeing-condition figure.

Source code in src/amical/tools.py

def plot_seeing_cond(cond, lim_seeing=None):  # pragma: no cover
    """Plot calibrator and target seeing conditions.

    Parameters
    ----------
    cond : object
        Seeing-condition data.
    lim_seeing : float or None, default=None
        Upper seeing limit.

    Returns
    -------
    matplotlib.figure.Figure
        Seeing-condition figure.
    """
    import matplotlib.pyplot as plt

    fig = plt.figure()
    ax1 = plt.gca()
    ax1.set_xlabel("mjd [days]")
    ax2 = ax1.twinx()
    ax2.set_ylabel('Seeing ["]', color="#c62d42")
    ax2.tick_params(axis="y", labelcolor="#c62d42")
    ax1.set_ylabel("Uncalibrated mean V$^2$")

    for x in cond:
        ax1.plot(x.mjd, x.vis2, ".", label=x.target)
        ax2.plot(x.mjd, x.seeing, "+", color="#c62d42")
    if lim_seeing is not None:
        ax2.axhline(lim_seeing, color="g", label="Seeing threshold")
    ax1.set_ylim(0, 1.2)
    ax2.set_ylim(0.6, 1.8)
    ax1.grid(alpha=0.1, color="grey")
    ax1.legend(loc="best", fontsize=9)
    plt.tight_layout()
    plt.show(block=False)
    return fig

rad2mas(rad)

Convert radians to milliarcseconds.

Parameters:

Name	Type	Description	Default
rad	float or array_like	Angle in radians.	required

Returns:

Type	Description
float or array_like	Angle in milliarcseconds.

Source code in src/amical/tools.py

def rad2mas(rad):
    """Convert radians to milliarcseconds.

    Parameters
    ----------
    rad : float or array_like
        Angle in radians.

    Returns
    -------
    float or array_like
        Angle in milliarcseconds.
    """
    mas = rad * (3600.0 * 180 / np.pi) * 10.0**3
    return mas

roundSciDigit(number)

Round a number to its significant-digit precision.

Parameters:

Name	Type	Description	Default
number	float	Number to round.	required

Returns:

Name	Type	Description
rounded	float	Rounded number.
sig_digit	int	Significant-digit count.

Source code in src/amical/tools.py

def roundSciDigit(number):
    """Round a number to its significant-digit precision.

    Parameters
    ----------
    number : float
        Number to round.

    Returns
    -------
    rounded : float
        Rounded number.
    sig_digit : int
        Significant-digit count.
    """
    ff = str(number).split(".")[0]
    d = str(number).split(".")[1]
    d, ff = m.modf(number)
    if ff == 0:
        res = str(d).split(".")[1]
        for i in range(len(res)):
            if float(res[i]) != 0.0:
                sig_digit = i + 1
                break
    else:
        sig_digit = 1

    return float(np.round(number, sig_digit)), sig_digit

sanitize_array(dic)

Convert NumPy boolean values in a nested dictionary to Python booleans.

Parameters:

Name	Type	Description	Default
dic	dict	Dictionary to sanitize.	required

Returns:

Type	Description
dict	Sanitized dictionary.

Source code in src/amical/tools.py

def sanitize_array(dic):  # pragma: no cover
    """Convert NumPy boolean values in a nested dictionary to Python booleans.

    Parameters
    ----------
    dic : dict
        Dictionary to sanitize.

    Returns
    -------
    dict
        Sanitized dictionary.
    """
    d2 = dic.copy()
    for k, v in dic.items():
        if isinstance(v, np.ndarray):
            d2[k] = sanitize_array(v)
        if isinstance(v, list):
            d2[k] = np.array(v)
    return d2

save_bs_hdf5(bs, filename)

Save extraction results to an HDF5 file.

Parameters:

Name	Type	Description	Default
bs	object	Extraction result.	required
filename	str or path - like	Output filename.	required

Returns:

Type	Description
None	The result is written to disk.

Source code in src/amical/tools.py

def save_bs_hdf5(bs, filename):
    """Save extraction results to an HDF5 file.

    Parameters
    ----------
    bs : object
        Extraction result.
    filename : str or path-like
        Output filename.

    Returns
    -------
    None
        The result is written to disk.
    """
    import h5py

    if ".h5" not in filename:
        filename += ".h5"

    # Step 1 - Create hdf5 and hierarchy tree/grp
    hf = h5py.File(filename, "w")

    grp_obs = hf.create_group("obs")
    grp_matrix = hf.create_group("matrix")
    grp_info = hf.create_group("infos")
    grp_mask = hf.create_group("mask")
    grp_hdr = hf.create_group("hdr")

    # Step 2 - Save observable (not in tree matrix, infos and mask).
    for key in bs:
        if key not in ["matrix", "mask", "infos"]:
            grp_obs.create_dataset(key, data=bs[key])

    # Step 3 - Save matrix (contains all individual observable
    # frame by frame and statistics matrix (covariance, variance, etc.)).
    matrix = bs.matrix
    for mat in matrix:
        if matrix[mat] is not None:
            grp_matrix.create_dataset(mat, data=matrix[mat])
        else:
            grp_matrix.create_dataset(mat, data=[0])

    # Step 4 - Save mask (contains all mask informations) as
    # well as u1, v1, u2, v2 coordinates (for the CP).
    mask = bs.mask
    for key in mask:
        if (mask[key] is not None) and (key != "t3_coord"):
            grp_mask.create_dataset(key, data=mask[key])

    t3_coord = mask["t3_coord"]
    for key in t3_coord:
        grp_mask.create_dataset(key, data=t3_coord[key])

    # Step 5 - Save informations (target, date, etc.).
    infos = bs.infos
    for key in infos:
        if key != "hdr":
            info = infos[key]
            if info is None:
                info = "Unknown"
            grp_info.attrs[key] = info

    # Step 6 - Save original header keywords.
    hdr = bs.infos.hdr
    for key in hdr:
        grp_hdr.attrs[key] = str(hdr[key])

    # Last step - close hdf5
    hf.close()
    return None

sphere_parang(hdr, n_dit_ifs=None)

Compute parallactic angles for SPHERE frames.

Parameters:

Name	Type	Description	Default
hdr	Header	FITS header.	required
n_dit_ifs	int or None	Number of IFS integrations.	None

Returns:

Type	Description
ndarray	Frame number, frame time, and parallactic angle.

Source code in src/amical/tools.py

def sphere_parang(hdr, n_dit_ifs=None):
    """Compute parallactic angles for SPHERE frames.

    Parameters
    ----------
    hdr : astropy.io.fits.Header
        FITS header.
    n_dit_ifs : int or None, default=None
        Number of IFS integrations.

    Returns
    -------
    numpy.ndarray
        Frame number, frame time, and parallactic angle.
    """
    from astropy.time import Time

    r2d = 180 / np.pi
    d2r = np.pi / 180

    detector = hdr["HIERARCH ESO DET ID"]
    if detector.strip() == "IFS":
        offset = 135.87 - 100.46  # from the SPHERE manual v4
    elif detector.strip() == "IRDIS":
        # correspond to the difference between the PUPIL tracking ant the FIELD tracking for IRDIS taken here: http://wiki.oamp.fr/sphere/AstrometricCalibration (PUPOFFSET)
        offset = 135.87
    else:
        offset = 0
        print(
            "WARNING: Unknown instrument in create_parang_list_sphere: " + str(detector)
        )

    try:
        # Get the correct RA and Dec from the header
        actual_ra = hdr["HIERARCH ESO INS4 DROT2 RA"]
        actual_dec = hdr["HIERARCH ESO INS4 DROT2 DEC"]

        # These values were in weird units: HHMMSS.ssss
        actual_ra_hr = np.floor(actual_ra / 10000.0)
        actual_ra_min = np.floor(actual_ra / 100.0 - actual_ra_hr * 100.0)
        actual_ra_sec = actual_ra - actual_ra_min * 100.0 - actual_ra_hr * 10000.0

        ra_deg = (
            (actual_ra_hr + actual_ra_min / 60.0 + actual_ra_sec / 60.0 / 60.0)
            * 360.0
            / 24.0
        )

        # the sign makes this complicated, so remove it now and add it back at the end
        sgn = np.sign(actual_dec)
        actual_dec *= sgn

        actual_dec_deg = np.floor(actual_dec / 10000.0)
        actual_dec_min = np.floor(actual_dec / 100.0 - actual_dec_deg * 100.0)
        actual_dec_sec = actual_dec - actual_dec_min * 100.0 - actual_dec_deg * 10000.0

        dec_deg = (
            actual_dec_deg + actual_dec_min / 60.0 + actual_dec_sec / 60.0 / 60.0
        ) * sgn
        geolat_rad = float(hdr["ESO TEL GEOLAT"]) * d2r
    except Exception:
        print("WARNING: No RA/Dec Keywords found in header")
        ra_deg = 0
        dec_deg = 0
        geolat_rad = 0

    if "NAXIS3" in hdr:
        if detector.strip() == "IFS":
            n_frames = n_dit_ifs
        else:
            n_frames = hdr["NAXIS3"]
    else:
        n_frames = 1

    # We want the exposure time per frame, derived from the total time from when the shutter
    # opens for the first frame until it closes at the end.
    # This is what ACC thought should be used
    # total_exptime = hdr['ESO DET SEQ1 EXPTIME']
    # This is what the SPHERE DC uses
    total_exptime = (
        Time(hdr["HIERARCH ESO DET FRAM UTC"]) - Time(hdr["HIERARCH ESO DET SEQ UTC"])
    ).sec
    # print total_exptime-total_exptime2
    delta_dit = total_exptime / n_frames
    dit = hdr["ESO DET SEQ1 REALDIT"]

    # Set up the array to hold the parangs
    parang_array = np.zeros(n_frames)

    # Output for debugging
    hour_angles = []

    if ("ESO DET SEQ UTC" in hdr.keys()) and ("ESO TEL GEOLON" in hdr.keys()):
        # The SPHERE DC method
        jd_start = Time(hdr["ESO DET SEQ UTC"]).jd
        lst_start = jd2lst(hdr["ESO TEL GEOLON"], jd_start) * 3600
        # Use the old method
        lst_start = float(hdr["LST"])
    else:
        lst_start = 0.0
        print("WARNING: No LST keyword found in header")

    # delta dit and dit are in seconds so we need to multiply them by this factor to add them to an LST
    time_to_lst = (24.0 * 3600.0) / (86164.1)

    if "ESO INS4 COMB ROT" in hdr.keys() and hdr["ESO INS4 COMB ROT"] == "PUPIL":
        for i in range(n_frames):
            ha_deg = (
                (lst_start + i * delta_dit * time_to_lst + time_to_lst * dit / 2.0)
                * 15.0
                / 3600
            ) - ra_deg
            hour_angles.append(ha_deg)

            # VLT TCS formula
            f1 = float(np.cos(geolat_rad) * np.sin(d2r * ha_deg))
            f2 = float(
                np.sin(geolat_rad) * np.cos(d2r * dec_deg)
                - np.cos(geolat_rad) * np.sin(d2r * dec_deg) * np.cos(d2r * ha_deg)
            )
            pa = -r2d * np.arctan2(-f1, f2)

            pa = pa + offset

            # Also correct for the derotator issues that were fixed on 12 July 2016 (MJD = 57581)
            if hdr["MJD-OBS"] < 57581:
                alt = hdr["ESO TEL ALT"]
                drot_begin = hdr["ESO INS4 DROT2 BEGIN"]
                # Formula from Anne-Lise Maire
                correction = (
                    np.arctan(np.tan((alt - 2 * drot_begin) * np.pi / 180))
                    * 180
                    / np.pi
                )
                pa += correction

            pa = (pa + 360) % 360
            parang_array[i] = pa

    else:
        if "ARCFILE" in hdr.keys():
            print(hdr["ARCFILE"] + " does seem to be taken in pupil tracking.")
        else:
            print("Data does not seem to be taken in pupil tracking.")

        for i in range(n_frames):
            parang_array[i] = 0

    # And a sanity check at the end
    try:
        # The parang start and parang end refer to the start and end of the sequence, not in the middle of the first and last frame.
        # So we need to correct for that
        expected_delta_parang = (
            (hdr["HIERARCH ESO TEL PARANG END"] - hdr["HIERARCH ESO TEL PARANG START"])
            * (n_frames - 1)
            / n_frames
        )
        delta_parang = parang_array[-1] - parang_array[0]
        if np.abs(expected_delta_parang - delta_parang) > 1.0:
            print(
                "WARNING! Calculated parallactic angle change is >1degree more than expected!"
            )

    except Exception:
        pass

    return parang_array

super_gaussian(x, sigma, m=3.0, amp=1.0, x0=0.0)

Function for creating a super-Gaussian window.

Parameters:

Name	Type	Description	Default
x	ndarray	2D array with the distances of each pixel to the image center.	required
sigma	float	Full width at half maximum (FWHM) of the super-Gaussian window. It is therefore not the standard deviation, as the parameter name would suggest.	required
m	float	Exponent used for the super-Gaussian function (default: 3.0).	3.0
amp	float	Amplitude of the Gaussian function (default: 1.0).	1.0
x0	float	Offset applied to the distances (default: 0.0)	0.0

Returns:

Type	Description
ndarray	2D array with the super-Gaussian window function.

Source code in src/amical/tools.py

def super_gaussian(
    x: np.ndarray, sigma: float, m: float = 3.0, amp: float = 1.0, x0: float = 0.0
) -> np.ndarray:
    """
    Function for creating a super-Gaussian window.

    Parameters
    ----------
    x : np.ndarray
        2D array with the distances of each pixel to the image center.
    sigma : float
        Full width at half maximum (FWHM) of the super-Gaussian
        window. It is therefore not the standard deviation, as the
        parameter name would suggest.
    m : float
        Exponent used for the super-Gaussian function (default: 3.0).
    amp : float
        Amplitude of the Gaussian function (default: 1.0).
    x0 : float
        Offset applied to the distances (default: 0.0)

    Returns
    -------
    np.ndarray
        2D array with the super-Gaussian window function.
    """

    return amp * (
        (np.exp(-(2 ** (2 * m - 1)) * np.log(2) * (((x - x0) ** 2) / (sigma**2)) ** m))
        ** 2
    )

wtmn(values, weights)

Compute an inverse-variance weighted mean and standard deviation.

Parameters:

Name	Type	Description	Default
values	ndarray	Values to average.	required
weights	ndarray	1-sigma uncertainties of values, with the same shape. The average is weighted by their inverse squares.	required

Returns:

Type	Description
mean, std : numpy.ndarray	Weighted mean and unbiased standard deviation.

Source code in src/amical/tools.py

def wtmn(values, weights):
    """Compute an inverse-variance weighted mean and standard deviation.

    Parameters
    ----------
    values : numpy.ndarray
        Values to average.
    weights : numpy.ndarray
        1-sigma uncertainties of `values`, with the same shape. The average is
        weighted by their inverse squares.

    Returns
    -------
    mean, std : numpy.ndarray
        Weighted mean and unbiased standard deviation.
    """
    inv_var = 1.0 / np.asarray(weights) ** 2
    mn = np.average(values, weights=inv_var, axis=0)

    ndim = values.ndim

    # Fast and numerically precise:
    variance = np.average((values - mn) ** 2, weights=inv_var, axis=0)
    std = np.sqrt(variance)

    if ndim == 2:
        e_max = []
        for i in range(values.shape[1]):
            e_max.append(np.max(weights[:, i]))
        std_unbias = []
        for i in range(values.shape[1]):
            std_unbias.append(np.max([e_max[i], std[i]]))
        std_unbias = np.array(std_unbias)
    else:
        e_max = weights.copy()
        std_unbias = np.max([np.max(e_max), std])
    return (mn, std_unbias)