From f11364c268f4d409595a11afd2c2465cf5f4f00a Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Thu, 30 Nov 2023 17:07:22 -0500 Subject: [PATCH 01/12] Start infrastructure for improved IFU sims: JWInstrument_with_IFU class framework --- webbpsf/tests/test_nirspec.py | 8 ++ webbpsf/webbpsf_core.py | 157 ++++++++++++++++++++++++++++------ 2 files changed, 139 insertions(+), 26 deletions(-) diff --git a/webbpsf/tests/test_nirspec.py b/webbpsf/tests/test_nirspec.py index 7b432c4d..b07de169 100644 --- a/webbpsf/tests/test_nirspec.py +++ b/webbpsf/tests/test_nirspec.py @@ -52,3 +52,11 @@ def test_calc_datacube_fast(): waves = np.linspace(3e-6, 5e-6, 3) cube = nrs.calc_datacube_fast(waves, fov_pixels=30, oversample=1, compare_methods=True) + + +def test_mode_switch(): + nrs = webbpsf_core.NIRSpec() + nrs.mode = 'IFU' + assert 'IFU' in nrs.aperturename + nrs.mode = 'imaging' + assert 'IFU' not in nrs.aperturename diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py index dae390a9..06b8dff4 100644 --- a/webbpsf/webbpsf_core.py +++ b/webbpsf/webbpsf_core.py @@ -83,7 +83,6 @@ class SpaceTelescopeInstrument(poppy.instrument.Instrument): configuration can be done by editing the :ref:`SpaceTelescopeInstrument.options` dictionary, either by passing options to ``__init__`` or by directly editing the dict afterwards. """ - telescope = 'Generic Space Telescope' options = {} # options dictionary """ A dictionary capable of storing other arbitrary options, for extensibility. The following are all optional, and @@ -156,10 +155,15 @@ def _get_filters(self): for filter_row in filter_table: filter_filename = os.path.join( - self._WebbPSF_basepath, self.name, 'filters', '{}_throughput.fits'.format(filter_row['filter']) + self._WebbPSF_basepath, + self.name, + 'filters', + '{}_throughput.fits'.format(filter_row['filter']) ) filter_info[filter_row['filter']] = Filter( - name=filter_row['filter'], filename=filter_filename, default_nlambda=filter_row['nlambda'] + name=filter_row['filter'], + filename=filter_filename, + default_nlambda=filter_row['nlambda'] ) filter_list.append(filter_row['filter']) return filter_list, filter_info @@ -172,8 +176,7 @@ def __init__(self, name='', pixelscale=0.064): self.name = name self._WebbPSF_basepath, self._data_version = utils.get_webbpsf_data_path( - data_version_min=DATA_VERSION_MIN, return_version=True - ) + data_version_min=DATA_VERSION_MIN, return_version=True) self._datapath = os.path.join(self._WebbPSF_basepath, self.name) self._image_mask = None @@ -390,7 +393,7 @@ def get_optical_system(self, fft_oversample=2, detector_oversample=None, fov_arc _log.info('Creating optical system model:') self._extra_keywords = OrderedDict() # Place to save info we later want to put - # into the FITS header for each PSF. + # into the FITS header for each PSF. if options is None: options = self.options @@ -399,7 +402,8 @@ def get_optical_system(self, fft_oversample=2, detector_oversample=None, fov_arc _log.debug('Oversample: %d %d ' % (fft_oversample, detector_oversample)) optsys = poppy.OpticalSystem( - name='{telescope}+{instrument}'.format(telescope=self.telescope, instrument=self.name), oversample=fft_oversample + name='{telescope}+{instrument}'.format(telescope=self.telescope, instrument=self.name), + oversample=fft_oversample ) # For convenience offsets can be given in cartesian or radial coords if 'source_offset_x' in options or 'source_offset_y' in options: @@ -502,17 +506,20 @@ def get_optical_system(self, fft_oversample=2, detector_oversample=None, fov_arc else: pass - optsys.add_detector( - self.pixelscale, fov_pixels=fov_pixels, oversample=detector_oversample, name=self.name + ' detector' - ) + optsys.add_detector(self.pixelscale, + fov_pixels=fov_pixels, + oversample=detector_oversample, + name=self.name + " detector") # --- invoke semi-analytic coronagraphic propagation if trySAM and not ('no_sam' in self.options and self.options['no_sam']): # if this flag is set, try switching to SemiAnalyticCoronagraph mode. _log.info('Trying to invoke switch to Semi-Analytic Coronagraphy algorithm') try: - SAM_optsys = poppy.SemiAnalyticCoronagraph(optsys, oversample=fft_oversample, occulter_box=SAM_box_size) - _log.info('SAC OK') + SAM_optsys = poppy.SemiAnalyticCoronagraph(optsys, + oversample=fft_oversample, + occulter_box=SAM_box_size) + _log.info("SAC OK") return SAM_optsys except ValueError as err: _log.warning( @@ -725,15 +732,10 @@ def psf_grid( if single_psf_centered is True: if isinstance(self._detector_npixels, tuple): - ( - det_npix_y, - det_npix_x, - ) = ( - self._detector_npixels - ) # A tuple has been provided for a non-square detector with different Y and X dimensions + det_npix_y, det_npix_x = self._detector_npixels # A tuple has been provided for a non-square detector with different Y and X dimensions else: det_npix_y = det_npix_x = self._detector_npixels # same dimensions in both X and Y - psf_location = (int(det_npix_x - 1) // 2, int(det_npix_y - 1) // 2) # center pt + psf_location = ( int(det_npix_x - 1) // 2, int(det_npix_y - 1) // 2) # center pt else: psf_location = self.detector_position[::-1] # (y,x) @@ -1885,9 +1887,33 @@ def calc_datacube_fast(self, wavelengths, compare_methods=False, *args, **kwargs return cubefast +class JWInstrument_with_IFU(JWInstrument): + """ Subclass which adds some additional infrastructure for IFU sims""" + + + def __init__(self, *args, **kwargs): + super().__init__( *args, **kwargs) + # dict of modes and default aperture names + self._modes_list = {'imaging': None, + 'IFU': None} + self._mode = 'imaging' + + @property + def mode(self): + """Currently selected instrument major mode, imaging or IFU""" + return self._mode + -class MIRI(JWInstrument): - """A class modeling the optics of MIRI, the Mid-InfraRed Instrument. + @mode.setter + def mode(self, value): + if value not in self._modes_list: + raise ValueError(f"'{value} is not an allowed mode for this instrument.") + self._mode = value + self.set_position_from_aperture_name(self._modes_list[value]) + + +class MIRI(JWInstrument_with_IFU): + """ A class modeling the optics of MIRI, the Mid-InfraRed Instrument. Relevant attributes include `filter`, `image_mask`, and `pupil_mask`. @@ -1907,7 +1933,7 @@ class MIRI(JWInstrument): def __init__(self): self.auto_pupil = True - JWInstrument.__init__(self, 'MIRI') + JWInstrument_with_IFU.__init__(self, 'MIRI') self.pixelscale = self._get_pixelscale_from_apername('MIRIM_FULL') self._rotation = 4.83544897 # V3IdlYAngle, Source: SIAF PRDOPSSOC-059 # This is rotation counterclockwise; when summed with V3PA it will yield the Y axis PA on sky @@ -2807,7 +2833,7 @@ def _get_fits_header(self, hdulist, options): hdulist[0].header['PILIN'] = ('False', 'Pupil imaging lens in optical path: T/F') -class NIRSpec(JWInstrument): +class NIRSpec(JWInstrument_with_IFU): """A class modeling the optics of NIRSpec, in **imaging** mode. This is not a substitute for a spectrograph model, but rather a way of simulating a PSF as it @@ -2825,12 +2851,16 @@ class NIRSpec(JWInstrument): """ def __init__(self): - JWInstrument.__init__(self, 'NIRSpec') + JWInstrument_with_IFU.__init__(self, 'NIRSpec') self.pixelscale = 0.10435 # Average over both detectors. SIAF PRDOPSSOC-059, 2022 Dec # Microshutters are 0.2x0.46 but we ignore that here. self._rotation = 138.5 # Average for both detectors in SIAF PRDOPSSOC-059 - # This is rotation counterclockwise; when summed with V3PA it will yield the Y axis PA on sky - self.filter_list.append('IFU') + # This is rotation counterclockwise; when summed with V3PA it will yield the Y axis PA on sky + + # Modes and default SIAF apertures for each + self._modes_list = {'imaging': 'NRS1_FULL', + 'IFU': 'NRS_FULL_IFU'} + self._IFU_pixelscale = 0.1043 # same. self.monochromatic = 3.0 self.filter = 'F110W' # or is this called F115W to match NIRCam?? @@ -2923,6 +2953,81 @@ def _get_fits_header(self, hdulist, options): hdulist[0].header['APERTURE'] = (str(self.image_mask), 'NIRSpec slit aperture name') + @JWInstrument.aperturename.setter + def aperturename(self, value): + """Set SIAF aperture name to new value, with validation. + + This also updates the pixelscale to the local value for that aperture, for a small precision enhancement. + + Similar to superclass function, but handles the more complex situation with NIRSpec apertures and detectors + """ + # Explicitly update detector reference coordinates to the default for the new selected aperture, + # otherwise old coordinates can persist under certain circumstances + + try: + ap = self.siaf[value] + except KeyError: + raise ValueError(f'Aperture name {value} not a valid SIAF aperture name for {self.name}') + + # NIRSpec apertures can either be per detector (i.e. "NRS1_FULL") or for the focal plane but not per detector (i.e. "NRS_FULL_IFU") + + if value[0:4] in ['NRS1', 'NRS2']: + # this is a regular per-detector aperture, so just call the regular code in the superclass + JWInstrument.aperturename.fset(self, value) + else: + # apertures that start with NRS define V2,V3 position, but not pixel coordinates and pixelscale. So we + # still have to use a full-detector aperturename for that. + detector_apername = self.detector + "_FULL" + + # Only update if new value is different + if self._aperturename != value: + # First, check some info from current settings, which we will use below as part of auto pixelscale code + # The point is to check if the pixel scale is set to a custom or default value, + # and if it's custom then don't override that. + # Note, check self._aperturename first to account for the edge case when this is called from __init__ before _aperturename is set + has_custom_pixelscale = self._aperturename and (self.pixelscale != self._get_pixelscale_from_apername(detector_apername)) + + # Now apply changes: + self._aperturename = value + # Update detector reference coordinates + # self.detector_position = (ap.XSciRef, ap.YSciRef) + + # Update DetectorGeometry class + self._detector_geom_info = DetectorGeometry(self.siaf, self._aperturename) + _log.info(f"{self.name} SIAF aperture name updated to {self._aperturename}") + + if not has_custom_pixelscale: + self.pixelscale = self._get_pixelscale_from_apername(detector_apername) + _log.debug(f"Pixelscale updated to {self.pixelscale} based on average X+Y SciScale at SIAF aperture {self._aperturename}") + + + if 'IFU' in self.aperturename: + self._mode = 'IFU' + else: + self._mode = 'imaging' # More to implement here later! + + + + def _tel_coords(self): + """ Convert from science frame coordinates to telescope frame coordinates using + SIAF transformations. Returns (V2, V3) tuple, in arcminutes. + + Note that the astropy.units framework is used to return the result as a + dimensional Quantity. + + Some extra steps for NIRSpec to handle the more complicated/flexible mapping between detector and sky coordinates + """ + + if self.aperturename.startswith("NRS_"): + # These apertures don't map directly to particular detector position in the usual way + # Return coords for center of the aperture reference location + return np.asarray((self._detector_geom_info.aperture.V2Ref, + self._detector_geom_info.aperture.V3Ref)) / 60 * units.arcmin + else: + return super()._tel_coords() + + + class NIRISS(JWInstrument): """A class modeling the optics of the Near-IR Imager and Slit Spectrograph (formerly TFI) From 4415329afcac182688ea3e14a56536003f30989b Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Tue, 20 Feb 2024 12:21:30 -0500 Subject: [PATCH 02/12] implement IFU mode for MIRI, and necessary updates to aperture handling methods --- webbpsf/tests/test_miri.py | 7 +++++ webbpsf/webbpsf_core.py | 61 ++++++++++++++++++++++++++++++++------ 2 files changed, 59 insertions(+), 9 deletions(-) diff --git a/webbpsf/tests/test_miri.py b/webbpsf/tests/test_miri.py index 0cb78997..0788d8a1 100644 --- a/webbpsf/tests/test_miri.py +++ b/webbpsf/tests/test_miri.py @@ -120,3 +120,10 @@ def test_miri_nonsquare_detector(): miri = webbpsf_core.MIRI() miri.detector_position = (1023, 1031) # recall this is X, Y order assert miri.detector_position == (1023, 1031) + +def test_mode_switch(): + miri = webbpsf_core.MIRI() + miri.mode = 'IFU' + assert 'IFU' in miri.aperturename + miri.mode = 'imaging' + assert 'IFU' not in miri.aperturename diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py index 06b8dff4..aacd43ac 100644 --- a/webbpsf/webbpsf_core.py +++ b/webbpsf/webbpsf_core.py @@ -961,9 +961,10 @@ def aperturename(self, value): # Only update if new value is different if self._aperturename != value: if ap.AperType == 'SLIT': - # Special case for SLIT apertures (NIRSpec and MIRI) + # Special case for SLIT apertures (NIRSpec or MIRI). Note, this includes all IFU apertures for NIRSpec # apertures of type SLIT define V2,V3 position, but not pixel coordinates and pixelscale. So we # still have to use a full-detector aperturename for that subset of apertures + # This code path also supports MIRI LRS detector_apername = self.detector + '_FULL' _log.info(f'Aperture {value} is of type SLIT; using {detector_apername} for detector geometry.') @@ -984,6 +985,15 @@ def aperturename(self, value): _log.debug( f'Pixelscale updated to {self.pixelscale} based on average X+Y SciScale at SIAF aperture {detector_apername}' ) + elif ap.AperType == 'COMPOUND' and self.name=='MIRI': + # For MIRI, many of the relevant IFU apertures are of COMPOUND type. + has_custom_pixelscale = False # custom scales not supported for MIRI IFU (yet?) + # Unlike NIRSpec, there simply do not exist full-detector SIAF apertures for the MIRI IFU detectors + _log.info(f'Aperture {value} is of type COMPOUND for MIRI; There do not exist corresponding SIAF apertures, so we ignore setting detector geometry.') + + # Now apply changes: + self._aperturename = value + else: if self.detector not in value: raise ValueError( @@ -1057,23 +1067,35 @@ def set_position_from_aperture_name(self, aperture_name): # setting the detector must happen -before- we set the position detname = aperture_name.split('_')[0] - if detname != 'NRS': # Many NIRSpec slit apertures are defined generally, not for a specific detector + if detname.startswith('MIRIFU'): + if 'CHANNEL1' in aperture_name or 'CHANNEL2' in aperture_name: + self.detector = 'MIRIFUSHORT' + else: + self.detector = 'MIRIFULONG' + + elif detname != 'NRS': # Many NIRSpec slit apertures are defined generally, not for a specific detector self.detector = detname # As a side effect this auto reloads SIAF info, see detector.setter self.aperturename = aperture_name - if self.name != 'NIRSpec' and ap.AperType != 'SLIT': - # Regular imaging apertures, so we can just look up the reference coords directly - self.detector_position = (ap.XSciRef, ap.YSciRef) # set this AFTER the SIAF reload - else: + if (self.name == 'NIRSpec' or self.name=='MIRI') and ap.AperType == 'SLIT': # NIRSpec slit apertures need some separate handling, since they don't map directly to detector pixels # In this case the detector position is not uniquely defined, but we ensure to get reasonable values by # using one of the full-detector NIRspec apertures - _log.debug('Inferring detector position using V coords for SLIT aperture: {ap.V2Ref, ap.V3Ref}') + # This code path also supports MIRI LRS SLIT + _log.debug(f'Inferring detector position using V coords for SLIT aperture: {ap.V2Ref, ap.V3Ref}') ref_in_tel = ap.V2Ref, ap.V3Ref nrs_full_aperture = self.siaf[self.detector + '_FULL'] ref_in_sci = nrs_full_aperture.tel_to_sci(*ref_in_tel) self.detector_position = ref_in_sci + elif self.name == 'MIRI' and ap.AperType == 'COMPOUND': + # MIRI IFU compound apertures need separate handling, since they don't map directoy to detector pixels + # In this case the detector position is not uniquely defined, and there do not exist + # in SIAF any full-detector MIRIFU apertures, so just set values to (512,512) as a placeholder. + self.detector_position = [512, 512] + else: + # Regular imaging apertures, so we can just look up the reference coords directly + self.detector_position = (ap.XSciRef, ap.YSciRef) # set this AFTER the SIAF reload _log.debug('From {} set det. pos. to {} {}'.format(aperture_name, detname, self.detector_position)) @@ -1938,6 +1960,10 @@ def __init__(self): self._rotation = 4.83544897 # V3IdlYAngle, Source: SIAF PRDOPSSOC-059 # This is rotation counterclockwise; when summed with V3PA it will yield the Y axis PA on sky + # Modes and default SIAF apertures for each + self._modes_list = {'imaging': 'MIRIM_FULL', + 'IFU': 'MIRIFU_CHANNEL1A'} + # Coordinate system note: The pupil shifts get applied at the instrument pupil, which is an image of the OTE exit pupil # and is thus flipped in Y relative to the V frame entrance pupil. Therefore flip sign of pupil_shift_y self.options['pupil_shift_x'] = -0.0068 # In flight measurement. See Wright, Sabatke, Telfer 2022, Proc SPIE @@ -1964,8 +1990,10 @@ def __init__(self): # The above tuples give the pixel resolution (perpendicular to the slice, along the slice). # The pixels are not square. - self._detectors = {'MIRIM': 'MIRIM_FULL'} # Mapping from user-facing detector names to SIAF entries. - self.detector = self.detector_list[0] + self._detectors = {'MIRIM': 'MIRIM_FULL', # Mapping from user-facing detector names to SIAF entries. + 'MIRIFUSHORT': 'MIRIFU_CHANNEL1A', # only applicable in IFU mode + 'MIRIFULONG': 'MIRIFU_CHANNEL3A'} # ditto + self.detector = 'MIRIM' self._detector_npixels = (1032, 1024) # MIRI detector is not square self.detector_position = (512, 512) @@ -2222,6 +2250,14 @@ def _get_fits_header(self, hdulist, options): if self.image_mask is not None: hdulist[0].header['TACQNAME'] = ('None', 'Target acquisition file name') + def _get_pixelscale_from_apername(self, apername): + """Simple utility function to look up pixelscale from apername""" + + if 'MIRIFU' in apername: + print('special case for MIRI IFU pixelscales') + return 0.1 + else: + return super()._get_pixelscale_from_apername(apername) class NIRCam(JWInstrument): """A class modeling the optics of NIRCam. @@ -3026,6 +3062,13 @@ def _tel_coords(self): else: return super()._tel_coords() + def _get_pixelscale_from_apername(self, apername): + """Simple utility function to look up pixelscale from apername""" + if 'IFU' in apername: + print('DEBUG - special case for NIRSpec IFU pixelscales') + return super()._get_pixelscale_from_apername('NRS1_FULL') + else: + return super()._get_pixelscale_from_apername(apername) class NIRISS(JWInstrument): From f3d38c9960b3b88152dd7bc7865077345dd049ab Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Wed, 21 Feb 2024 00:23:19 -0500 Subject: [PATCH 03/12] implement MIRI MRS band attribute, and mode switching including aperturename functionality (and tests) --- webbpsf/tests/test_miri.py | 32 +++++++++++ webbpsf/webbpsf_core.py | 113 ++++++++++++++++++++++++++++++++++++- 2 files changed, 144 insertions(+), 1 deletion(-) diff --git a/webbpsf/tests/test_miri.py b/webbpsf/tests/test_miri.py index 0788d8a1..18080375 100644 --- a/webbpsf/tests/test_miri.py +++ b/webbpsf/tests/test_miri.py @@ -123,7 +123,39 @@ def test_miri_nonsquare_detector(): def test_mode_switch(): miri = webbpsf_core.MIRI() + imager_rotation = miri._rotation + + # Explicitly switch mode to IFU miri.mode = 'IFU' assert 'IFU' in miri.aperturename + assert miri.detector =='MIRIFUSHORT' + assert miri.aperturename.startswith('MIRIFU_CH') + assert miri._rotation != imager_rotation + # Explicitly switch back to imaging miri.mode = 'imaging' assert 'IFU' not in miri.aperturename + assert miri.detector =='MIRIM' + assert miri.aperturename.startswith('MIRIM_') + assert miri._rotation == imager_rotation + + # Implicitly switch to IFU + miri.set_position_from_aperture_name('MIRIFU_CHANNEL3B') + assert 'IFU' in miri.aperturename + assert miri.detector =='MIRIFULONG' + assert miri.aperturename == 'MIRIFU_CHANNEL3B' + assert miri._rotation != imager_rotation + + # implicitly switch to imaging + # LRS is an odd case, SLIT aper type but operates like in imaging mode + miri.set_position_from_aperture_name('MIRIM_SLIT') + assert 'IFU' not in miri.aperturename + assert miri.detector =='MIRIM' + assert miri.aperturename.startswith('MIRIM_') + assert miri._rotation == imager_rotation + + # And back to IFU again: + miri.mode = 'IFU' + assert 'IFU' in miri.aperturename + assert miri.detector =='MIRIFUSHORT' + assert miri.aperturename.startswith('MIRIFU_CH') + assert miri._rotation != imager_rotation diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py index aacd43ac..e030b542 100644 --- a/webbpsf/webbpsf_core.py +++ b/webbpsf/webbpsf_core.py @@ -993,6 +993,8 @@ def aperturename(self, value): # Now apply changes: self._aperturename = value + # Update DetectorGeometry class + self._detector_geom_info = DetectorGeometry(self.siaf, self._aperturename) else: if self.detector not in value: @@ -1068,10 +1070,14 @@ def set_position_from_aperture_name(self, aperture_name): detname = aperture_name.split('_')[0] if detname.startswith('MIRIFU'): + self._mode = 'IFU' if 'CHANNEL1' in aperture_name or 'CHANNEL2' in aperture_name: self.detector = 'MIRIFUSHORT' else: self.detector = 'MIRIFULONG' + elif detname.startswith('MIRIM'): + self._mode = 'imaging' + self.detector = detname # As a side effect this auto reloads SIAF info, see detector.setter elif detname != 'NRS': # Many NIRSpec slit apertures are defined generally, not for a specific detector self.detector = detname # As a side effect this auto reloads SIAF info, see detector.setter @@ -1990,6 +1996,25 @@ def __init__(self): # The above tuples give the pixel resolution (perpendicular to the slice, along the slice). # The pixels are not square. + # Mappings between alternate names used for MRS subbands + self._MRS_dichroic_to_subband = {"SHORT": "A", "MEDIUM": "B", "LONG": "C"} + self._MRS_subband_to_dichroic = {"A": "SHORT", "B": "MEDIUM", "C": "LONG"} + + self._band = None + self._MRS_bands = {"1A": [4.887326748103221, 5.753418963216559], # Values provided by Polychronis Patapis + "1B": [5.644625711181792, 6.644794583147869], # To-do: obtain from CRDS pipeline refs + "1C": [6.513777066360325, 7.669147994055998], + "2A": [7.494966046398437, 8.782517027772244], + "2B": [8.651469658142522, 10.168811217793243], + "2C": [9.995281242621394, 11.73039280033565], + "3A": [11.529088518317131, 13.491500288051483], + "3B": [13.272122736770127, 15.550153182343314], + "3C": [15.389530615108631, 18.04357852656418], + "4A": [17.686540162850203, 20.973301482912323], + "4B": [20.671069749545193, 24.476094964546686], + "4C": [24.19608171436692, 28.64871057821349]} + + self._detectors = {'MIRIM': 'MIRIM_FULL', # Mapping from user-facing detector names to SIAF entries. 'MIRIFUSHORT': 'MIRIFU_CHANNEL1A', # only applicable in IFU mode 'MIRIFULONG': 'MIRIFU_CHANNEL3A'} # ditto @@ -2228,8 +2253,12 @@ def _update_aperturename(self): apname = 'MIRIM_FULL' # LRS slit uses full array readout else: apname = self._image_mask_apertures[self._image_mask] - else: + elif self.mode == 'imaging': apname = 'MIRIM_FULL' + else: + # IFU mode is complex, don't try to set a different apname here + # (This gracefully works around an edge case in mode switching) + apname = self.aperturename # Call aperturename.setter to update ap ref coords and DetectorGeometry class self.aperturename = apname @@ -2259,6 +2288,88 @@ def _get_pixelscale_from_apername(self, apername): else: return super()._get_pixelscale_from_apername(apername) + def _get_aperture_rotation(self, apername): + """Get the rotation angle of a given aperture, using values from SIAF. + + Returns ~ position angle counterclockwise from the V3 axis, in degrees + (i.e. consistent with SIAF V3IdlYangle) + + Note, MIRIFU aperture geometry is extremely complex, and this oversimplifies. + See https://jwst-docs.stsci.edu/jwst-mid-infrared-instrument/miri-instrumentation/miri-mrs-field-and-coordinates + Consistent with that reference, we compute the angle of the along-slice (alpha) direction relative to the + horizontal (V2). Because the apertuers are skewed, this yields different values by several degrees + than an equivalent calculation for beta. + """ + if apername.startswith('MIRIM'): + return self.siaf['MIRIM_FULL'].V3IdlYAngle + elif apername.startswith('MIRIFU'): + # These SLIT or COMPOUND apertures do not have a V3IdlYangle param defined in SIAF + # But we can work out the angles from the aperture corner vertices which are defined. + cx, cy = self.siaf[apername].corners('tel', rederive=False) + # The aperture shapes are irregular quadrilaterals, not squares or rectangles + # So, take the angles of both alpha-axis (~horizontal) sides, and average them to get an average rotation angle + dx = cx[0] - cx[1] + dy = cy[0] - cy[1] + dx2 = cx[3] - cx[2] + dy2 = cy[3] - cy[2] + # take the average, and convert to degrees. + # The signs and parity of the arctan2 are atypical here, to match the expected output convention of + # rotation counterclockwise relative to the V2V3 frame. + avg_V3IdlYangle = np.rad2deg((np.arctan2(dy, -dx) + np.arctan2(dy2, -dx2)) / 2) + return avg_V3IdlYangle + else: + raise ValueError(f"Unexpected/invalid apername for MIRI: {apername}") + + + @JWInstrument_with_IFU.aperturename.setter + def aperturename(self, value): + """Set aperturename, also update the rotation for MIRIM vs. IFU channel """ + # apply the provided aperture name + # Note, the syntax for calling a parent class property setter is... painful: + super(MIRI, type(self)).aperturename.fset(self, value) + # Update the rotation angle + self._rotation = self._get_aperture_rotation(self.aperturename) + + # if it's an IFU aperture, we're now in IFU mode: + self._mode = 'IFU' if value.startswith('MIRIFU') else 'imaging' + # if in IFU mode, we probably also want to update the IFU band + + if self.mode == 'IFU': + if self.band != value[-2:]: + self.band = value[-2:] + + + @property + def band(self): + """MRS IFU spectral band. E.g. '1A', '3B'. Only applicable in IFU mode.""" + if self.mode == 'IFU': + return self._band + else: + return None + + @band.setter + def band(self, value): + if self.mode != 'IFU': + if value is not None: + raise RuntimeError("The 'band' property is only valid for IFU mode simulations.") + return + + if value in self._MRS_bands.keys(): + self._band = value + #self._slice_width = self._IFU_pixelscale[f"Ch{self._band[0]}"][0] + + self.aperturename = 'MIRIFU_CHANNEL' + value + # setting aperturename will also auto update self._rotation + #self._rotation = self.MRS_rotation[self._band] + # update filter, image_mask and detector + #self._filter = "D"+ self.subband_to_dichroic[self._band[1]] + #self._image_mask = "MIRI-IFU_" + self._band[0] + #self._update_detector() + #if not (self.MRSbands[self.band][0] <= self._wavelength <= self.MRSbands[self.band][1]): + # self._wavelength = np.mean(self.MRSbands[self.band]) + else: + raise ValueError(f"Not a valid MRS band: {value}") + class NIRCam(JWInstrument): """A class modeling the optics of NIRCam. From d342faf62e7a75b133f2484737e339111e06ad41 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Wed, 21 Feb 2024 00:36:27 -0500 Subject: [PATCH 04/12] switching MRS band should update detector attribute too --- webbpsf/tests/test_miri.py | 9 +++++++++ webbpsf/webbpsf_core.py | 1 + 2 files changed, 10 insertions(+) diff --git a/webbpsf/tests/test_miri.py b/webbpsf/tests/test_miri.py index 18080375..5f4ef520 100644 --- a/webbpsf/tests/test_miri.py +++ b/webbpsf/tests/test_miri.py @@ -159,3 +159,12 @@ def test_mode_switch(): assert miri.detector =='MIRIFUSHORT' assert miri.aperturename.startswith('MIRIFU_CH') assert miri._rotation != imager_rotation + + # band switching should toggle detector and aper name + miri.band = '4C' + assert miri.detector == 'MIRIFULONG' + assert miri.aperturename == 'MIRIFU_CHANNEL4C' + + miri.band = '2A' + assert miri.detector == 'MIRIFUSHORT' + assert miri.aperturename == 'MIRIFU_CHANNEL2A' diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py index e030b542..8119f327 100644 --- a/webbpsf/webbpsf_core.py +++ b/webbpsf/webbpsf_core.py @@ -2337,6 +2337,7 @@ def aperturename(self, value): if self.mode == 'IFU': if self.band != value[-2:]: self.band = value[-2:] + self._detector = 'MIRIFULONG' if self.band[0] in ['3', '4'] else 'MIRIFUSHORT' @property From f8fe3646568567f9fdca539d6d3843b46f73a33a Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Wed, 21 Feb 2024 00:53:12 -0500 Subject: [PATCH 05/12] get PSF calcs at least running, for MIRI IFU mode --- webbpsf/utils.py | 2 +- webbpsf/webbpsf_core.py | 22 +++++++++++----------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/webbpsf/utils.py b/webbpsf/utils.py index 2dfee5d8..644694d9 100644 --- a/webbpsf/utils.py +++ b/webbpsf/utils.py @@ -976,7 +976,7 @@ def determine_inst_name_from_v2v3(v2v3): _log.debug('Field coordinates determined to be in NIRSpec field') elif ( (v2v3[0] <= -6.2 * u.arcmin) - and (v2v3[0] >= -8.3 * u.arcmin) + and (v2v3[0] >= -8.5 * u.arcmin) and (v2v3[1] <= -5.2 * u.arcmin) and (v2v3[1] >= -7.3 * u.arcmin) ): diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py index 8119f327..1d88d3d0 100644 --- a/webbpsf/webbpsf_core.py +++ b/webbpsf/webbpsf_core.py @@ -1039,11 +1039,11 @@ def _tel_coords(self): if self._detector_geom_info.aperture.AperType == 'SLIT': # These apertures don't map directly to particular detector position in the usual way # Return coords for center of the aperture reference location - return ( - np.asarray((self._detector_geom_info.aperture.V2Ref, self._detector_geom_info.aperture.V3Ref)) - / 60 - * units.arcmin - ) + return np.asarray((self._detector_geom_info.aperture.V2Ref, + self._detector_geom_info.aperture.V3Ref)) / 60 * units.arcmin + elif self._detector_geom_info.aperture.AperType == 'COMPOUND': + # handle MIRI MRS apertures, which don't have V2Ref,V3Ref defined, but this works: + return np.asarray(self.siaf[self.aperturename].reference_point('tel') ) / 60 * units.arcmin else: return self._detector_geom_info.pix2angle(self.detector_position[0], self.detector_position[1]) @@ -1247,14 +1247,14 @@ def _calc_psf_format_output(self, result, options): elif self.name == 'MIRI': # Apply distortion effects to MIRI psf: Distortion and MIRI Scattering _log.debug('MIRI: Adding optical distortion and Si:As detector internal scattering') - if self.image_mask != 'LRS slit' and self.pupil_mask != 'P750L': + if self.mode != 'IFU': psf_siaf = distortion.apply_distortion(result) # apply siaf distortion + psf_siaf_rot = detectors.apply_miri_scattering(psf_siaf) # apply scattering effect + psf_distorted = detectors.apply_detector_charge_diffusion(psf_siaf_rot,options) # apply detector charge transfer model else: - psf_siaf = result # distortion model in SIAF not available for LRS - psf_siaf_rot = detectors.apply_miri_scattering(psf_siaf) # apply scattering effect (cruciform artifact) - psf_distorted = detectors.apply_detector_charge_diffusion( - psf_siaf_rot, options - ) # apply detector charge transfer model + # there is not yet any distortion calibration for the IFU, and + # we don't want to apply charge diffusion directly here + psf_distorted = result elif self.name == 'NIRSpec': # Apply distortion effects to NIRSpec psf: Distortion only # (because applying detector effects would only make sense after simulating spectral dispersion) From 7a5461f75fa3a0a72444acf64e88272d5cb660a9 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Wed, 21 Feb 2024 15:53:48 -0500 Subject: [PATCH 06/12] implement NIRSpec and MIRI IFUs consistently with band attributes and get_IFU_wavelengths functions --- webbpsf/__init__.py | 2 +- webbpsf/tests/test_nirspec.py | 9 ++ webbpsf/webbpsf_core.py | 154 +++++++++++++++++++++++++++++----- 3 files changed, 141 insertions(+), 24 deletions(-) diff --git a/webbpsf/__init__.py b/webbpsf/__init__.py index c202f4bb..100a5bfa 100644 --- a/webbpsf/__init__.py +++ b/webbpsf/__init__.py @@ -40,7 +40,7 @@ class UnsupportedPythonError(Exception): # required. If changes to the code and data mean WebbPSF won't work # properly with an old data package, increment this version number. # (It's checked against $WEBBPSF_DATA/version.txt) -DATA_VERSION_MIN = (1, 2, 1) +DATA_VERSION_MIN = (1, 3, 0) class Conf(_config.ConfigNamespace): diff --git a/webbpsf/tests/test_nirspec.py b/webbpsf/tests/test_nirspec.py index b07de169..0ff58194 100644 --- a/webbpsf/tests/test_nirspec.py +++ b/webbpsf/tests/test_nirspec.py @@ -56,7 +56,16 @@ def test_calc_datacube_fast(): def test_mode_switch(): nrs = webbpsf_core.NIRSpec() + # check mode swith to IFU nrs.mode = 'IFU' assert 'IFU' in nrs.aperturename + assert nrs.band == 'PRISM/CLEAR' + + # check switch of which IFU band + nrs.grating = 'G395H' + nrs.filter = 'F290LP' + assert nrs.band == 'G395H/F290LP' + + # check mode switch back to imaging nrs.mode = 'imaging' assert 'IFU' not in nrs.aperturename diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py index 1d88d3d0..ba416ba3 100644 --- a/webbpsf/webbpsf_core.py +++ b/webbpsf/webbpsf_core.py @@ -23,6 +23,7 @@ import glob from collections import namedtuple, OrderedDict import functools +from abc import ABC, abstractmethod import numpy as np import scipy.interpolate, scipy.ndimage @@ -1825,6 +1826,10 @@ def calc_datacube_fast(self, wavelengths, compare_methods=False, *args, **kwargs # Allow up to 10,000 wavelength slices. The number matters because FITS # header keys can only have up to 8 characters. Backward-compatible. + + if isinstance(wavelengths, units.Quantity): + wavelengths = wavelengths.to_value(units.m) + nwavelengths = len(wavelengths) if nwavelengths < 100: label_wl = lambda i: 'WAVELN{:02d}'.format(i) @@ -1915,7 +1920,7 @@ def calc_datacube_fast(self, wavelengths, compare_methods=False, *args, **kwargs return cubefast -class JWInstrument_with_IFU(JWInstrument): +class JWInstrument_with_IFU(JWInstrument, ABC): """ Subclass which adds some additional infrastructure for IFU sims""" @@ -1926,6 +1931,8 @@ def __init__(self, *args, **kwargs): 'IFU': None} self._mode = 'imaging' + self._IFU_bands_cubepars = {} # placeholder, subclass should implement + @property def mode(self): """Currently selected instrument major mode, imaging or IFU""" @@ -1939,6 +1946,26 @@ def mode(self, value): self._mode = value self.set_position_from_aperture_name(self._modes_list[value]) + @property + @abstractmethod + def band(self): + # Subclass must implement this + pass + + def get_IFU_wavelengths(self, nlambda=None): + """Return an array of wavelengths spanning the currently selected IFU sub-band + """ + if self.mode != 'IFU': + raise RuntimeError("This method only applies in IFU mode") + spaxelsize, wavestep, minwave, maxwave = self._IFU_bands_cubepars[self.band] + if nlambda: + # Return the specified number of wavelengths, across that band + return np.linspace(minwave, maxwave, nlambda) * units.micron + else: + # Return wavelength across that band using the same spectral sampling + # as the instrument and pipeline + return np.arange(minwave, maxwave, wavestep) * units.micron + class MIRI(JWInstrument_with_IFU): """ A class modeling the optics of MIRI, the Mid-InfraRed Instrument. @@ -1988,32 +2015,48 @@ def __init__(self): self.monochromatic = 8.0 self._IFU_pixelscale = { - 'Ch1': (0.18, 0.19), - 'Ch2': (0.28, 0.19), - 'Ch3': (0.39, 0.24), - 'Ch4': (0.64, 0.27), + 'Ch1': (0.176, 0.196), + 'Ch2': (0.277, 0.196), + 'Ch3': (0.387, 0.245), + 'Ch4': (0.645, 0.273), } - # The above tuples give the pixel resolution (perpendicular to the slice, along the slice). - # The pixels are not square. + # The above tuples give the pixel resolution (first the 'alpha' direction, perpendicular to the slice, + # then the 'beta' direction, along the slice). + # The pixels are not square. See https://jwst-docs.stsci.edu/jwst-mid-infrared-instrument/miri-observing-modes/miri-medium-resolution-spectroscopy # Mappings between alternate names used for MRS subbands self._MRS_dichroic_to_subband = {"SHORT": "A", "MEDIUM": "B", "LONG": "C"} self._MRS_subband_to_dichroic = {"A": "SHORT", "B": "MEDIUM", "C": "LONG"} self._band = None - self._MRS_bands = {"1A": [4.887326748103221, 5.753418963216559], # Values provided by Polychronis Patapis - "1B": [5.644625711181792, 6.644794583147869], # To-do: obtain from CRDS pipeline refs - "1C": [6.513777066360325, 7.669147994055998], - "2A": [7.494966046398437, 8.782517027772244], - "2B": [8.651469658142522, 10.168811217793243], - "2C": [9.995281242621394, 11.73039280033565], - "3A": [11.529088518317131, 13.491500288051483], - "3B": [13.272122736770127, 15.550153182343314], - "3C": [15.389530615108631, 18.04357852656418], - "4A": [17.686540162850203, 20.973301482912323], - "4B": [20.671069749545193, 24.476094964546686], - "4C": [24.19608171436692, 28.64871057821349]} - +# self._MRS_bands = {"1A": [4.887326748103221, 5.753418963216559], # Values provided by Polychronis Patapis +# "1B": [5.644625711181792, 6.644794583147869], # To-do: obtain from CRDS pipeline refs +# "1C": [6.513777066360325, 7.669147994055998], +# "2A": [7.494966046398437, 8.782517027772244], +# "2B": [8.651469658142522, 10.168811217793243], +# "2C": [9.995281242621394, 11.73039280033565], +# "3A": [11.529088518317131, 13.491500288051483], +# "3B": [13.272122736770127, 15.550153182343314], +# "3C": [15.389530615108631, 18.04357852656418], +# "4A": [17.686540162850203, 20.973301482912323], +# "4B": [20.671069749545193, 24.476094964546686], +# "4C": [24.19608171436692, 28.64871057821349]} + self._IFU_bands_cubepars = { # pipeline data cube parameters + # Taken from ifucubepars_table in CRDS file 'jwst_miri_cubepar_0014.fits', current as of 2023 December + # Each tuple gives pipeline spaxelsize, spectralstep, wave_min, wave_max + '1A': (0.13, 0.0008, 4.90, 5.74), + '1B': (0.13, 0.0008, 5.66, 6.63), + '1C': (0.13, 0.0008, 6.53, 7.65), + '2A': (0.17, 0.0013, 7.51, 8.77), + '2B': (0.17, 0.0013, 8.67, 10.13), + '2C': (0.17, 0.0013, 10.01, 11.70), + '3A': (0.20, 0.0025, 11.55, 13.47), + '3B': (0.20, 0.0025, 13.34, 15.57), + '3C': (0.20, 0.0025, 15.41, 17.98), + '4A': (0.35, 0.0060, 17.70, 20.95), + '4B': (0.35, 0.0060, 20.69, 24.48), + '4C': (0.35, 0.0060, 24.40, 28.70), + } self._detectors = {'MIRIM': 'MIRIM_FULL', # Mapping from user-facing detector names to SIAF entries. 'MIRIFUSHORT': 'MIRIFU_CHANNEL1A', # only applicable in IFU mode @@ -2355,7 +2398,7 @@ def band(self, value): raise RuntimeError("The 'band' property is only valid for IFU mode simulations.") return - if value in self._MRS_bands.keys(): + if value in self._IFU_bands_cubepars.keys(): self._band = value #self._slice_width = self._IFU_pixelscale[f"Ch{self._band[0]}"][0] @@ -2371,6 +2414,17 @@ def band(self, value): else: raise ValueError(f"Not a valid MRS band: {value}") + def _calc_psf_format_output(self, result, options): + """ Format output HDUList. In particular, add some extra metadata if in IFU mode""" + super()._calc_psf_format_output(result, options) + if self.mode =='IFU': + n_exts = len(result) + for ext in np.arange(n_exts): + result[ext].header['MODE'] = ("IFU", 'This is a MIRI MRS IFU mode simulation') + result[ext].header['FILTER'] = ('MIRIFU_CHANNEL'+self.band, 'MIRI IFU sub-band simulated') + result[ext].header['BAND'] = (self.band, 'MIRI IFU sub-band simulated') + + class NIRCam(JWInstrument): """A class modeling the optics of NIRCam. @@ -3032,6 +3086,20 @@ def __init__(self): self.image_mask = 'MSA all open' self.pupil_mask = self.pupil_mask_list[-1] + self.disperser_list = ['PRISM', 'G140M', 'G140H', 'G235M', 'G235H', 'G394M', 'G395H'] + self._disperser = None + self._IFU_bands_cubepars = { + 'PRISM/CLEAR': (0.10, 0.0050, 0.60, 5.30), + 'G140M/F070LP': (0.10, 0.0006, 0.70, 1.27), + 'G140M/F100LP': (0.10, 0.0006, 0.97, 1.89), + 'G140H/F070LP': (0.10, 0.0002, 0.70, 1.27), + 'G140H/F100LP': (0.10, 0.0002, 0.97, 1.89), + 'G235M/F170LP': (0.10, 0.0011, 1.66, 3.17), + 'G235H/F170LP': (0.10, 0.0004, 1.66, 3.17), + 'G395M/F290LP': (0.10, 0.0018, 2.87, 5.27), + 'G395H/F290LP': (0.10, 0.0007, 2.87, 5.27), + } + det_list = ['NRS1', 'NRS2'] self._detectors = dict() for name in det_list: @@ -3148,14 +3216,15 @@ def aperturename(self, value): self.pixelscale = self._get_pixelscale_from_apername(detector_apername) _log.debug(f"Pixelscale updated to {self.pixelscale} based on average X+Y SciScale at SIAF aperture {self._aperturename}") - if 'IFU' in self.aperturename: self._mode = 'IFU' + if self._disperser is None: + self.disperser = 'PRISM' # Set some default spectral mode + self.filter = 'CLEAR' else: self._mode = 'imaging' # More to implement here later! - def _tel_coords(self): """ Convert from science frame coordinates to telescope frame coordinates using SIAF transformations. Returns (V2, V3) tuple, in arcminutes. @@ -3182,6 +3251,45 @@ def _get_pixelscale_from_apername(self, apername): else: return super()._get_pixelscale_from_apername(apername) + @property + def disperser(self): + """NIRSpec spectral dispersing element (grating or prism). + Only applies for IFU mode sims, currently; used to help set the + wavelength range to simulate + """ + if self.mode =='IFU': + return self._disperser + else: + return None + + @disperser.setter + def disperser(self, value): + if (value is None) or (value in self.disperser_list): + self._disperser = value + else: + raise RuntimeError(f'Not a valid NIRSpec disperser name: {value}') + + def _calc_psf_format_output(self, result, options): + """ Format output HDUList. In particular, add some extra metadata if in IFU mode""" + super()._calc_psf_format_output(result, options) + if self.mode == 'IFU': + n_exts = len(result) + for ext in np.arange(n_exts): + result[ext].header['MODE'] = ("IFU", 'This is a NIRSpec IFU mode simulation') + result[ext].header['GRATING'] = (self.disperser, 'Name of the grating (or prism) element simulated.') + + + @property + def band(self): + if self.mode != 'IFU': + return None + + return self.disperser + "/" + self.filter + + @band.setter + def band(self, value): + raise RuntimeError("This is a read-only property. Set grating and/or filter attributes instead.") + class NIRISS(JWInstrument): """A class modeling the optics of the Near-IR Imager and Slit Spectrograph From c160eca8ebd9c2cc4da2839e7da958b098528b30 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Wed, 21 Feb 2024 16:47:48 -0500 Subject: [PATCH 07/12] implement pixelscale switching for MRS sims --- webbpsf/tests/test_miri.py | 27 +++++++++++++++++++++++++++ webbpsf/tests/test_nirspec.py | 18 +++++++++++++++++- webbpsf/webbpsf_core.py | 13 +++++++++---- 3 files changed, 53 insertions(+), 5 deletions(-) diff --git a/webbpsf/tests/test_miri.py b/webbpsf/tests/test_miri.py index 5f4ef520..9f2e276f 100644 --- a/webbpsf/tests/test_miri.py +++ b/webbpsf/tests/test_miri.py @@ -122,8 +122,13 @@ def test_miri_nonsquare_detector(): assert miri.detector_position == (1023, 1031) def test_mode_switch(): + """Test switching between imaging and IFU modes, and switching IFU bands + Also checks this works to switch aperturenane, and conversely setting aperturename switches mode if needed. + Also checks that this automatically changes the rotation and pixelscale properties, as expected. + """ miri = webbpsf_core.MIRI() imager_rotation = miri._rotation + imager_pixelscale = miri.pixelscale # Explicitly switch mode to IFU miri.mode = 'IFU' @@ -131,12 +136,15 @@ def test_mode_switch(): assert miri.detector =='MIRIFUSHORT' assert miri.aperturename.startswith('MIRIFU_CH') assert miri._rotation != imager_rotation + assert miri.pixelscale > imager_pixelscale # Explicitly switch back to imaging miri.mode = 'imaging' assert 'IFU' not in miri.aperturename assert miri.detector =='MIRIM' assert miri.aperturename.startswith('MIRIM_') assert miri._rotation == imager_rotation + assert miri.pixelscale == imager_pixelscale + # Implicitly switch to IFU miri.set_position_from_aperture_name('MIRIFU_CHANNEL3B') @@ -144,6 +152,7 @@ def test_mode_switch(): assert miri.detector =='MIRIFULONG' assert miri.aperturename == 'MIRIFU_CHANNEL3B' assert miri._rotation != imager_rotation + assert miri.pixelscale > imager_pixelscale # implicitly switch to imaging # LRS is an odd case, SLIT aper type but operates like in imaging mode @@ -152,6 +161,7 @@ def test_mode_switch(): assert miri.detector =='MIRIM' assert miri.aperturename.startswith('MIRIM_') assert miri._rotation == imager_rotation + assert miri.pixelscale == imager_pixelscale # And back to IFU again: miri.mode = 'IFU' @@ -159,12 +169,29 @@ def test_mode_switch(): assert miri.detector =='MIRIFUSHORT' assert miri.aperturename.startswith('MIRIFU_CH') assert miri._rotation != imager_rotation + assert miri.pixelscale > imager_pixelscale # band switching should toggle detector and aper name miri.band = '4C' assert miri.detector == 'MIRIFULONG' assert miri.aperturename == 'MIRIFU_CHANNEL4C' + assert miri.pixelscale > 3*imager_pixelscale miri.band = '2A' assert miri.detector == 'MIRIFUSHORT' assert miri.aperturename == 'MIRIFU_CHANNEL2A' + assert imager_pixelscale < miri.pixelscale < 2*imager_pixelscale + +def test_IFU_wavelengths(): + """ Test computing the wqvelength sampling for a sim IFU cube """ + miri = webbpsf_core.MIRI() + # check mode swith to IFU + miri.mode = 'IFU' + miri.band = '2A' + waves = miri.get_IFU_wavelengths() + assert isinstance(waves, u.Quantity) + + assert len(waves) > 900 # there are lots of wavelengths in MRScubes + # and test we can specify a reduced wavelength sampling: + for n in (10, 100): + assert len(miri.get_IFU_wavelengths(n)) == n diff --git a/webbpsf/tests/test_nirspec.py b/webbpsf/tests/test_nirspec.py index 0ff58194..17ab9deb 100644 --- a/webbpsf/tests/test_nirspec.py +++ b/webbpsf/tests/test_nirspec.py @@ -55,6 +55,7 @@ def test_calc_datacube_fast(): def test_mode_switch(): + """ Test switch between IFU and imaging modes """ nrs = webbpsf_core.NIRSpec() # check mode swith to IFU nrs.mode = 'IFU' @@ -62,10 +63,25 @@ def test_mode_switch(): assert nrs.band == 'PRISM/CLEAR' # check switch of which IFU band - nrs.grating = 'G395H' + nrs.disperser = 'G395H' nrs.filter = 'F290LP' assert nrs.band == 'G395H/F290LP' # check mode switch back to imaging nrs.mode = 'imaging' assert 'IFU' not in nrs.aperturename + +def test_IFU_wavelengths(): + """ Test computing the wqvelength sampling for a sim IFU cube """ + nrs = webbpsf_core.NIRSpec() + # check mode swith to IFU + nrs.mode = 'IFU' + nrs.disperser = 'G235H' + nrs.filter = 'F170LP' + waves = nrs.get_IFU_wavelengths() + assert isinstance(waves, u.Quantity) + + assert len(waves) > 3000 # there are lots of wavelengths in the high-resolution grating cubes + # and test we can specify a reduced wavelength sampling: + for n in (10, 100): + assert len(nrs.get_IFU_wavelengths(n)) == n diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py index ba416ba3..6065dd5c 100644 --- a/webbpsf/webbpsf_core.py +++ b/webbpsf/webbpsf_core.py @@ -996,6 +996,9 @@ def aperturename(self, value): self._aperturename = value # Update DetectorGeometry class self._detector_geom_info = DetectorGeometry(self.siaf, self._aperturename) + if not has_custom_pixelscale: + self.pixelscale = self._get_pixelscale_from_apername(value) + _log.debug( f"Pixelscale updated to {self.pixelscale} based on IFU cubepars for {value}") else: if self.detector not in value: @@ -2326,8 +2329,12 @@ def _get_pixelscale_from_apername(self, apername): """Simple utility function to look up pixelscale from apername""" if 'MIRIFU' in apername: - print('special case for MIRI IFU pixelscales') - return 0.1 + if apername.startswith('MIRIFU_CHANNEL'): + band = apername[-2:] + spaxelsize, _, _, _= self._IFU_bands_cubepars[band] + return spaxelsize + else: + raise RuntimeError(f"Not sure how to determine pixelscale for {apername}") else: return super()._get_pixelscale_from_apername(apername) @@ -2401,7 +2408,6 @@ def band(self, value): if value in self._IFU_bands_cubepars.keys(): self._band = value #self._slice_width = self._IFU_pixelscale[f"Ch{self._band[0]}"][0] - self.aperturename = 'MIRIFU_CHANNEL' + value # setting aperturename will also auto update self._rotation #self._rotation = self.MRS_rotation[self._band] @@ -3246,7 +3252,6 @@ def _tel_coords(self): def _get_pixelscale_from_apername(self, apername): """Simple utility function to look up pixelscale from apername""" if 'IFU' in apername: - print('DEBUG - special case for NIRSpec IFU pixelscales') return super()._get_pixelscale_from_apername('NRS1_FULL') else: return super()._get_pixelscale_from_apername(apername) From 4fbc614395c4338661e871e3eeb5c5c2c4959dc0 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Wed, 24 Apr 2024 15:03:59 -0400 Subject: [PATCH 08/12] add NIRSpec IFUalign rotation; start relnotes re IFU mode --- docs/relnotes.rst | 12 ++++++++++++ webbpsf/webbpsf_core.py | 18 +++++++++++++++++- 2 files changed, 29 insertions(+), 1 deletion(-) diff --git a/docs/relnotes.rst b/docs/relnotes.rst index 088c693a..1129039b 100644 --- a/docs/relnotes.rst +++ b/docs/relnotes.rst @@ -27,12 +27,24 @@ See https://github.com/spacetelescope/webbpsf/issues for currently open issues a Version History and Change Log ------------------------------- +Version 1.3.0 +============= +*Date TBD* + +This release comes with new features and improvements including but not limited to: + +X. Improved support for NIRSpec and MIRI IFU PSF calculations, including addition of a ``mode`` attribute for toggling between imaging mode and IFU mode simulations; an option for much faster (but slightly simplified) IFU datacube PSF calculations. Spectral bandpass information added for the IFU bands for both NIRSpec and MIRI. For NIRSpec, IFU mode PSF outputs are rotated by an additional 90 degrees to match the convention used in pipeline-output s3d datacubes made using the IFUAlign orientation. This extra rotation can be optionally disabled if so desired; see the NIRSpec class docstring. + + Version 1.2.1 ============= Minor documentation updates Version 1.2.0 ============= + +*2023 August* + We are pleased to announce the release of the latest version of WebbPSF version 1.2.0, now available on PyPi and GitHub. This release comes with new features and improvements including but not limited to: 1. The addition of detector effects for JWST simulations. H2RG detector effects are included in two flavors, a simple ad hoc Gaussian convolution to model charge diffusion effects and a set of convolution kernels to model interpixel capacitance (IPC) and post-pixel coupling effects. We have found that these effects greatly improve the agreement between observations and simulations. See `JWST Detector Effects for more details. `_ diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py index 6065dd5c..2a11286c 100644 --- a/webbpsf/webbpsf_core.py +++ b/webbpsf/webbpsf_core.py @@ -3055,7 +3055,18 @@ class NIRSpec(JWInstrument_with_IFU): is added to the optical system. This is an estimate of the pupil stop imposed by the outer edge of the grating clear aperture, estimated based on optical modeling by Erin Elliot and Marshall Perrin. - **Note: IFU to be implemented later** + Notes on IFU support: + Additional features for modeling NRS IFU PSFs are enabled by setting the .mode attribute to 'IFU'. + + The pipeline-output data products, assuming the 'ifualign' frame is used in the cube build step, which + is rotated relative to the typical 'sci' output frame used in all other webbpsf sim outputs. + For convenience, for IFU-mode simulations an extra rotation is included in the PSF calculation + such that the output product orientation matches the IFUalign s3d cube orientation. This happens + automatically and transparently to the user (and source offset parameters, e.g. options['source_offset_x'] + will automatically be interpreted as X and Y position in that output frame, including the effects of the + rotation). If the rotation to ifualign frame is for some reason not desired, it can be disabled by + setting nrs.options['ifualign_rotation'] = False + """ def __init__(self): @@ -3166,6 +3177,11 @@ def _addAdditionalOptics(self, optsys, oversample=2): if self.include_si_wfe: optsys.add_pupil(optic=self._si_wfe_class(self, where='spectrograph')) + if self.mode == 'IFU' and self.options.get('ifualign_rotation',True): + optsys.add_rotation(90, hide=True) # Rotate by 90 degrees clockwise to match the IFUalign output convention, with slices horizontal. + optsys.planes[-1].wavefront_display_hint = 'intensity' + + return (optsys, trySAM, SAM_box_size) def _get_fits_header(self, hdulist, options): From 569d51e25992d14c907f0cc0f7a898143ab33d06 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Sun, 12 May 2024 09:57:25 -0400 Subject: [PATCH 09/12] handle LRS edge case as part of this PR too; for compatibility with changes in PR #787 --- webbpsf/webbpsf_core.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py index 2a11286c..ca7b5ba6 100644 --- a/webbpsf/webbpsf_core.py +++ b/webbpsf/webbpsf_core.py @@ -1252,7 +1252,12 @@ def _calc_psf_format_output(self, result, options): # Apply distortion effects to MIRI psf: Distortion and MIRI Scattering _log.debug('MIRI: Adding optical distortion and Si:As detector internal scattering') if self.mode != 'IFU': - psf_siaf = distortion.apply_distortion(result) # apply siaf distortion + if self._detector_geom_info.aperture.AperType != 'SLIT': + psf_siaf = distortion.apply_distortion(result) # apply siaf distortion + else: + # slit type aperture, specifically LRS SLIT, does not have distortion polynomials + # therefore omit apply_distortion if a SLIT aperture is selected. + psf_siaf = result psf_siaf_rot = detectors.apply_miri_scattering(psf_siaf) # apply scattering effect psf_distorted = detectors.apply_detector_charge_diffusion(psf_siaf_rot,options) # apply detector charge transfer model else: From c69245325dc7479de0bf9dc26773b93ec6b529bb Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Mon, 13 May 2024 15:13:01 -0400 Subject: [PATCH 10/12] update minimal-webbpsf-data to version 1.3.0 --- .github/workflows/ci_workflows.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci_workflows.yml b/.github/workflows/ci_workflows.yml index 198be84f..b825a6f8 100644 --- a/.github/workflows/ci_workflows.yml +++ b/.github/workflows/ci_workflows.yml @@ -58,7 +58,7 @@ jobs: - name: Get WebbPSF Data run: | # Get WebbPSF data files (just a subset of the full dataset!) and set up environment variable - wget https://stsci.box.com/shared/static/0ojjfg3cieqdpd18vl1bjnpe63r82dx8.gz -O /tmp/minimal-webbpsf-data.tar.gz + wget https://stsci.box.com/shared/static/gkbnn3jwcootr9riwv0rv1t8rpvxl3vr.gz -O /tmp/minimal-webbpsf-data.tar.gz tar -xzvf /tmp/minimal-webbpsf-data.tar.gz echo "WEBBPSF_PATH=${{github.workspace}}/webbpsf-data" >> $GITHUB_ENV From a6804b13207c68386ec15fa71c9a88eb540201fd Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Wed, 15 May 2024 00:49:50 -0400 Subject: [PATCH 11/12] implement outfile parameter for calc_datacube_fast --- webbpsf/webbpsf_core.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py index b48b2b89..441714cd 100644 --- a/webbpsf/webbpsf_core.py +++ b/webbpsf/webbpsf_core.py @@ -1784,7 +1784,7 @@ def load_wss_opd_by_date(self, date=None, choice='closest', verbose=True, plot=F opd_fn = webbpsf.mast_wss.get_opd_at_time(date, verbose=verbose, choice=choice, **kwargs) self.load_wss_opd(opd_fn, verbose=verbose, plot=plot, **kwargs) - def calc_datacube_fast(self, wavelengths, compare_methods=False, *args, **kwargs): + def calc_datacube_fast(self, wavelengths, compare_methods=False, outfile=None, *args, **kwargs): """Calculate a spectral datacube of PSFs: Simplified, much MUCH faster version. This is adapted from poppy.Instrument.calc_datacube, optimized and simplified @@ -1920,6 +1920,11 @@ def calc_datacube_fast(self, wavelengths, compare_methods=False, *args, **kwargs return cube, cubefast, waves, waves2 # return extra stuff for compariosns + if outfile is not None: + cubefast[0].header["FILENAME"] = (os.path.basename(outfile), "Name of this file") + cubefast.writeto(outfile, overwrite=True) + _log.info("Saved result to " + outfile) + return cubefast class JWInstrument_with_IFU(JWInstrument, ABC): From a8416a13a8dc1fc8b115d1604a0fc7a3b2ad6457 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Wed, 15 May 2024 13:30:02 -0400 Subject: [PATCH 12/12] for calc_datacube_fast, improve choice of ref wavelength; and clean up some redundant FITS header keywords --- webbpsf/webbpsf_core.py | 29 ++++++++++++++++++++++++----- 1 file changed, 24 insertions(+), 5 deletions(-) diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py index 441714cd..ef26df11 100644 --- a/webbpsf/webbpsf_core.py +++ b/webbpsf/webbpsf_core.py @@ -1797,12 +1797,15 @@ def calc_datacube_fast(self, wavelengths, compare_methods=False, outfile=None, * many times changing only the wavelength for just the last DFT step to the detector. 2) Assumes we do not need the binned-to-detector-resolution nor distorted versions; we just want the oversampled PSF datacube at many wavelengths as fast as possible. + (If the binned output is also desired, it can be computed post facto. + TODO: A future revision of this function may also add here an option for computing those + derived versions as well.) Testing for NIRSpec IFU indicates this achieves ~150x speedup, and the differences in computed oversampled PSF are typically ~1/100th or less relative to the local PSF values in any given pixel. - A consequence of the above iassumption 1 is that this methiod is not well applicable + A consequence of the above assumption 1 is that this method is not well applicable for cases that have image plane masks, nor for NIRCam in general. It does seem to be reasonably applicable for NIRSpec IFU calculations within the current limited fidelity of webbpsf for that mode, IF we also neglect the image plane stop around the IFU FOV. @@ -1840,16 +1843,30 @@ def calc_datacube_fast(self, wavelengths, compare_methods=False, outfile=None, * else: raise ValueError('Maximum number of wavelengths exceeded. ' 'Cannot be more than 10,000.') - # Set up cube and initialize structure based on PSF at first wavelength - poppy.poppy_core._log.info('Starting multiwavelength data cube calculation.') - REF_WAVE = 2e-6 # This must not be too short, to avoid phase wrapping for the C3 bump + # Set up cube and initialize structure based on PSF at a representative wavelength + _log.info('Starting fast/simplified multiwavelength data cube calculation.') + ref_wave = np.mean(wavelengths) + MIN_REF_WAVE = 2e-6 # This must not be too short, to avoid phase wrapping for the C3 bump + if ref_wave < MIN_REF_WAVE: + ref_wave = MIN_REF_WAVE + _log.info(f"Performing initial propagation at minimum wavelength {MIN_REF_WAVE*1e6:.2f} microns; minimum set to avoid phase wrap of segment C3 surface.") + else: + _log.info(f"Performing initial propagation at average wavelength {ref_wave*1e6:.2f} microns.") - psf, waves = self.calc_psf(*args, monochromatic=REF_WAVE, return_intermediates=True, **kwargs) + psf, waves = self.calc_psf(*args, monochromatic=ref_wave, return_intermediates=True, **kwargs) from copy import deepcopy # Setup arrays to save data # Copy the first (oversampled) HDU only cubefast = astropy.io.fits.HDUList(deepcopy(psf[0])) + try: + # This is cosmetic only. Delete some exteraneous/redundant header keywords from poppy. + # This function will below add a complete set of wavelength keywords + del cubefast[0].header['WAVE0'] + del cubefast[0].header['WGHT0'] + except KeyError: + pass + ext = 0 cubefast[ext].data = np.zeros((nwavelengths, psf[ext].data.shape[0], psf[ext].data.shape[1])) cubefast[ext].data[0] = psf[ext].data @@ -1884,6 +1901,8 @@ def calc_datacube_fast(self, wavelengths, compare_methods=False, outfile=None, * wl = wavelengths[i] psfw = quickosys.calc_psf(wavelength=wl, normalize='None') cubefast[0].data[i] = psfw[0].data + cubefast[ext].header[label_wl(i)] = wavelengths[i] + cubefast[0].header['NWAVES'] = nwavelengths ### OPTIONAL