roman.py

"""
=================
Roman Instruments
=================

WARNING: This model has not yet been validated against other PSF
         simulations, and uses several approximations (e.g. for
         mirror polishing errors, which are taken from HST).
"""

import os.path
import poppy
import numpy as np

from scipy.interpolate import griddata, RegularGridInterpolator
from astropy.io import fits
import astropy.units as u
import logging

from . import utils
from . import webbpsf_core
from . import distortion
from .optics import _fix_zgrid_NaNs

_log = logging.getLogger('webbpsf')
import pprint

GRISM_FILTERS = ('GRISM0', 'GRISM1')
PRISM_FILTERS = ('PRISM',)

class WavelengthDependenceInterpolator(object):
    """WavelengthDependenceInterpolator can be configured with
    `n_zernikes` worth of Zernike coefficients at up to `n_wavelengths`
    wavelengths, and will let you `get_aberration_terms` for any
    wavelength in range interpolated linearly between measured/known
    points
    """

    def __init__(self, n_wavelengths=16, n_zernikes=22):
        self._n_wavelengths = n_wavelengths
        self._n_zernikes = n_zernikes
        self._aberration_terms = np.zeros((n_wavelengths, n_zernikes), dtype=np.float64)
        self._wavelengths = []

    def set_aberration_terms(self, wavelength, zernike_array):
        """Supply a reference `wavelength` and a `zernike_array`
        (of length `n_zernikes`) where the aberration is known
        """
        n_wavelengths_set = len(self._wavelengths)
        if wavelength not in self._wavelengths and n_wavelengths_set < self._n_wavelengths:
            self._wavelengths.append(wavelength)
            aberration_row_idx = n_wavelengths_set  # which is now index of last row
        elif wavelength in self._wavelengths:
            aberration_row_idx = self._wavelengths.index(wavelength)
        else:
            # can't add more wavelengths without allocating new _aberration_terms array
            raise ValueError("Already have information at {} wavelengths "
                             "(pass larger n_wavelengths to __init__?)".format(self._n_wavelengths))
        if len(zernike_array) != self._n_zernikes:
            raise ValueError("Expected {} aberration terms (pass different "
                             "n_zernikes to __init__?)".format(self._n_zernikes))
        self._aberration_terms[aberration_row_idx] = zernike_array

    def get_aberration_terms(self, wavelength):
        """Return the Zernike coefficients as interpolated for this
        `wavelength`"""
        # return array of length n_zernikes interpolated for this wavelength
        if wavelength in self._wavelengths:
            # aberration known exactly for this wavelength
            aberration_row_idx = self._wavelengths.index(wavelength)
            return self._aberration_terms[aberration_row_idx]
        else:
            # we have to interpolate @ this wavelength
            aberration_terms = griddata(self._wavelengths, self._aberration_terms, wavelength, method='linear')
            if np.any(np.isnan(aberration_terms)):
                if isinstance(wavelength, u.Quantity):
                    wavelength = wavelength.to(u.m).value
                wavelength_closest = np.clip(wavelength, np.min(self._wavelengths), np.max(self._wavelengths))
                _log.warn("Attempted to get aberrations at wavelength {:.2g} "
                          "outside the range of the reference data; clipping to closest wavelength {:.2g}".format(
                    wavelength, wavelength_closest))

                aberration_terms = griddata(self._wavelengths, self._aberration_terms, wavelength_closest,
                                            method='linear')
            return aberration_terms

class FieldDependentAberration(poppy.ZernikeWFE):
    """FieldDependentAberration incorporates aberrations that
    are interpolated in wavelength, x, and y pixel positions by
    computing the Zernike coefficients for a particular wavelength
    and position.
    """

    """By default, `get_aberration_terms` will zero out Z1, Z2, and Z3
    (piston, tip, and tilt) as they are not meaningful for telescope
    PSF calculations (the former is irrelevant, the latter two would
    be handled by a distortion solution). Change
    `_omit_piston_tip_tilt` to False to include the Z1-3 terms."""
    _omit_piston_tip_tilt = True
    _field_position = None

    def __init__(self, pixel_width, pixel_height,
                 name="Field-dependent Aberration", radius=1.0, oversample=1, interp_order=3):
        self.pixel_width, self.pixel_height = pixel_width, pixel_height
        self.field_position = pixel_width // 2, pixel_height // 2
        self._wavelength_interpolators = {}
        self.pupil_diam = radius * 2.0
        super().__init__(
            name=name,
            verbose=True,
            radius=radius,
            oversample=oversample,
            interp_order=interp_order
        )

    def get_opd(self, wave):
        """Set the Zernike coefficients (for ZernikeWFE.getOPD) based
        on the wavelength of the incoming wavefront and the pixel
        position
        """
        if not isinstance(wave, poppy.Wavefront):
            wavelength = wave
        else:
            wavelength = wave.wavelength
        self.coefficients = wavelength * self.get_aberration_terms(wavelength)
        return super().get_opd(wave)

    @property
    def field_position(self):
        return self._field_position

    @field_position.setter
    def field_position(self, position):
        """Set the x and y pixel position on the detector for which to
        interpolate aberrations"""
        x_pixel, y_pixel = position
        if x_pixel > self.pixel_width or x_pixel < 0:
            raise ValueError("Requested pixel_x position lies outside "
                             "the detector width ({})".format(x_pixel))
        if y_pixel > self.pixel_height or y_pixel < 0:
            raise ValueError("Requested pixel_y position lies outside "
                             "the detector height ({})".format(y_pixel))

        self._field_position = x_pixel, y_pixel

    def add_field_point(self, x_pixel, y_pixel, interpolator):
        """Supply a wavelength-space interpolator for a pixel position
        on the detector"""
        self._wavelength_interpolators[(x_pixel, y_pixel)] = interpolator

    def get_aberration_terms(self, wavelength):
        """Supply the Zernike coefficients for the aberration based on
        the wavelength and pixel position on the detector"""
        if self.field_position in self._wavelength_interpolators:
            # short path: this is a known point
            interpolator = self._wavelength_interpolators[self.field_position]
            coefficients = interpolator.get_aberration_terms(wavelength)
        else:
            # get aberrations at all field points
            field_points, aberration_terms = [], []
            for field_point_coords, point_interpolator in self._wavelength_interpolators.items():
                field_points.append(field_point_coords)
                aberration_terms.append(point_interpolator.get_aberration_terms(wavelength))
            aberration_array = np.asarray(aberration_terms)
            assert len(aberration_array.shape) == 2, "computed aberration array is not 2D " \
                                                     "(inconsistent number of Zernike terms " \
                                                     "at each point?)"
            field_position = tuple(self.field_position)
            coefficients = griddata(
                np.asarray(field_points),
                np.asarray(aberration_terms),
                field_position,
                method='linear'
            )
            if np.any(np.isnan(coefficients)):
                # FIND TWO CLOSEST INPUT GRID POINTS:
                dist = []
                corners = field_points[1:]  # use only the corner points
                for i, ip in enumerate(corners):
                    dist.append(np.sqrt(((ip[0] - field_position[0]) ** 2) + ((ip[1] - field_position[1]) ** 2)))
                min_dist_indx = np.argsort(dist)[:2]  # keep two closest points
                # DEFINE LINE B/W TWO POINTS, FIND ORTHOGONAL LINE AT POINT OF INTEREST,
                # AND FIND INTERSECTION OF THESE TWO LINES.
                x1, y1 = corners[min_dist_indx[0]]
                x2, y2 = corners[min_dist_indx[1]]
                dx = x2 - x1
                dy = y2 - y1
                a = (dy * (field_position[1] - y1) + dx * (field_position[0] - x1)) / (dx * dx + dy * dy)
                closest_interp_point = (x1 + a * dx, y1 + a * dy)
                # INTERPOLATE ABERRATIONS TO CLOSEST INTERPOLATED POINT:
                coefficients = griddata(
                    np.asarray(field_points),
                    np.asarray(aberration_terms),
                    closest_interp_point,
                    method='linear')
                # IF CLOSEST INTERPOLATED POINT IS STILL OUTSIDE THE INPUT GRID,
                # THEN USE NEAREST GRID POINT INSTEAD:
                if np.any(np.isnan(coefficients)):
                    coefficients = aberration_terms[min_dist_indx[0] + 1]
                    _log.warn("Attempted to get aberrations at field point {} which is outside the range "
                              "of the reference data; approximating to nearest input grid point".format(field_position))
                else:
                    _log.warn("Attempted to get aberrations at field point {} which is outside the range "
                              "of the reference data; approximating to nearest interpolated point {}".format(
                        field_position, closest_interp_point))
                assert not np.any(np.isnan(coefficients)), "Could not compute aberration " \
                                                           "at field point {}".format(field_position)
        if self._omit_piston_tip_tilt:
            _log.debug("Omitting piston/tip/tilt")
            coefficients[:3] = 0.0  # omit piston, tip, and tilt Zernikes
        return coefficients


def _load_wfi_detector_aberrations(filename):
    from astropy.io import ascii
    zernike_table = ascii.read(filename, encoding='utf-8-sig')
    detectors = {}

    def build_detector_from_table(number, zernike_table):
        """Build a FieldDependentAberration optic for a detector using
        Zernikes Z1-Z22 at various wavelengths and field points"""
        single_detector_info = zernike_table[zernike_table['sca'] == number]
        field_points = set(single_detector_info['field_point'])
        interpolators = {}
        detector = FieldDependentAberration(
            4096,
            4096,
            radius=RomanInstrument.PUPIL_RADIUS,
            name="Field Dependent Aberration (SCA{:02})".format(number)
        )
        for field_id in field_points:
            field_point_rows = single_detector_info[single_detector_info['field_point'] == field_id]
            local_x, local_y = field_point_rows[0]['local_x'], field_point_rows[0]['local_y']
            interpolator = build_wavelength_dependence(field_point_rows)

            midpoint_pixel = 4096 / 2
            # (local_x in mm / 10 um pixel size) -> * 1e2
            # local_x and _y range from -20.44 to +20.44, so adding to the midpoint pixel
            # makes sense to place (-20.44, -20.44) at (4, 4)
            pixx, pixy = (round(midpoint_pixel - local_x * 1e2),
                          round(midpoint_pixel + local_y * 1e2))

            detector.add_field_point(pixx, pixy, interpolator)
        return detector

    def build_wavelength_dependence(rows):
        """Build an interpolator object that interpolates Z1-Z22 in
        wavelength space"""
        wavelengths = set(rows['wavelength'])
        interpolator = WavelengthDependenceInterpolator(n_wavelengths=len(wavelengths),
                                                        n_zernikes=22)
        for row in rows:
            z = np.zeros(22)
            for idx in range(22):
                z[idx] = row['Z{}'.format(idx + 1)]
            interpolator.set_aberration_terms(row['wavelength'] * 1e-6, z)

        return interpolator

    detector_ids = set(zernike_table['sca'])
    for detid in detector_ids:
        detectors["SCA{:02}".format(detid)] = build_detector_from_table(detid, zernike_table)

    return detectors


@utils.combine_docstrings
class RomanInstrument(webbpsf_core.SpaceTelescopeInstrument):
    PUPIL_RADIUS = 2.4 / 2.0
    """
    RomanInstrument contains data and functionality common to Roman
    instruments, such as setting the pupil shape
    """
    telescope = "Roman"

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.options['jitter'] = 'gaussian'
        self.options['jitter_sigma'] = 0.012 # arcsec/axis, see https://roman.ipac.caltech.edu/sims/Param_db.html#telescope

    def calc_psf(self, outfile=None, source=None, nlambda=None, monochromatic=None,
                 fov_arcsec=None, fov_pixels=None, oversample=None, detector_oversample=None, fft_oversample=None,
                 overwrite=True, display=False, save_intermediates=False, return_intermediates=False,
                 normalize='first', add_distortion=True, crop_psf=False):
        """
        Compute a PSF

        Parameters
        ----------
        add_distortion : bool
            If True, will add 2 new extensions to the PSF HDUlist object. The 2nd extension
            will be a distorted version of the over-sampled PSF and the 3rd extension will
            be a distorted version of the detector-sampled PSF.
        crop_psf : bool
            Included for API compatibility with the JWST instrument classes, but has no
            effect on the results for Roman WFI PSF calculations.

        """

        # Save new keywords to the options dictionary
        self.options['add_distortion'] = add_distortion
        self.options['crop_psf'] = crop_psf

        # add_distortion keyword is not implemented for RomanCoronagraph Class
        if self.name == "RomanCoronagraph" and add_distortion==True:
            self.options['add_distortion'] = False
            self.options['crop_psf'] = False
            _log.info("Geometric distortions are not implemented in WebbPSF for Roman CGI. The add_distortion keyword must be set to False for this case.")
        
        # Run poppy calc_psf
        psf = webbpsf_core.SpaceTelescopeInstrument.calc_psf(self, outfile=outfile, source=source, nlambda=nlambda,
                                                monochromatic=monochromatic, fov_arcsec=fov_arcsec,
                                                fov_pixels=fov_pixels, oversample=oversample,
                                                detector_oversample=detector_oversample, fft_oversample=fft_oversample,
                                                overwrite=overwrite, display=display,
                                                save_intermediates=save_intermediates, 
                                                return_intermediates=return_intermediates, normalize=normalize)


        return psf

    # slightly different versions of the following two functions
    # from the parent superclass
    # in order to interface with the FieldDependentAberration class
    @property
    def detector_position(self):
        """The pixel position in (X, Y) on the detector"""
        return self._detectors[self._detector].field_position

    @detector_position.setter
    def detector_position(self, position):
        # exact copy of superclass function except we save the
        # into a different location.
        try:
            x, y = map(int, position)
        except ValueError:
            raise ValueError("Detector pixel coordinates must be pairs of nonnegative numbers, "
                             "not {}".format(position))
        if x < 0 or y < 0:
            raise ValueError("Detector pixel coordinates must be nonnegative integers")
        if x > self._detector_npixels - 1 or y > self._detector_npixels - 1:
            raise ValueError("The maximum allowed detector pixel "
                             "coordinate value is {}".format(self._detector_npixels - 1))

        self._detectors[self._detector].field_position = (int(position[0]), int(position[1]))

    def _get_aberrations(self):
        """Get the OpticalElement that applies the field-dependent
        optical aberrations. (Called in get_optical_system.)"""
        return self._detectors[self._detector]

    def _get_fits_header(self, result, options):
        """Populate FITS Header keywords"""
        super()._get_fits_header(result, options)
        result[0].header['DETXPIXL'] = (self.detector_position[0],
                                        'X pixel position (for field dependent aberrations)')
        result[0].header['DETYPIXL'] = (self.detector_position[1],
                                        'Y pixel position (for field dependent aberrations)')
        result[0].header['DETECTOR'] = (self.detector, 'Detector selected')


    def _calc_psf_format_output(self, result, options):
        """
        Add distortion to the created 1-extension PSF

        Apply desired formatting to output file:
                 - rebin to detector pixel scale if desired
                 - set up FITS extensions if desired
                 - output either the oversampled, rebinned, or both
        Which image(s) get output depends on the value of the options['output_mode']
        parameter. It may be set to 'Oversampled image' to output just the oversampled image,
        'Detector sampled image' to output just the image binned down onto detector pixels, or
        'Both as FITS extensions' to output the oversampled image as primary HDU and the
        rebinned image as the first image extension. For convenience, the option can be set
        to just 'oversampled', 'detector', or 'both'.

        Modifies the 'result' HDUList object.

        """
       # Pull values from options dictionary
        add_distortion = options.get('add_distortion', True)
        crop_psf = options.get('crop_psf', True)
        # Add distortion if set in calc_psf
        if add_distortion:
            _log.debug("Adding PSF distortion(s)")
            
            # Set up new extensions to add distortion to:
            n_exts = len(result)
            for ext in np.arange(n_exts):
                hdu_new = fits.ImageHDU(result[ext].data, result[ext].header)  # these will be the PSFs that are edited
                result.append(hdu_new)
                ext_new = ext + n_exts
                result[ext_new].header["EXTNAME"] = result[ext].header["EXTNAME"][0:4] + "DIST"  # change extension name
                _log.debug("Appending new extension {} with EXTNAME = {}".format(ext_new, result[ext_new].header["EXTNAME"]))

            _log.debug("WFI: Adding optical distortion")
            psf_distorted = distortion.apply_distortion(result)  # apply siaf distortion model

            # Edit the variable to match if input didn't request distortion
            # (cannot set result = psf_distorted due to return method)
            [result.append(fits.ImageHDU()) for i in np.arange(len(psf_distorted) - len(result))]
            for ext in np.arange(len(psf_distorted)): result[ext] = psf_distorted[ext]

        # Rewrite result variable based on output_mode set:
        webbpsf_core.SpaceTelescopeInstrument._calc_psf_format_output(self, result, options)

class WFIPupilController:
    """
    This is a helper class for the WFI and is used to swap in
    the correct pupil each time the detector is changed.
    The pupil depends on pupil_mask, detector, and filter;
    pupil_mask is set automatically upon the receipt of the
    detector and filter values selected by the user in the WFI class.
    The user should not interact with this class directly, only
    through the API provided through the WFI class.
    """

    def __init__(self):
        self._datapath = None
        self._pupil_basepath = None

        self._pupil = None
        self._pupil_mask = None # new for webbpsf 1.0

        # Flag to en-/disable automatic selection of the appropriate pupil_mask
        self._auto_pupil = True

        # Flag to en-/disable automatic selection of the appropriate pupil file
        self._auto_pupil_mask = True

        # Save formattable pupil file names for each pupil mask
        self.pupil_file_formatters = {
            'SKINNY': 'RST_WIM_Filter_skinny_{0}.fits.gz',
            'WIDE': 'RST_WIM_Filter_F184_{0}.fits.gz',
            'GRISM': 'RST_WSM_Grism_Grism_{0}.fits.gz',
            'PRISM': 'RST_WSM_Prism_Prism_{0}.fits.gz'}

    @property
    def pupil(self):
        """
        The path to the FITS file containing pupil information for the
        detector/filter combination sent from the WFI class. Cannot be
        directly set by the user.
        """
        return self._pupil

    @pupil.setter
    def pupil(self, value):
        raise AttributeError('Pupil cannot be directly specified. '
                             'Use lock_pupil() instead.')

    @property
    def pupil_mask(self):
        """
        The corresponding mask for the filter sent from the WFI class.
        (See WFI.pupil_mask_list for a list of valid filters.) Cannot
        be directly set by the user.
        """
        return self._pupil_mask

    @pupil_mask.setter
    def pupil_mask(self, name):
        raise AttributeError('Pupil mask cannot be directly specified. '
                             'Use lock_pupil_mask() instead.')

    def _get_filter_mask(self, wfi_filter):
        """
        Returns the appropriate mask for a given WFI filter.

        Parameters
        ----------
        wfi_filter : string
            See WFI.filter_list for a list of valid filters.
        """
        wfi_filter = wfi_filter.upper()

        if wfi_filter in GRISM_FILTERS:
            return 'GRISM'
        elif wfi_filter in PRISM_FILTERS:
            return 'PRISM'
        elif wfi_filter in ['F184', 'F213']:
            return 'WIDE'
        else:
            # this method should only be called after WFI.filter was validated,
            # so we assume all inputs are valid and direct those that don't pass
            # preceding cases to skinny
            return 'SKINNY'

    def set_base_path(self, datapath):
        """
        Sets the root directory of the path to WebbPSF's data files.
        This should be set before this class is used.

        Parameters
        ----------
        datapath : string
            Path to WebbPSF-WFI data files
        """
        self._datapath = datapath
        self._pupil_basepath = os.path.join(self._datapath, "pupils")

    def update_pupil(self, filter, detector):
        """
        Selects the specific pupil file corresponding with a detector
        and filter combination sent from the WFI class. Also finds and
        indirectly sets the proper pupil_mask in the process.

        Parameters
        ----------
        filter : string
            See WFI.filter_list for a list of valid filters.

        detector : string
            See WFI.detector_list for a list of valid detectors.
        """
        if not self._auto_pupil:
            _log.info('Automatic pupil selection was locked; '
                      'using user-provided pupil.')
            return

        if self._pupil_basepath is None:
           raise Exception('update_pupil called before setting pupil file path')

        # change detector string to match file format (e.g., "SCA01" -> "SCA_1")
        det_substr = f"{detector[:3]}_{str(int((detector[3:])))}"

        # figure out proper mask based on filter (or use locked mask if enabled)
        pupil_mask = (self._get_filter_mask(filter) if self._auto_pupil_mask
                      else self.pupil_mask)

        path_formatter = self.pupil_file_formatters[pupil_mask]
        pupil = os.path.join(self._pupil_basepath,
                             path_formatter.format(det_substr))

        self._pupil_mask = pupil_mask
        self._pupil = pupil

        _log.info(f"Using {'' if self._auto_pupil_mask else 'locked '}"
                  f"pupil mask '{pupil_mask}' and detector '{detector}'.")

    def lock_pupil(self, pupil_path):
        """
        Prevents the WFIPupilController class from dynamically updating
        the path to the pupil on any changes to the detector or filter
        selected in the WFI class. Instead, the path remains locked on
        whichever `pupil_path` was provided to this method.

        CAUTION: This is non-standard usage of the WFI class and may
        lead to unexpected behavior.

        Parameters
        ----------
        pupil_path : string
            The custom path to your pupil file.
        """
        self._pupil_mask = None
        self._pupil = pupil_path
        self._auto_pupil = False

    def unlock_pupil(self):
        """
        Undoes the effects of lock_pupil() and resets WFIPupilController
        to its default state of updating the pupil whenever a detector
        or filter is changed in the WFI class.
        """
        self._auto_pupil = True

    def lock_pupil_mask(self, pupil_mask):
        """
        Prevents the WFIPupilController class from dynamically updating
        the pupil mask on any changes to the filter selected in the WFI
        class. Instead, the pupil mask remains locked on whichever
        `pupil_mask` was provided to this method.

        CAUTION: This is non-standard usage of the WFI class and may
        lead to unexpected behavior.

        Parameters
        ----------
        filter : string
            See WFI.pupil_mask_list for a list of valid pupil masks.
        """
        if pupil_mask not in self.pupil_file_formatters.keys():
            raise Exception('invalid pupil mask')
        elif not self._auto_pupil:
            raise Exception('Pupil is locked. Unlock pupil before locking pupil mask.')
        else:
            self._pupil_mask = pupil_mask
            self._auto_pupil_mask = False

    def unlock_pupil_mask(self):
        """
        Undoes the effects of lock_pupil_mask() and resets
        WFIPupilController to its default state of updating the pupil
        mask whenever filter is changed in the WFI class.
        """
        self._auto_pupil_mask = True


class WFI(RomanInstrument):
    """
    WFI represents the Roman mission's Wide Field Imager.

    WARNING: This model has not yet been validated against other PSF
             simulations, and uses several approximations (e.g. for
             mirror polishing errors, which are taken from HST).
    """

    def __init__(self):
        # pixel scale is from Roman-AFTA SDT report final version (p. 91)
        # https://roman.ipac.caltech.edu/sims/Param_db.html
        pixelscale = 110e-3 # arcsec/px

        # Initialize the aberrations for super().__init__
        self._aberration_files = {}
        self._is_custom_aberration = False
        self._current_aberration_file = ""

        super().__init__("WFI", pixelscale=pixelscale)

        # Initialize the pupil controller
        self._pupil_controller = WFIPupilController()
        self._pupil_controller.set_base_path(self._datapath)

        self.pupil_mask_list = list(self._pupil_controller.pupil_file_formatters.keys())

        # Define default aberration files for WFI modes
        self._aberration_files = {
            'imaging': os.path.join(self._datapath, 'wim_zernikes_cycle9.csv'),
            'prism': os.path.join(self._datapath,
                                  'wsm_prism_zernikes_cycle9.csv'),
            'grism': os.path.join(self._datapath,
                                  'wsm_grism_zernikes_cycle9.csv'),
            'custom': None}

        # Load and set default detector from aberration file
        self._detector_npixels = 4096
        self._load_detector_aberrations(self._aberration_files[self.mode])
        self.detector = 'SCA01'

        self.opd_list = [os.path.join(self._WebbPSF_basepath,
                                      'upscaled_HST_OPD.fits')]
        self.pupilopd = self.opd_list[-1]

    def _addAdditionalOptics(self, optsys, **kwargs):
        _log.debug("   No optics added for WFI")
        return optsys, False, None

    def _load_detector_aberrations(self, path):
        """
        Helper function that, given a path to a file containing detector
        aberrations, loads the Zernike values and populates the class'
        dictator list with `FieldDependentAberration` detectors. This
        function achieves this by calling the
        `webbpsf.roman._load_wfi_detector_aberrations` function.

        Users should use the `override_aberrations` function to override
        current aberrations.

        Parameters
        ----------
        path : string
            The path to the file containing detector aberrations.
        """
        detectors = _load_wfi_detector_aberrations(path)
        assert len(detectors.keys()) > 0

        self._detectors = detectors
        self._current_aberration_file = path

    def _validate_config(self, **kwargs):
        """
        Validates that the WFI is configured sensibly.

        This mainly consists of selecting the masked or unmasked pupil
        appropriately based on the wavelengths requested.
        """
        assert self.filter is not None, 'filter is None'
        assert self.detector is not None, 'detector is None'
        self._update_pupil()

        assert self.pupil is not None, 'pupil is None'
        super()._validate_config(**kwargs)

    def _get_filter_mode(self, wfi_filter):
        """
        Given a filter name, returns the WFI mode.

        Parameters
        ----------
        wfi_filter : string
            Name of WFI filter. See WFI.filter_list for valid values.

        Returns
        -------
        mode : string
            Returns 'imaging', 'grism', or 'prism' depending on filter.

        Raises
        ------
        ValueError
            ...if the input filter is not found in the WFI filter list.
        """

        wfi_filter = wfi_filter.upper()
        if wfi_filter in GRISM_FILTERS:
            return 'grism'
        elif wfi_filter in PRISM_FILTERS:
            return 'prism'
        elif wfi_filter in self.filter_list:
            return 'imaging'
        else:
            raise ValueError(f"Instrument {self.name} doesn't have a filter "
                             f"called {wfi_filter}.")

    def _update_pupil(self, filter=None, detector=None):
        if detector is None:
            detector = self.detector
        if filter is None:
            filter = self.filter

        if detector is not None and filter is not None:
            self._pupil_controller.update_pupil(filter=filter,detector=detector)

    @RomanInstrument.detector.setter
    def detector(self, value):
        """
        The current WFI detector. See WFI.detector_list for valid values.
        """
        if value.upper() not in self.detector_list:
            raise ValueError("Invalid detector. Valid detector names are: {}".format(', '.join(self.detector_list)))

        self._detector = value.upper()
        if self._detector is not None:
            self._update_pupil(detector=self._detector)

    @RomanInstrument.filter.setter
    def filter(self, value):
        """
        The current WFI filter. See WFI.filter_list for valid values.
        """
        # Update filter
        value = value.upper()

        if value not in self.filter_list:
            raise ValueError(f"Instrument {self.name} doesn't have a "
                             f"filter called {value}.")

        self._filter = value

        # Update aberrations if self._aberration_files has been initiated (not
        # empty) and if they haven't been locked by user
        if self._aberration_files and not self._is_custom_aberration:

            # identify aberration file for new mode
            mode = self._get_filter_mode(self._filter)
            aberration_file = self._aberration_files[mode]

            # if aberrations are not already loaded for the new mode,
            # load and replace detectors using the new mode's aberration file
            if not os.path.samefile(self._current_aberration_file,
                                    aberration_file):
                self._load_detector_aberrations(aberration_file)

        # Update pupil only if detector was previously loaded
        # ( i.e., skip this step when called by super() )
        if self.detector is not None:
            self._update_pupil(filter=self._filter)

    @property
    def pupil(self):
        """
        The path to the FITS file containing pupil information for the
        detector/filter combination sent from the WFI class. Cannot be
        directly set by the user.
        """
        return self._pupil_controller.pupil

    @pupil.setter
    def pupil(self, value):
        # don't allow pupil to be set until the pupil controller is active. (a
        # parent class tries to set it to None in WFI's preceding super() call)
        if hasattr(self, '_pupil_controller'):
            raise AttributeError('Pupil cannot be directly specified. '
                                 'Use lock_pupil() instead.')

    @property
    def pupil_mask(self):
        """
        The corresponding mask for the current filter. Cannot be
        directly set by the user.
        """
        return self._pupil_controller.pupil_mask

    @pupil_mask.setter
    def pupil_mask(self, name):
        raise AttributeError('Pupil mask cannot be directly specified. '
                             'Use lock_pupil_mask() instead.')

    @property
    def mode(self):
        """
        The current WFI mode. Cannot be directly set by the user.
        """
        return self._get_filter_mode(self.filter)

    @mode.setter
    def mode(self, value):
        raise AttributeError("WFI mode cannot be directly specified; "
                             "it is set by changing filters.")

    def lock_aberrations(self, aberration_path):
        """
        This function loads user provided aberrations from a file and
        locks this instrument to only use the provided aberrations (even
        if the filter or mode change).

        To release the lock and load the default aberrations, use
        unlock_aberrations(). To load new user provided
        aberrations, call this function with the new path.

        To load custom aberrations, please provide a csv file
        containing the detector names, field point positions and Zernike
        values. The file should contain the following column
        names/values (comments in parentheses should not be included):
            - sca (Detector number)
            - wavelength (µm)
            - field_point (field point number/ID for SCA and wavelength,
            starts with 1)
            - local_x (mm, local detector coords)
            - local_y (mm, local detector coords)
            - global_x (mm, global instrument coords)
            - global_y (mm, global instrument coords)
            - axis_local_angle_x (XAN)
            - axis_local_angle_y (YAN)
            - wfe_rms_waves (nm)
            - wfe_pv_waves (waves)
            - Z1 (Zernike phase NOLL coefficients)
            - Z2 (Zernike phase NOLL coefficients)
            - Z3 (Zernike phase NOLL coefficients)
            - Z4 (Zernike phase NOLL coefficients)
              .
              .
              .

        Please refer to the default aberration files for examples. If
        you have the WebbPSF data installed and defined, you can get the
        path to that file by running the following:
            >>> from webbpsf import roman
            >>> wfi = roman.WFI()
            >>> print(wfi._aberration_files["imaging"])

        Warning: You should not edit the default files!
        """
        self._load_detector_aberrations(aberration_path)
        self._aberration_files['custom'] = aberration_path
        self._is_custom_aberration = True

    def unlock_aberrations(self):
         """
         Releases the lock on the detector aberration file location
         and loads the default file.
         """
         aberration_path = self._aberration_files[self.mode]
         self._load_detector_aberrations(aberration_path)
         self._aberration_files['custom'] = None
         self._is_custom_aberration = False

    def lock_pupil(self, pupil_path):
        """
        Prevents dynamic updates of the path to the proper pupil file on
        any changes to the selected detector or filter. Instead, the
        path remains locked on whichever `pupil_path` was provided here.

        WARNING: This is non-standard usage of the WFI class and may
        lead to unexpected behavior.

        Parameters
        ----------
        pupil_path : string
            The custom path to your pupil file.
        """
        if os.path.isfile(pupil_path):
            self._pupil_controller.lock_pupil(pupil_path)
        else:
            raise FileNotFoundError(f"{pupil_path} not found.")

        _log.warning("Disabling default pupil selection behavior.")

    def unlock_pupil(self):
        """
        Undoes the effects of lock_pupil() by resetting the class to
        its default state of updating the pupil whenever a detector or
        filter is changed. If necessary, it also sets the proper pupil
        for the current detector/filter combination.
        """
        self._pupil_controller.unlock_pupil()
        self._update_pupil() # reset pupil

    def lock_pupil_mask(self, pupil_mask):
        """
        Prevents dynamic updates of the pupil mask on any change to the
        selected filter. Instead, the pupil mask remains locked on
        whichever `pupil_mask` was provided here.

        WARNING: This is non-standard usage of the WFI class and may
        lead to unexpected behavior.

        Parameters
        ----------
        filter : string
            See WFI.pupil_mask_list for a list of valid pupil masks.
        """
        self._pupil_controller.lock_pupil_mask(pupil_mask)
        self._update_pupil()

    def unlock_pupil_mask(self):
        """
        Undoes the effects of lock_pupil_mask() and resets the class to
        its default state of updating the pupil mask whenever the filter
        is changed.
        """
        self._pupil_controller.unlock_pupil_mask()
        self._update_pupil() # reset pupil mask


class RomanCoronagraph(RomanInstrument):
    """
    Roman Coronagraph Instrument

    Simulates the PSF of the Roman coronagraph.

    Current functionality is limited to the Shaped Pupil Coronagraph (SPC)
    observing modes, and these modes are only simulated with static, unaberrated
    wavefronts, without relay optics and without DM control. The design
    represented here is an approximation to a baseline concept, and will be
    subject to change based on trades studies and technology development.

    Parameters
    ----------
    mode : str
        Roman Coronagraph Instrument observing mode. If not specified, the
        __init__ function will set this to a default mode 'CHARSPC_F660'
    pixelscale : float
        Detector pixelscale. If not specified, the pixelscale will default to
        0.02 arcsec for configurations usint the IMAGER camera and 0.025 arcsec
        for the IFS.
    fov_arcsec : float
        Field of view in arcseconds. If not specified, the field of view will
        default to 3.20 arcsec for the IMAGER camera and 1.76 arcsec for the IFS.

    """

    camera_list = ['IMAGER', 'IFS']
    filter_list = ['F660', 'F721', 'F770', 'F890']
    apodizer_list = ['CHARSPC', 'DISKSPC']
    fpm_list = ['CHARSPC_F660_BOWTIE', 'CHARSPC_F770_BOWTIE', 'CHARSPC_F890_BOWTIE', 'DISKSPC_F721_ANNULUS']
    lyotstop_list = ['LS30D88']

    _mode_table = {  # MODE             CAMERA    FILTER  APODIZER   FPM             LYOT STOP
        'CHARSPC_F660': ('IFS', 'F660', 'CHARSPC', 'CHARSPC_F660_BOWTIE', 'LS30D88'),
        'CHARSPC_F770': ('IFS', 'F770', 'CHARSPC', 'CHARSPC_F770_BOWTIE', 'LS30D88'),
        'CHARSPC_F890': ('IFS', 'F890', 'CHARSPC', 'CHARSPC_F890_BOWTIE', 'LS30D88'),
        'DISKSPC_F721': ('IMAGER', 'F721', 'DISKSPC', 'DISKSPC_F721_ANNULUS', 'LS30D88')}

    def __init__(self, mode=None, pixelscale=None, fov_arcsec=None, apply_static_opd=False):
        super().__init__("RomanCoronagraph", pixelscale=pixelscale)

        self._detector_npixels = 1024
        self._detectors = {camera: 'placeholder' for camera in self.camera_list}

        self.pupil_mask_list = self.lyotstop_list  # alias for use in webbpsf_core
        self.image_mask_list = self.fpm_list  # alias for use in webbpsf_core
        self.pupil = os.path.join(self._WebbPSF_basepath, 'AFTA_CGI_C5_Pupil_onax_256px_flip.fits')
        if apply_static_opd:
            self.pupilopd = os.path.join(self._WebbPSF_basepath, 'CGI', 'OPD', 'CGI_static_OPD.fits')
        else:
            self.pupilopd = None
        self.aberration_optic = None
        self.options = {'force_coron': True}
        # Allow the user to preemptively override the default instrument FoV and pixel scale
        if fov_arcsec is not None:
            self.fov_arcsec = fov_arcsec
            self._override_fov = True
        else:
            self._override_fov = False
        if pixelscale is not None:
            self._pixelscale = pixelscale
            self._override_pixelscale = True
        else:
            self._override_pixelscale = False

        if mode is None:
            self.print_mode_table()
            _log.info("Since the mode was not specified at instantiation, defaulting to CHARSPC_F660")
            self.mode = 'CHARSPC_F660'
        else:
            self.mode = mode

    @property
    def camera(self):
        """Currently selected camera name"""
        return self._camera

    @camera.setter
    def camera(self, value):
        value = value.upper()  # force to uppercase
        if value not in self.camera_list:
            raise ValueError("Instrument {0} doesn't have a camera called {1}.".format(self.name, value))
        self._camera = value
        if value == 'IMAGER':
            if not hasattr(self, 'fov_arcsec') or not self._override_fov:
                self.fov_arcsec = 3.2
            if not hasattr(self, 'pixelscale') or not self._override_pixelscale:
                self.pixelscale = 0.020  # Nyquist at 465 nm
        else:  # default to 'IFS'
            if not hasattr(self, 'fov_arcsec') or not self._override_fov:
                self.fov_arcsec = 2 * 0.82  # 2015 SDT report, Section 3.4.1.1.1:
                                            # IFS has 76 lenslets across the (2 x 0.82) arcsec FoV.
            if not hasattr(self, 'pixelscale') or not self._override_pixelscale:
                self.pixelscale = 0.025  # Nyquist at 600 nm

    # for coronagraph, there is one detector per camera and it should be set automatically.
    @property
    def detector(self):
        return self.camera

    @detector.setter
    def detector(self, value):
        raise RuntimeError("Can't set detector directly for RomanCoronagraph; set camera instead.")

    @property
    def filter(self):
        """Currently selected filter name"""
        return self._filter

    @filter.setter
    def filter(self, value):
        value = value.upper()  # force to uppercase
        if value not in self.filter_list:
            raise ValueError("Instrument {0} doesn't have a filter called {1}.".format(self.name, value))
        self._filter = value

    @property
    def apodizer(self):
        """Currently selected apodizer name"""
        return self._apodizer

    @apodizer.setter
    def apodizer(self, value):
        value = value.upper()  # force to uppercase
        if value not in self.apodizer_list:
            raise ValueError("Instrument {0} doesn't have a apodizer called {1}.".format(self.name, value))
        self._apodizer = value
        if value == 'DISKSPC':
            self._apodizer_fname = \
                os.path.join(self._datapath, "optics/DISKSPC_SP_256pix.fits.gz")
        else:  # for now, default to CHARSPC
            self._apodizer_fname = \
                os.path.join(self._datapath, "optics/CHARSPC_SP_256pix.fits.gz")

    @property
    def fpm(self):
        """Currently selected FPM name"""
        return self._fpm

    @fpm.setter
    def fpm(self, value):
        value = value.upper()  # force to uppercase
        if value not in self.fpm_list:
            raise ValueError("Instrument {0} doesn't have a FPM called {1}.".format(self.name, value))
        self._fpm = value
        if value.startswith('DISKSPC'):
            self._fpmres = 3
            self._owa = 20.
            self._Mfpm = int(np.ceil(self._fpmres * self._owa))
            self._fpm_fname = \
                os.path.join(self._datapath,
                             "optics/DISKSPC_FPM_65WA200_360deg_-_FP1res{0:d}_evensamp_D{1:03d}_{2:s}.fits.gz".format(
                                 self._fpmres, 2 * self._Mfpm, self.filter))
        else:
            self._fpmres = 4
            self._owa = 9.
            self._Mfpm = int(np.ceil(self._fpmres * self._owa))
            self._fpm_fname = \
                os.path.join(self._datapath,
                             "optics/CHARSPC_FPM_25WA90_2x65deg_-_FP1res{0:d}_evensamp_D{1:03d}_{2:s}.fits.gz".format(
                                 self._fpmres, 2 * self._Mfpm, self.filter))

    @property
    def lyotstop(self):
        """Currently selected Lyot stop name"""
        return self._lyotstop

    @lyotstop.setter
    def lyotstop(self, value):
        # preserve case for this one since we're used to that with the lyot mask names
        if value not in self.lyotstop_list:
            raise ValueError("Instrument {0} doesn't have a Lyot mask called {1}.".format(self.name, value))
        self._lyotstop = value
        self._lyotstop_fname = os.path.join(self._datapath, "optics/SPC_LS_30D88_256pix.fits.gz")

    @property
    def mode_list(self):
        """Available Observation Modes"""
        keys = self._mode_table.keys()
        keys = sorted(keys)
        return keys

    # mode works differently since it's a meta-property that affects the other ones:
    @property
    def mode(self):
        """Currently selected mode name"""
        for modename, settings in self._mode_table.items():
            if (self.camera == settings[0].upper() and self.filter == settings[1].upper() and
                    self.apodizer == settings[2].upper() and self.fpm == settings[3].upper() and
                    self.lyotstop == settings[4]):
                return modename
        return 'Custom'

    @mode.setter
    def mode(self, value):
        if value not in self.mode_list:
            raise ValueError("Instrument {0} doesn't have a mode called {1}.".format(self.name, value))
        settings = self._mode_table[value]
        self.camera = settings[0]
        self.filter = settings[1]
        self.apodizer = settings[2]
        self.fpm = settings[3]
        self.lyotstop = settings[4]
        _log.info('Set the following optical configuration:')
        _log.info('camera = {0}, filter = {1}, apodizer = {2}, fpm = {3}, lyotstop = {4}'.format(\
                  self.camera, self.filter, self.apodizer, self.fpm, self.lyotstop))

    def print_mode_table(self):
        """Print the table of observing mode options and their associated optical configuration"""
        _log.info("Printing the table of Roman Coronagraph Instrument observing modes supported by WebbPSF.")
        _log.info("Each is defined by a combo of camera, filter, apodizer, "
                  "focal plane mask (FPM), and Lyot stop settings:")
        _log.info(pprint.pformat(self._mode_table))

    @property
    def detector_position(self):
        """The pixel position in (X, Y) on the detector"""
        return 512, 512

    @detector_position.setter
    def detector_position(self, position):
        raise RuntimeError("Detector position not adjustable for RomanCoronagraph")

    def _validate_config(self, **kwargs):
        super()._validate_config(**kwargs)

    def _addAdditionalOptics(self, optsys, oversample=4):
        """Add coronagraphic or spectrographic optics for RomanCoronagraph."""

        trySAM = False

        if ('pupil_shift_x' in self.options and self.options['pupil_shift_x'] != 0) or \
                ('pupil_shift_y' in self.options and self.options['pupil_shift_y'] != 0):
            shift = (self.options['pupil_shift_x'], self.options['pupil_shift_y'])
        else:
            shift = None

        # Add the shaped pupil apodizer
        optsys.add_pupil(transmission=self._apodizer_fname, name=self.apodizer, shift=None)

        # Add the FPM
        optsys.add_image(transmission=self._fpm_fname, name=self.fpm)

        # Add Lyot stop
        self.pupil_mask = self.lyotstop
        optsys.add_pupil(transmission=self._lyotstop_fname, name=self.lyotstop, shift=shift)

        # Cast as MatrixFTCoronagraph; this configures the detector
        occ_box_size = 1.
        mft_optsys = poppy.MatrixFTCoronagraph(optsys, oversample=oversample, occulter_box=occ_box_size)

        return mft_optsys, trySAM, occ_box_size

    def _get_aberrations(self):
        """Get the OpticalElement that applies the field-dependent
        optical aberrations. (Called in get_optical_system.)"""
        return None

    def _get_fits_header(self, result, options):
        """Populate FITS Header keywords"""
        super()._get_fits_header(result, options)
        pupil_hdr = fits.getheader(self.pupil)
        apodizer_hdr = fits.getheader(self._apodizer_fname)
        fpm_hdr = fits.getheader(self._fpm_fname)
        lyotstop_hdr = fits.getheader(self._lyotstop_fname)

        result[0].header.set('MODE', self.mode, comment='Observing mode')
        result[0].header.set('CAMERA', self.camera, comment='Imager or IFS')
        result[0].header.set('APODIZER', self.apodizer, comment='Apodizer')
        result[0].header.set('APODTRAN', os.path.basename(self._apodizer_fname),
                             comment='Apodizer transmission')
        result[0].header.set('PUPLSCAL', apodizer_hdr['PUPLSCAL'],
                             comment='Apodizer pixel scale in m/pixel')
        result[0].header.set('PUPLDIAM', apodizer_hdr['PUPLDIAM'],
                             comment='Full apodizer array size, incl padding.')
        result[0].header.set('FPM', self.fpm, comment='Focal plane mask')
        result[0].header.set('FPMTRAN', os.path.basename(self._fpm_fname),
                             comment='FPM transmission')
        result[0].header.set('FPMSCAL', fpm_hdr['PIXSCALE'], comment='FPM spatial sampling, arcsec/pix')
        result[0].header.set('LYOTSTOP', self.lyotstop, comment='Lyot stop')
        result[0].header.set('LSTRAN', os.path.basename(self._lyotstop_fname),
                             comment='Lyot stop transmission')
        result[0].header.set('PUPLSCAL', lyotstop_hdr['PUPLSCAL'],
                             comment='Lyot stop pixel scale in m/pixel')
        result[0].header.set('PUPLDIAM', lyotstop_hdr['PUPLDIAM'],
                             comment='Lyot stop array size, incl padding.')