def get_indices(w1, w2, header):
    """Returns wavelength layer indices for two given wavelengths in Angstrom.

    Args:
        w1 (float): Lower wavelength, in Angstrom.
        w2 (float): Upper wavelength, in Angstrom.
        header (astropy.io.fits.Header): FITS header for this data cube.

    Returns:
        int tuple: The lower and upper wavelength indices for this range.

    Example:

        To find the indices corresponding to the wavelength range 4100-4200A:

        >>> from astropy.io import fits
        >>> from cwitools import cubes
        >>> mydata,header = fits.getdata("mydata.fits",header=True)
        >>> a,b = cubes.get_band(4100,4200,header)
        >>> mydata_trimmed = mydata[a:b]

    """
    w0, dw, p0 = header['CRVAL3'], header['CD3_3'], header['CRPIX3']
    w0 -= p0 * dw
    return int((w1 - w0) / dw), int((w2 - w0) / dw)
