API References

def extractCutout(cube: SpectralCube,
                  glon: float, glat: float,
                  width: float,
                  v_min: Optional[float] = None, v_max: Optional[float] = None) -> SpectralCube:
    """
    Extract a square cutout from a larger cube.

    Parameters
    ----------
    cube:
        The original SpectralCube to extract from.
    glon:
        Galactic longitude of the center.
    glat:
        Galactic latitude of the center.
    width:
        Width in degrees.
    v_min:
        Minimum velocity in km/s (default = -49.5 km/s).
    v_max:
        Maxiumum velocity in km/s (default = 150 km/s).

    Returns
    -------
    subcube: spectral_cube.SpectralCube
        Returns the extracted subcube.
    """
    glon_min = glon - width / 2
    glon_max = glon + width / 2
    glat_min = glat - width / 2
    glat_max = glat + width / 2

    if v_min is None:
        v_min = -49.5  # [km/s]
    if v_max is None:
        v_max = 150  # [km/s]

    return cube.subcube(xlo=glon_min * u.deg, xhi=glon_max * u.deg,
                        ylo=glat_min * u.deg, yhi=glat_max * u.deg,
                        zlo=v_min * u.km / u.s, zhi=v_max * u.km / u.s)

def getRgalFromVlsr(vlsr, glon, glat=0., R_0=8.34, t_0=240., precision=0.001):
    '''
    Calculate Galactocentric radius from line of sight velocity and Galactic longitude.

    @param vlsr: Line of sight velocity with respect to local standard of rest [km s^-1]
    @type vlsr: float
    @param glon: Galactic longitude [deg]
    @type glon: float
    @param glat: (optional; default glat=0) Galactic latitude [deg]
    @type glat: float
    @param R_0: (optional; default R_0=8.34) Distance to Galactic center [kpc]
    @type R_0: float
    @param t_0: (optional; default Theta_0=248.) Circular rotation velocity at orbit of
                the sun [km s^-1]
    @type t_0: float
    @param precission: (optional; default prec=0.001) Precission of the result [km s^-1]
    @type precission: float
    @return: Galactocentric radius [kpc]
    @rtype: float

    '''
    # make it 2D and column vector (i.e. transpose) if array is given
    if hasattr(vlsr, '__len__'):
        # solution for input arrays
        if type(vlsr) != np.array:
            vlsr = np.array(vlsr)
        vlsr = vlsr[np.newaxis].T
    if hasattr(glon, '__len__'):
        # solution for input arrays
        if type(glon) != np.array:
            glon = np.array(glon)
        glon = glon[np.newaxis].T

    distances = np.arange(0.1, 30. + precision, precision)
    residuals = _residual_brand(distances, vlsr, glon, glat=glat, t_0=t_0, R_0=R_0)
    min_index = np.argmin(np.abs(residuals), axis=len(residuals.shape) - 1)
    R_gal_min = distances[min_index]
    print(distances)

    try:
        # solution for input arrays
        # assume no solution is found if in the extreme values
        is_invalid = (min_index == 0) + (min_index == len(distances) - 1)
        R_gal_min[is_invalid] = np.nan
    except TypeError:
        # solution for single float input
        if min_index in (0, len(distances)):
            R_gal_min = np.nan

    return np.round(R_gal_min, int(np.abs(np.log10(precision))))

def getVlsrFromRgal(R_gal, glon, glat=0., R_0=8.34, t_0=240, precision=0.001):
    '''
    Calculate line of sight velocity from Galactocentric radius.

    @param R_gal: Galactocentric radius [kpc]
    @type R_gal: float
    @param glon: Galactic longitude [deg]
    @type glon: float
    @param glat: (optional; default glat=0) Galactic latitude [deg]
    @type glat: float
    @param R_0: (optional; default R_0=8.34) Distance to Galactic center [kpc]
    @type R_0: float
    @param t_0: (optional; default Theta_0=248.) Circular rotation velocity at orbit of
                the sun [km s^-1]
    @type t_0: float
    @param precission: (optional; default prec=0.001) Precission of the result [km s^-1]
    @type precission: float
    @return: Line of sight velocity with respect to local standard of rest [km s^-1]
    @rtype: float

    '''
    # make it 2D and column vector (i.e. transpose) if array is given
    if hasattr(R_gal, '__len__'):
        # solution for input arrays
        if type(R_gal) != np.array:
            R_gal = np.array(R_gal)
        R_gal = R_gal[np.newaxis].T
    if hasattr(glon, '__len__'):
        # solution for input arrays
        if type(glon) != np.array:
            glon = np.array(glon)
        glon = glon[np.newaxis].T

    vlsrs = np.arange(-200., 200. + precision, precision)
    residuals = _residual_brand(R_gal, vlsrs, glon, glat=glat, t_0=t_0, R_0=R_0)
    min_index = np.argmin(np.abs(residuals), axis=len(residuals.shape) - 1)
    vlsr_min = vlsrs[min_index]

    try:
        # solution for input arrays
        # assume no solution is found if in the extreme values
        is_invalid = (min_index == 0) + (min_index == len(vlsrs) - 1)
        vlsr_min[is_invalid] = np.nan
    except TypeError:
        # solution for single float input
        if vlsr_min in (-200., 200.):
            vlsr_min = np.nan

    return np.round(vlsr_min, int(np.abs(np.log10(precision))))

def getGalacticBar(length:Optional[float]=6., angle:Optional[float]=-50) -> Galactocentric:
    """
    Get the position of the Galactic bar of the Milky Way.

    Parameters
    ----------
    length:float
        Length of the bar in kpc.
    angle:float
        Angle of the bar in degrees.

    Returns
    -------
    bar:astropy.coordinates.Galactocentric
        The position of the bar.
    """
    x = np.linspace(-length/2., length/2., int(length*10)) * np.cos(np.radians(angle)) * u.kpc
    y = np.linspace(-length/2., length/2., int(length*10)) * np.sin(np.radians(angle)) * u.kpc
    return Galactocentric(x=x, y=y, z=np.zeros(len(x)) * u.kpc,
                          galcen_distance=8.3*u.kpc, z_sun=27*u.pc)

def getSpiralArm(name:str, model:str='reid2014',
                 glon_min:Optional[float]=-np.inf, glon_max:Optional[float]=np.inf,
                 r_gal_min:Optional[float]=0., r_gal_max:Optional[float]=np.inf,
                 resolution:Optional[float]=0.0001, R_0:Optional[float]=8.34) -> Galactocentric:
    """
    Get the spiral arm model for a given arm name.

    Parameters
    ----------
    name: str
        Name of the spiral arm.
    model: str
        Name of the spiral arm model to be used.
        Available models: reid2014 (default), ck, vallee2015.
    glon_min: float
        Minimum galactic longitude [deg] of the spiral arm in degrees.
    glon_max: float
        Maximum galactic longitude [deg] of the spiral arm in degrees.
    r_gal_min: float
        Minimum galactocentric distance [kpc] of the spiral arm in kpc.
    r_gal_max: float
        Maximum galactocentric distance [kpc] of the spiral arm in kpc.
    resolution: float
        Resolution of the spiral arm model in degrees.
    R_0: float
        Galactocentric distance of the sun in kpc.

    Returns
    -------
    coordinates: astropy.coordinates.Galactocentric
        Coordinates of the spiral arm in the galactocentric frame.
    """
    spiral_arm_params = _getSpiralParameters(model)

    R_ref = spiral_arm_params[name]['R_ref']
    beta_ref = spiral_arm_params[name]['beta_ref']
    pitch = spiral_arm_params[name]['pitch']
    beta_min = spiral_arm_params[name]['beta_min']
    beta_max = spiral_arm_params[name]['beta_max']
    beta = np.arange(beta_max, beta_min, -1. * resolution)  # [deg]

    R_gal = R_ref * np.exp(-1 * (beta - beta_ref) * np.pi / 180. * np.tan(np.pi / 180. * pitch))  # kpc

    #     x = R_gal * -1. * np.cos(beta * np.pi/180.)
    #     y = R_gal * np.sin(beta * np.pi/180.)
    x = -1 * R_gal * np.sin((90 - beta) * np.pi / 180.)
    y = R_gal * np.cos((90 - beta) * np.pi / 180.)

    arm = SkyCoord(frame='galactocentric', x=x*u.kpc, y=y*u.kpc, z=0.0*u.kpc,
                   galcen_distance=R_0*u.kpc,
                   galcen_coord=SkyCoord(frame='galactic', l=0*u.deg, b=0*u.deg).icrs,
                   z_sun=0.*u.kpc,
                   roll=0.).galactocentric
    arm = arm[arm.galactic.l.degree > glon_min]
    arm = arm[arm.galactic.l.degree < glon_max]
    arm = arm[arm.galactic.distance.kpc > r_gal_min]
    arm = arm[arm.galactic.distance.kpc < r_gal_max]
    return arm

def getSpiralArms(model:Optional[str]='reid2014') -> dict:
    """
    Returns the spiral arms in a dictionary.

    Parameters
    ----------
    model: str
        The model to use

    Returns
    -------
    arms: dict
        A dictionary with the spiral arms.
    """
    sagittarius, scutum, outer, perseus = getSpiralArmsDetail(model=model)

    return {'sagittarius': sagittarius,
            'scutum': scutum,
            'perseus': perseus,
            'outer': outer}

def data_is_same(data1, data2, replace_nans=None):
    """
    Compare two data cubes.

    Parameters
    ----------
    data1: numpy.ndarray
        The first data cube.
    data2: numpy.ndarray
        The second data cube.
    replace_nans: float, optional
        Value to replace NaNs with before comparison. Default is None, which means no replacement is done.

    Returns
    -------
    bool
        True if the data cubes are the same, False otherwise.

    """
    if replace_nans is not None:
        data1_test = copy(data1)
        data1_test[np.isnan(data1_test)] = replace_nans

        data2_test = copy(data2)
        data2_test[np.isnan(data2_test)] = replace_nans

        residual = data1_test - data2_test
    else:
        residual = data1 - data2

    if np.min(residual) == np.max(residual) == 0:
        print("Data is the same")
        return True
    else:
        print("Data is not the same.")
        return False

def fits_cubes_are_same(filename1, filename2):
    """
    Compoare two FITS files.

    Parameters
    ----------
    filename1: str
        The first FITS file.
    filename2: str
        The second FITS file.

    Returns
    -------
    bool
        True if the files are the same, False otherwise.
    """
    cube1 = SpectralCube.read(filename1)
    header1 = cube1.hdu.header
    data1 = cube1.hdu.data

    cube2 = SpectralCube.read(filename2)
    header2 = cube2.hdu.header
    data2 = cube2.hdu.data

    headers_same = headers_are_same(header1, header2, exclude=["HISTORY"])
    data_same = data_is_same(data1, data2)

    if headers_same and data_same:
        print("-> Files represent identical astronomical data.")
        return True
    else:
        print("-> Files are not the same.")
        return False

def headers_are_same(header1, header2, exclude=None):
    """
    Compare two FITS headers.

    Parameters
    ----------
    header1: astropy.io.fits.Header
        The first FITS header.
    header2: astropy.io.fits.Header
        The second FITS header.
    exclude: list of str, optional
        List of header keywords to exclude from the comparison. Default is None, which means no keywords are excluded.

    Returns
    -------
    bool
        True if the headers are the same, False otherwise.

    """
    if exclude is None:
        exclude = []

    error_count = 0
    for k in header1.keys():
        if k not in exclude:
            if k not in header2:
                error_count += 1
            else:
                if header1[k] != header2[k]:
                    error_count += 1
                    print(f"found difference [{k}]")
                    print(f"  old: {header1[k]}")
                    print(f"  new: {header2[k]}")

    if error_count == 0:
        print("Headers are the same.")
        return True
    else:
        print("Headers are not the same (see above).")
        return False