DESI Instrument Throughput

In [1]:

    

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

import os
import os.path
import re


    

In [3]:

    

import astropy.table
import astropy.units as u
from astropy.io import fits
from astropy import wcs


    

In [4]:

    

def read_raytracing(num_sections=14, num_groups=12):
    
    # Initialize the result arrays.
    wavelength = np.zeros(num_sections)
    band_size = np.zeros(num_sections)
    field_angle = np.empty(num_groups)
    spot_centroid = np.empty((num_sections, num_groups, 2))
    rms_spot_radius = np.empty((num_sections, num_groups))
    rms_spot_size = np.empty((num_sections, num_groups, 2))
    
    # Read the file into memory.
    filename = os.path.join(os.environ['DESIMODEL'],
                            'data', 'inputs', 'throughput', 'raytracing.txt')
    with open(filename, 'r') as f:
        lines = f.readlines()
    print('Read {0} lines from {1}.'.format(len(lines), filename))
    
    for section in range(num_sections):
        i = 86 * section
        # First two section headers are "Whole band", "360-375nm".
        # Remaining headers have the form: 360-400 band
        section_header = lines[i].rstrip()
        if section > 0:
            wlen_lo, wlen_hi = float(section_header[:3]), float(section_header[4:7])
            wavelength[section] = 0.5 * (wlen_lo + wlen_hi)
            band_size[section] = wlen_hi - wlen_lo
        for group in range(num_groups):
            j = i + 7 * group + 2
            # Extract the field angle [deg] from a line of the form:
            # Field coordinate  :  0.00000000E+000  4.50000000E-001
            fa = float(lines[j].rstrip().split()[-1])
            if section == 0:
                field_angle[group] = fa
            else:
                assert field_angle[group] == fa
            # Extract the spot centroid [mm] from a line of the form:
            # Image coordinate  : -1.01677290E-014  1.09734671E+002
            spot_centroid[section, group, :] = [float(x) for x in lines[j + 1].split()[-2:]]
            # Extract the RMS spot radius from a line of the form:
            # RMS Spot Radius   :  1.91079098E+001 µm
            rms_spot_radius[section, group] = float(lines[j+2].split()[-2])
            # Extract the RMS spot X, Y sizes from lines of the form:
            # RMS Spot X Size   :  1.06094918E+001 µm
            rms_spot_size[section, group, 0] = float(lines[j+3].split()[-2])
            rms_spot_size[section, group, 1] = float(lines[j+4].split()[-2])
            
    # Apply units.
    wavelength = wavelength * u.nm
    band_size = band_size * u.nm
    field_angle = field_angle * u.deg
    spot_centroid = spot_centroid * u.mm
    rms_spot_radius = rms_spot_radius * u.um
    rms_spot_size = rms_spot_size * u.um

    return wavelength, band_size, field_angle, spot_centroid, rms_spot_radius, rms_spot_size

wavelength, band_size, field_angle, \
spot_centroid, rms_spot_radius, rms_spot_size = read_raytracing()


    

    




Read 1203 lines from /Users/david/Cosmo/DESI/code/desimodel-git/data/inputs/throughput/raytracing.txt.

In [5]:

    

def plot_raytracing(wavelength, band_size, field_angle, spot_centroid, rms_spot_radius,
                    rms_spot_size, save=None):

    nwlen, npos = rms_spot_radius.shape[:2]
    angle = field_angle.to(u.deg).value
    
    # Initialize plot.
    fig, (top, btm) = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
    
    for i in range(1, nwlen):
        lo = (wavelength[i] - 0.5 * band_size[i]).to(u.nm).value
        hi = (wavelength[i] + 0.5 * band_size[i]).to(u.nm).value
        label = '{:.0f}-{:.0f}nm'.format(lo, hi)
        top.plot(angle, rms_spot_radius[i].to(u.um).value / np.sqrt(2), label=label)
        btm.plot(angle, (spot_centroid[i, :, 1] - spot_centroid[0, :, 1])
                 .to(u.um).value, label=label)
    
    top.plot(angle, rms_spot_radius[0].to(u.um).value / np.sqrt(2),
             'bo-', label='Whole band')
    top.legend(bbox_to_anchor=(0.5, 1), loc='upper center',
               ncol=nwlen // 3, fontsize='small')
    top.set_ylabel('RMS spot radius / sqrt(2) [um]')
    top.set_ylim(6., 20.)
    top.grid()
    btm.set_xlabel('Field Angle [deg]')
    btm.set_ylabel('Lateral color shift [um]')
    btm.set_ylim(-20, +20)
    btm.set_xlim(0., 1.6)
    btm.grid()
    plt.tight_layout()
    if save:
        plt.savefig(save)

plot_raytracing(wavelength, band_size, field_angle, spot_centroid,
                rms_spot_radius, rms_spot_size, save='raytracing.png')


    

In [6]:

    

def get_achromatic_blur(version='v13'):
    # Achromatic blur contributions from cells B14:B26 of the throughput sheet.
    if version == 'v13':
        achromatic_blur = np.array([
            4.728, 2.431, 1.049, 3.990, 0.000, 0.000, 0.001,
            0.422, 0.205, 0.400, 1.561, 3.544]) * u.um
    elif version == 'v10':
        # Achromatic blur contributions from cells B43:B72 of the throughput sheet.
        achromatic_blur = np.array([
            7.629, 2.431, 1.049, 0.667, 0.000, 0.000, 0.001,
            0.422, 0.205, 0.400, 1.561, 4.035]) * u.um
    else:
        raise ValueError('Invalid version: {0}.'.format(version))
    return np.sqrt(np.sum(achromatic_blur ** 2))


    

In [7]:

    

print(f"v10: {get_achromatic_blur('v10'):.3f}  v13: {get_achromatic_blur('v13'):.3f}")


    

    




v10: 9.206 um  v13: 7.789 um

In [8]:

    

def get_achromatic_rms_offset():
    # Achromatic offset contributions from cells B45:B74, omitting B52, of the
    # Throughput sheet from 
    achromatic_offsets = np.array([
        0.000, 0.000, 4.000, 0.000, 2.000, 2.000, 3.000, 1.000, 2.000,
        2.000, 1.500, 2.000, 1.500, 1.420, 0.110, 0.090, 0.000, 0.400, 0.090,
        0.062, 1.173, 0.002, 0.006, 0.195, 0.030, 0.012, 0.391, 0.469]) * u.um
    return np.sqrt(np.sum(achromatic_offsets ** 2))


    

In [9]:

    

print(f'sigma1d = {get_achromatic_rms_offset() / np.sqrt(2):.3f}')


    

    




sigma1d = 5.220 um

In [10]:

    

def extrapolate(x_in, y_in, x_lo, x_hi):
    
    assert x_in.shape == y_in.shape
    
    try:
        x_unit = x_in.unit
    except AttributeError:
        x_unit = 1
    try:
        y_unit = y_in.unit
    except AttributeError:
        y_unit = 1
        
    dx = x_in[1:] - x_in[:-1]
    dy = y_in[1:] - y_in[:-1]
    y_lo = y_in[0] + dy[0] / dx[0] * (x_lo - x_in[0])
    y_hi = y_in[-1] + dy[-1] / dx[-1] * (x_hi - x_in[-1])
    
    x_out = np.hstack(
        (x_lo.to(x_unit).value,
         x_in.to(x_unit).value,
         x_hi.to(x_unit).value)) * x_unit
    y_out = np.hstack(
        (y_lo.to(y_unit).value,
         y_in.to(y_unit).value,
         y_hi.to(y_unit).value)) * y_unit

    return x_out, y_out


    

In [11]:

    

def save_blur(wavelength, field_angle, raytracing_blur,
              wlen_lo=355*u.nm, wlen_hi=985*u.nm,
              filename='DESI-0347_blur.ecsv'):
    
    description = """\
RMS spot size as a function of wavelength and field angle derived from \
DESI-0347-v10. Chromatic contributions are based on a text file prepared \
by Mike Sholl that is also plotted in the geometric_blur tab of the \
throughput spreadsheet in DESI-0347-v13. Achromatic contributions are taken \
from cells B14:B26 added in quadrature."""
    table = astropy.table.Table(meta=dict(description=description))
    
    # Add wavelength points at each end.
    wlen_ext, _ = extrapolate(wavelength, raytracing_blur[:, 0], wlen_lo, wlen_hi)
    
    # Convert wavelength values from nm to Angstrom.
    table['wavelength'] = astropy.table.Column(
        wlen_ext, unit='Angstrom', format='{:.1f}',
        description='Wavelength')
    
    for i, r in enumerate(field_angle):
        column_name = 'r={0:.2f}'.format(r)
        # Extrapolate values at each end.
        _, blur_ext = extrapolate(wavelength, raytracing_blur[:, i], wlen_lo, wlen_hi)
        # Add achromatic blur in quadrature.
        blur = np.sqrt(blur_ext ** 2 + get_achromatic_blur() ** 2)
        table[column_name] = astropy.table.Column(
            blur, unit='micron', format='{:.4f}',
            description='RMS spot size at field angle {0}.'.format(column_name))

    table.write(filename, format='ascii.ecsv', overwrite=True)

save_blur(wavelength[1:], field_angle, rms_spot_radius[1:] / np.sqrt(2))


    

In [12]:

    

def save_offset(wavelength, field_angle, raytracing_offset,
                wlen_lo=355*u.nm, wlen_hi=985*u.nm,
                filename='DESI-0347_offset.ecsv'):
    
    description = """\
Radial centroid offset as a function of wavelength and field angle derived \
from DESI-0347-v10. Values are based on a text file prepared \
by Mike Sholl that is also plotted in the geometric_blur tab of the \
throughput spreadsheet in DESI-0347-v13.  Values are derived from ray \
tracing of the ideal optics and do not include the random achromatic \
errors in cells B43:B72 of the throughput sheet."""
    table = astropy.table.Table(meta=dict(description=description))
    
    # Add wavelength points at each end.
    wlen_ext, _ = extrapolate(wavelength, raytracing_offset[:, 0], wlen_lo, wlen_hi)
    
    table['wavelength'] = astropy.table.Column(
        wlen_ext, unit='Angstrom', format='{:.1f}',
        description='Center of wavelength band')
    
    for i, r in enumerate(field_angle):
        column_name = 'r={0:.2f}'.format(r)
        # Extrapolate values at each end.
        _, offset_ext = extrapolate(wavelength, raytracing_offset[:, i], wlen_lo, wlen_hi)
        table[column_name] = astropy.table.Column(
            offset_ext, unit='micron', format='{:.4f}',
            description='Radial centroid offset at field angle {0}.'.format(column_name))

    table.write(filename, format='ascii.ecsv', overwrite=True)

save_offset(wavelength[1:], field_angle, spot_centroid[1:, :, 1] - spot_centroid[0, :, 1])


    

In [13]:

    

import desimodel.focalplane.sim


    

In [14]:

    

def save_static_offsets(seed=1, fov=3.2 * u.deg):

    # Generate random dx, dy offsets on a square grid.
    dx, dy = desimodel.focalplane.sim.generate_random_centroid_offsets(seed=seed)
    npix = len(dx)
    scale = fov.to(u.deg).value / npix
    
    # Save offsets to a FITS file. HDU0 just has a header.
    header = fits.Header()
    header['COMMENT'] = 'Random focal plane centroid offsets in microns.'
    header['SEED'] = seed
    hdu0 = fits.PrimaryHDU(header=header)
    
    # HDU1, 2 contain dx, dy offsets in microns as 2D image data.
    w = wcs.WCS(naxis=2)
    w.wcs.ctype = ['x', 'y']
    w.wcs.crpix = [npix / 2. + 0.5, npix / 2. + 0.5]
    w.wcs.cdelt = [scale, scale]
    w.wcs.crval = [0., 0.]
    header = w.to_header()
    header['BUNIT'] = 'um'
    header['COMMENT'] = '+x component of focal plane offset.'
    hdu1 = fits.ImageHDU(data=dx.to(u.um).value.astype(np.float32),
                         name='XOFFSET', header=header)
    header['COMMENT'] = '+y component of focal plane offset.'
    hdu2 = fits.ImageHDU(data=dy.to(u.um).value.astype(np.float32),
                         name='YOFFSET', header=header)
    
    hdus = fits.HDUList([hdu0, hdu1, hdu2])
    hdus.writeto('DESI-0347_static_offset_{0}.fits'.format(seed), overwrite=True)


    

In [15]:

    

save_static_offsets(seed=1)
save_static_offsets(seed=2)
save_static_offsets(seed=3)


    

In [16]:

    

def parse_quantity(string, dimensions=None):
    _float_pattern = re.compile(
        '\s*([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)\s*')
    # Look for a valid number starting the string.
    found_number = _float_pattern.match(string)
    if not found_number:
        raise ValueError('Unable to parse quantity.')
    value = float(found_number.group(1))
    unit = string[found_number.end():]
    quantity = astropy.units.Quantity(value, unit)
    if dimensions is not None:
        try:
            if not isinstance(dimensions, astropy.units.Unit):
                dimensions = astropy.units.Unit(dimensions)
            quantity = quantity.to(dimensions)
        except (ValueError, astropy.units.UnitConversionError):
            raise ValueError('Quantity "{0}" is not convertible to {1}.'
                             .format(string, dimensions))
    return quantity


    

In [17]:

    

def load_table_2d(filename, x_column_name, y_column_prefix):

    table = astropy.table.Table.read(filename, format='ascii.ecsv')
    nx = len(table)
    x = table[x_column_name].copy()
    
    # Look for columns whose name has the specified prefix.
    y_value, y_index = [], []
    y_unit, data_unit = 1, 1
    for i, colname in enumerate(table.colnames):
        if colname.startswith(y_column_prefix):
            # Parse the column name as a value.
            y = parse_quantity(colname[len(y_column_prefix):])
            if y_unit == 1:
                y_unit = y.unit
            elif y_unit != y.unit:
                raise RuntimeError('Column unit mismatch: {0} != {1}.'
                                   .format(y_unit, y.unit))
            if data_unit == 1:
                data_unit = table[colname].unit
            elif data_unit != table[colname].unit:
                raise RuntimeError('Data unit mismatch: {0} != {1}.'
                                   .format(data_unit, table[colname].unit))
            y_value.append(y.value)
            y_index.append(i)

    # Prepare the array of y values.
    ny = len(y_value)
    y = np.array(y_value) * y_unit
    
    # Extract values for each x,y pair.
    data = np.empty((nx, ny))
    for j, i in enumerate(y_index):
        data[:, j] = table.columns[i][:]
        
    return x, y, data * data_unit


    

In [18]:

    

def plot_static_offsets(seed=1, subsampling=4, save=None):

    hdus = fits.open('DESI-0347_static_offset_{0}.fits'.format(seed))

    # Extract the offset images.
    dx = hdus['XOFFSET'].data
    dy = hdus['YOFFSET'].data
    dr = np.sqrt(dx ** 2 + dy ** 2)
    
    # Reconstruct the linear WCS.
    hdr = hdus[1].header
    w = wcs.WCS(hdr)
    nxy = len(dx)
    pix = np.arange(nxy)
    xlo, ylo, xhi, yhi = w.wcs_pix2world([[0, 0], [nxy - 1, nxy - 1]], 0).flatten()
    assert xlo == ylo and xhi == yhi
    xy = np.linspace(xlo, xhi, nxy)
    xy_pad = hdr['CDELT1']
    unit = u.Unit(hdr['BUNIT'])
    
    # Downsample.
    dx = dx[::subsampling, ::subsampling]
    dy = dy[::subsampling, ::subsampling]
    dr = dr[::subsampling, ::subsampling]
    xy = xy[::subsampling]

    plt.figure(figsize=(12, 10))
    plt.quiver(xy, xy, dx, dy, dr, headwidth=1.5, headlength=2.)
    plt.xlim(xy[0] - xy_pad, xy[-1] + xy_pad)
    plt.ylim(xy[0] - xy_pad, xy[-1] + xy_pad)
    plt.xlabel('x-offset from plate center [deg]')
    plt.ylabel('y-offset from plate center [deg]')
    plt.gca().set_aspect(1)
    plt.colorbar(pad=0.01, shrink=0.95).set_label(
        'Centroid offset [{0}]'.format(unit))
    plt.tight_layout()
    if save:
        plt.savefig(save)


    

In [19]:

    

plot_static_offsets(2, save='static_offset_2.png')


    

In [20]:

    

def plot_validation(wgt_exp=1.0, save=None):
    
    # Load the new blur and offset data.
    wlen, fangle, blur = load_table_2d('DESI-0347_blur.ecsv', 'wavelength', 'r=')
    wlen, fangle, offset = load_table_2d('DESI-0347_offset.ecsv', 'wavelength', 'r=')

    # Get the DESI-0347 achromatic blurs and centroid offsets.
    blur_v13 = get_achromatic_blur(version='v13')
    blur_v10 = get_achromatic_blur(version='v10')
    offset_rms = np.sqrt(get_achromatic_rms_offset() ** 2 + desimodel.focalplane.sim.default_offset ** 2)
    
    # Calculate field-angle-weighted averages, to account for the larger number of targets
    # at larger field angles.
    wgt = fangle.value ** wgt_exp
    wgt_sum = np.sum(wgt)
    weighted_mean = lambda x: np.sum(wgt * x, axis=1) / wgt_sum
    weighted_rms = lambda x: np.sqrt(
        np.sum(wgt * x ** 2, axis=-1) * wgt_sum - np.sum(wgt * x, axis=-1) ** 2) / wgt_sum
    
    avg_blur_v13 = weighted_mean(blur)
    avg_blur_v10 = weighted_mean(np.sqrt(blur ** 2 - blur_v13 ** 2 + blur_v10 ** 2))
    avg_offset = np.sqrt(weighted_rms(offset) ** 2 + offset_rms ** 2)
    
    # Plot the new field-angle-weighted throughput data.
    plt.plot(wlen.to(u.Angstrom).value,
             avg_blur_v13.to(u.um).value, 'b.-', lw=2, alpha=0.25, label='DESIMODEL blur v13')
    plt.plot(wlen.to(u.Angstrom).value,
             avg_blur_v10.to(u.um).value, 'b.--', lw=2, alpha=0.25, label='DESIMODEL blur v10')
    plt.plot(wlen.to(u.Angstrom).value,
             avg_offset.to(u.um).value, 'r.-', lw=2, alpha=0.25, label='DESIMODEL offset')
    
    # Compare with the field-averaged blur and RMS offsets from DESI-347-v12 C27:P27, 
    wlen347 = np.array([360, 375, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 980])
    blur347 = np.array([14.78, 14.51, 13.91, 13.44, 12.86, 12.51,
                        12.62, 13.03, 13.53, 14.00, 14.37, 14.62, 14.73, 14.83])
    rmsoff347 = np.array([13.475, 13.017, 12.412, 11.715, 10.979, 11.050, 11.701,
                          12.384, 12.874, 13.085, 12.989, 12.759, 12.843, 12.928])
    plt.scatter(10 * wlen347, blur347, c='b', label='DESI347-v13')
    plt.scatter(10 * wlen347, rmsoff347, c='r', label='DESI347-v13')

    plt.legend(bbox_to_anchor=(0.5, 1.0), loc='upper center',
               ncol=2, fontsize='small')
    plt.xlabel('Wavelength [Angstrom]')
    plt.ylabel('Blur, RMS offset [micron]')
    plt.xlim(3400, 10000)
    plt.grid()

    plt.tight_layout()
    if save:
        plt.savefig(save)
    
plot_validation(save='validation.png')