-
Notifications
You must be signed in to change notification settings - Fork 58
/
webbpsf_core.py
2880 lines (2350 loc) · 138 KB
/
webbpsf_core.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
"""
============
WebbPSF Core
============
An object-oriented modeling system for the JWST instruments.
Classes:
* SpaceTelescopeInstrument
* JWInstrument
* MIRI
* NIRCam
* NIRSpec
* NIRISS
* FGS
WebbPSF makes use of python's ``logging`` facility for log messages, using
the logger name "webbpsf".
Code by Marshall Perrin <mperrin@stsci.edu>
"""
import os
import glob
from collections import namedtuple, OrderedDict
import numpy as np
import scipy.interpolate, scipy.ndimage
import astropy
import astropy.io.fits as fits
import astropy.io.ascii as ioascii
import astropy.units as units
import poppy
import pysiaf
from . import conf
from . import utils
from . import optics
from . import DATA_VERSION_MIN
from . import detectors
from . import distortion
from . import gridded_library
from . import opds
from . import constants
import webbpsf.mast_wss
try:
from .version import version
except ImportError:
version = ''
try:
_HAS_SYNPHOT = poppy.instrument._HAS_SYNPHOT
except AttributeError:
_HAS_SYNPHOT = False
if _HAS_SYNPHOT:
import synphot
import logging
_log = logging.getLogger('webbpsf')
Filter = namedtuple('Filter', ['name', 'filename', 'default_nlambda'])
class SpaceTelescopeInstrument(poppy.instrument.Instrument):
""" A generic Space Telescope Instrument class.
*Note*: Do not use this class directly - instead use one of the specific instrument subclasses!
This class provides a simple interface for modeling PSF formation through the instrument,
with configuration options and software interface loosely resembling the configuration of the instrument
hardware mechanisms.
This module currently only provides a modicum of error checking, and relies on the user
being knowledgable enough to avoid trying to simulate some physically impossible or just plain silly
configuration (such as trying to use a FQPM with the wrong filter).
The instrument constructors do not take any arguments. Instead, create an instrument object and then
configure the `filter` or other attributes as desired. The most commonly accessed parameters are
available as object attributes: `filter`, `image_mask`, `pupil_mask`, `pupilopd`. More advanced
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
may or may not be meaningful depending on which instrument is selected.
This is a superset of the options provided in :py:attr:`poppy.Instrument.options`.
Parameters
----------
source_offset_r : float
Radial offset of the target from the center, in arcseconds
source_offset_theta : float
Position angle for that offset, in degrees CCW.
pupil_shift_x, pupil_shift_y : float
Relative shift of the intermediate (coronagraphic) pupil in X and Y
relative to the telescope entrance pupil, expressed as a decimal between -1.0-1.0
Note that shifting an array too much will wrap around to the other side unphysically, but
for reasonable values of shift this is a non-issue. This option only has an effect for optical models that
have something at an intermediate pupil plane between the telescope aperture and the detector.
pupil_rotation : float
Relative rotation of the intermediate (coronagraphic) pupil relative to
the telescope entrance pupil, expressed in degrees counterclockwise.
This option only has an effect for optical models that have something at
an intermediate pupil plane between the telescope aperture and the detector.
rebin : bool
For output files, write an additional FITS extension including a version of the output array
rebinned down to the actual detector pixel scale?
jitter : string "gaussian" or None
Type of jitter model to apply. Currently only convolution with a Gaussian kernel of specified
width `jitter_sigma` is implemented. (default: None)
jitter_sigma : float
Width of the jitter kernel in arcseconds (default: 0.006 arcsec, 1 sigma per axis)
parity : string "even" or "odd"
You may wish to ensure that the output PSF grid has either an odd or even number of pixels.
Setting this option will force that to be the case by increasing npix by one if necessary.
Note that this applies to the number detector pixels, rather than the subsampled pixels if oversample > 1.
force_coron : bool
Set this to force full coronagraphic optical propagation when it might not otherwise take place
(e.g. calculate the non-coronagraphic images via explicit propagation to all optical surfaces, FFTing
to intermediate pupil and image planes whether or not they contain any actual optics, rather than
taking the straight-to-MFT shortcut)
no_sam : bool
Set this to prevent the SemiAnalyticMethod coronagraph mode from being
used when possible, and instead do the brute-force FFT calculations.
This is usually not what you want to do, but is available for comparison tests.
The SAM code will in general be much faster than the FFT method,
particularly for high oversampling.
"""
_detectors = {}
"""
Dictionary mapping detector names to detector or wavefront information in some manner.
The specific meaning of this mapping must be defined by subclasses as part of their
implementation.
(Subclasses must populate this at `__init__`.)
"""
_detector = None
"""
The name of the currently selected detector. Must be a key in _detectors, as validated by the
`detector` property setter.
(Subclasses must populate this at `__init__`.)
"""
def _get_filters(self):
filter_table = ioascii.read(os.path.join(self._WebbPSF_basepath, self.name, 'filters.tsv'))
filter_info = {}
filter_list = [] # preserve the order from the table
for filter_row in filter_table:
filter_filename = os.path.join(
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']
)
filter_list.append(filter_row['filter'])
return filter_list, filter_info
def _get_default_nlambda(self, filtername):
""" Return the default # of wavelengths to be used for calculation by a given filter """
return self._filters[filtername].default_nlambda
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)
self._datapath = os.path.join(self._WebbPSF_basepath, self.name)
self._image_mask = None
self._pupil_mask = None
self.pupil = None
"Filename *or* fits.HDUList for the pupil mask."
self.pupilopd = None # This can optionally be set to a tuple indicating (filename, slice in datacube)
"""Filename *or* fits.HDUList for pupil OPD.
This can be either a full absolute filename, or a relative name in which case it is
assumed to be within the instrument's `data/OPDs/` directory, or an actual
fits.HDUList object corresponding to such a file. If the file contains a
datacube, you may set this to a tuple (filename, slice) to select a
given slice, or else the first slice will be used."""
self.pupil_radius = None # Set when loading FITS file in get_optical_system
self.options = {} # dict for storing other arbitrary options.
# filter_list available filter names in order by wavelength for public api
# _filters a dict of named tuples with name, filename, & default_nlambda
# with the filter name as the key
self.filter_list, self._filters = self._get_filters()
# choose a default filter, in case the user doesn't specify one
self.filter = self.filter_list[0]
self._rotation = None
self._image_mask = None
self.image_mask_list = []
"List of available image_masks"
self._pupil_mask = None
self.pupil_mask_list = []
"List of available pupil_masks"
self.pixelscale = pixelscale
"Detector pixel scale, in arcsec/pixel"
self._spectra_cache = {} # for caching synphot results.
# n.b.STInstrument subclasses must set these
self._detectors = {}
self._detector = None
self._detector_npixels = 2048
@property
def image_mask(self):
"""Currently selected image plane mask, or None for direct imaging"""
return self._image_mask
@image_mask.setter
def image_mask(self, name):
if name == "": name = None
if name is not None:
if name in self.image_mask_list:
pass # there's a perfect match, this is fine.
else:
name = name.upper() # force to uppercase
if name not in self.image_mask_list: # if still not found, that's an error.
raise ValueError("Instrument %s doesn't have an image mask called '%s'." % (self.name, name))
self._image_mask = name
if hasattr(self, '_image_mask_apertures') and name in self._image_mask_apertures:
self.set_position_from_aperture_name(self._image_mask_apertures[name])
@property
def pupil_mask(self):
"""Currently selected Lyot pupil mask, or None for direct imaging"""
return self._pupil_mask
@pupil_mask.setter
def pupil_mask(self, name):
if name == "":
name = None
if name is not None:
if name in self.pupil_mask_list:
pass # there's a perfect match, this is fine.
else:
name = name.upper() # force to uppercase
if name not in self.pupil_mask_list:
raise ValueError("Instrument %s doesn't have a pupil mask called '%s'." % (self.name, name))
self._pupil_mask = name
def __str__(self):
return "<{telescope}: {instrument_name}>".format(telescope=self.telescope, instrument_name=self.name)
@property
def detector(self):
"""Detector selected for simulated PSF
Used in calculation of field-dependent aberrations. Must be
selected from detectors in the `detector_list` attribute.
"""
return self._detector
@detector.setter
def detector(self, value):
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()
self._update_aperturename() # automatically set an appropriate aperture name
@property
def detector_list(self):
"""Detectors on which the simulated PSF could lie"""
return sorted(self._detectors.keys())
@property
def detector_position(self):
"""The pixel position in (X, Y) on the detector, relative to the currently-selected SIAF aperture subarray.
By default the SIAF aperture will correspond to the full-frame detector, so (X,Y) will in that case be
absolute (X,Y) pixels on the detector. But if you select a subarray aperture name from the SIAF, then
the (X,Y) are interpreted as (X,Y) within that subarray.
Please note, this is X,Y order - **not** a Pythonic y,x axes ordering.
"""
return self._detector_position
@detector_position.setter
def detector_position(self, position):
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._detector_position = (int(position[0]), int(position[1]))
@property
def aperturename(self):
""" SIAF aperture name for detector pixel to sky coords transformations"""
return self._aperturename
@aperturename.setter
def aperturename(self, value):
# Override in subclass to provide more specific functionality
self._aperturename = value
def _update_aperturename(self):
""" Update SIAF aperture name after change in detector or other relevant properties
"""
self.aperturename = self._detectors[self._detector]
def _get_fits_header(self, result, options):
""" populate FITS Header keywords """
super(SpaceTelescopeInstrument, self)._get_fits_header(result, options)
result[0].header['FILTER'] = (self.filter, 'Filter name')
if self.image_mask is not None:
result[0].header['CORONMSK'] = (self.image_mask, "Image plane mask")
if self.pupil_mask is not None:
result[0].header['PUPIL'] = (self.pupil_mask, "Pupil plane mask")
result[0].header['VERSION'] = (version, "WebbPSF software version")
result[0].header['DATAVERS'] = (self._data_version, "WebbPSF reference data files version")
result[0].header['DET_NAME'] = (self.detector, "Name of detector on this instrument")
# Correct detector pixel coordinates to allow for even arrays to be centered on half pixel boundary
dpos = np.asarray(self.detector_position, dtype=float)
oversamp = result[0].header['OVERSAMP']
size = result[0].data.shape[0]
if size / oversamp % 2 == 0: dpos += 0.5 # even arrays must be at a half pixel
result[0].header['DET_X'] = (dpos[0], "Detector X pixel position of array center")
result[0].header['DET_Y'] = (dpos[1], "Detector Y pixel position of array center")
for key in self._extra_keywords:
result[0].header[key] = self._extra_keywords[key]
def get_optical_system(self, fft_oversample=2, detector_oversample=None,
fov_arcsec=2, fov_pixels=None, options=None):
""" Return an OpticalSystem instance corresponding to the instrument as currently configured.
When creating such an OpticalSystem, you must specify the parameters needed to define the
desired sampling, specifically the oversampling and field of view.
Parameters
----------
fft_oversample : int
Oversampling factor for intermediate plane calculations. Default is 2
detector_oversample: int, optional
By default the detector oversampling is equal to the intermediate calculation oversampling.
If you wish to use a different value for the detector, set this parameter.
Note that if you just want images at detector pixel resolution you will achieve higher fidelity
by still using some oversampling (i.e. *not* setting `oversample_detector=1`) and instead rebinning
down the oversampled data.
fov_pixels : float
Field of view in pixels. Overrides fov_arcsec if both set.
fov_arcsec : float
Field of view, in arcseconds. Default is 2
Returns
-------
osys : poppy.OpticalSystem
an optical system instance representing the desired configuration.
"""
_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.
if options is None: options = self.options
if detector_oversample is None: detector_oversample = fft_oversample
_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
)
# For convenience offsets can be given in cartesian or radial coords
if 'source_offset_x' in options or 'source_offset_y' in options:
offx = options.get('source_offset_x', 0)
offy = options.get('source_offset_y', 0)
optsys.source_offset_r = np.sqrt(offx ** 2 + offy ** 2)
optsys.source_offset_theta = np.rad2deg(np.arctan2(-offx, offy))
_log.debug("Source offset from X,Y = ({}, {}) is (r,theta) = {},{}".format(
offx, offy, optsys.source_offset_r, optsys.source_offset_theta))
if 'source_offset_r' in options:
optsys.source_offset_r = options['source_offset_r']
if 'source_offset_theta' in options:
optsys.source_offset_theta = options['source_offset_theta']
# Telescope entrance pupil
pupil_optic = self._get_telescope_pupil_and_aberrations()
optsys.add_pupil(pupil_optic)
pupil_rms_wfe_nm = np.sqrt(np.mean(pupil_optic.opd[pupil_optic.amplitude == 1] ** 2)) * 1e9
self._extra_keywords['TEL_WFE'] = (pupil_rms_wfe_nm, '[nm] Telescope pupil RMS wavefront error')
if hasattr(pupil_optic, 'header_keywords'):
self._extra_keywords.update(pupil_optic.header_keywords())
self.pupil_radius = pupil_optic.pupil_diam / 2.0
# add coord transform from entrance pupil to exit pupil
optsys.add_inversion(axis='y', name='OTE exit pupil', hide=True)
# add rotation at this point, if present - needs to be after the
# exit pupil inversion.
# Sign convention: Here we are rotating the *wavefront* so the sign is opposite the _rotation attribute,
# which gives the V3IdlYangle for the detector rotation.
if self._rotation is not None:
optsys.add_rotation(-self._rotation, hide=True)
optsys.planes[-1].wavefront_display_hint = 'intensity'
# Allow instrument subclass to add field-dependent aberrations
aberration_optic = self._get_aberrations()
if aberration_optic is not None:
optsys.add_pupil(aberration_optic)
try:
# Calculate SI WFE over just the OTE entrance pupil aperture,
# though with a flip in the Y axis to account for entrance vs. exit pupil conventions
exit_pupil_mask = pupil_optic.amplitude[::-1] == 1
inst_rms_wfe_nm = np.sqrt(np.mean(aberration_optic.opd[exit_pupil_mask] ** 2)) * 1e9
self._extra_keywords['SI_WFE'] = (inst_rms_wfe_nm, '[nm] instrument pupil RMS wavefront error')
except (TypeError, IndexError):
# Currently the above does not work for Roman, but fixing this is deferred to future work
pass
if hasattr(aberration_optic, 'header_keywords'):
self._extra_keywords.update(aberration_optic.header_keywords())
# ---- Add defocus if requested
if 'defocus_waves' in options:
defocus_waves = options['defocus_waves']
defocus_wavelength = float(options['defocus_wavelength']) if 'defocus_wavelength' in options else 2.0e-6
_log.info(f"Adding defocus of {defocus_waves:.3f} waves at {defocus_wavelength*1e6:.3f} microns" )
lens = poppy.ThinLens(
name='Defocus',
nwaves=defocus_waves,
reference_wavelength=defocus_wavelength,
radius=self.pupil_radius
)
optsys.add_pupil(optic=lens)
self._extra_keywords['DEFOCUS'] = (defocus_waves, '# of waves of defocus added')
self._extra_keywords['DEFOC_WL'] = (defocus_wavelength, 'Wavelength reference for defocus added')
# ---- add coronagraph or spectrograph optics if requested,
# and possibly flag to invoke semi-analytic coronagraphic propagation
# first error check for null strings, which should be considered like None
if self.image_mask == "": self.image_mask = None
if self.pupil_mask == "": self.pupil_mask = None
if (self.image_mask is not None or self.pupil_mask is not None or
'WL' in self.filter or # special case handling for NIRCam WLP4 filter that is also a lens
('force_coron' in options and options['force_coron'])):
_log.debug("Adding coronagraph/spectrograph optics...")
optsys, trySAM, SAM_box_size = self._addAdditionalOptics(optsys, oversample=fft_oversample)
else:
trySAM = False
# --- add the detector element.
if fov_pixels is None:
if not np.isscalar(fov_arcsec): fov_arcsec = np.asarray(fov_arcsec) # cast to ndarray if 2D
fov_pixels = np.round(fov_arcsec / self.pixelscale)
if 'parity' in options:
if options['parity'].lower() == 'odd' and np.remainder(fov_pixels, 2) == 0: fov_pixels += 1
if options['parity'].lower() == 'even' and np.remainder(fov_pixels, 2) == 1: fov_pixels += 1
else:
pass
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")
return SAM_optsys
except ValueError as err:
_log.warning(
"Could not switch to Semi-Analytic Coronagraphy mode; invalid set of optical planes? "
"Using default propagation instead.")
_log.warning(str(err))
# _log.warn("ERROR ({0}): {1}".format(errno, strerror))
pass
return optsys
def _get_telescope_pupil_and_aberrations(self):
"""return OpticalElement modeling wavefront aberrations for the telescope.
See also get_aberrations for the SI aberrations.
"""
# ---- set pupil OPD
if isinstance(self.pupilopd, str): # simple filename
opd_map = self.pupilopd if os.path.exists(self.pupilopd) else \
os.path.join(self._datapath, "OPD", self.pupilopd)
elif hasattr(self.pupilopd, '__getitem__') and isinstance(self.pupilopd[0], str):
# tuple with filename and slice
opd_map = (self.pupilopd[0] if os.path.exists(self.pupilopd[0])
else os.path.join(self._datapath, "OPD", self.pupilopd[0]),
self.pupilopd[1])
elif isinstance(self.pupilopd, (fits.HDUList, poppy.OpticalElement)):
opd_map = self.pupilopd # not a path per se but this works correctly to pass it to poppy
elif self.pupilopd is None:
opd_map = None
else:
raise TypeError("Not sure what to do with a pupilopd of that type:" + str(type(self.pupilopd)))
# ---- set pupil intensity
if self.pupil is None:
raise RuntimeError("The pupil shape must be specified in the "
"instrument class or by setting self.pupil")
if isinstance(self.pupil, poppy.OpticalElement):
# supply to POPPY as-is
pupil_optic = self.pupil
else:
# wrap in an optic and supply to POPPY
if isinstance(self.pupil, str): # simple filename
if os.path.exists(self.pupil):
pupil_transmission = self.pupil
else:
pupil_transmission = os.path.join(
self._WebbPSF_basepath,
self.pupil
)
elif isinstance(self.pupil, fits.HDUList):
# POPPY can use self.pupil as-is
pupil_transmission = self.pupil
else:
raise TypeError("Not sure what to do with a pupil of "
"that type: {}".format(type(self.pupil)))
# ---- apply pupil intensity and OPD to the optical model
pupil_optic = poppy.FITSOpticalElement(
name='{} Entrance Pupil'.format(self.telescope),
transmission=pupil_transmission,
opd=opd_map,
planetype=poppy.poppy_core.PlaneType.pupil
# rotation=self._rotation
)
return pupil_optic
def _addAdditionalOptics(self, optsys, oversample=2):
"""Add instrument-internal optics to an optical system, typically coronagraphic or
spectrographic in nature. This method must be provided by derived instrument classes.
Returns
--------
optsys : OpticalSystem
modified to add coronagraph optics
useSAM : bool
flag that, after adding the Detector, the whole thing should be converted to
a SemiAnalyticCoronagraph model
SAM_box_size : float
size of box that entirely encloses the image plane occulter, in arcsec.
"""
raise NotImplementedError("needs to be subclassed.")
def _get_synphot_bandpass(self, filtername):
""" Return a synphot.spectrum.SpectralElement object for the given desired band.
By subclassing this, you can define whatever custom bandpasses are appropriate for
your instrument
"""
# use our local throughput files and create a synphot
# transmission object.
try:
filter_info = self._filters[filtername]
except KeyError:
msg = "Couldn't find filter '{}' for {} in PySynphot or local throughput files"
raise RuntimeError(msg.format(filtername, self.name))
# The existing FITS files all have wavelength in ANGSTROMS since that is
# the pysynphot convention...
filterfits = fits.open(filter_info.filename)
waveunit = filterfits[1].header.get('WAVEUNIT')
if waveunit is None:
_log.warning('The supplied file, {}, does not have a WAVEUNIT keyword. Assuming it '
'is Angstroms.'.format(filter_info.filename))
waveunit = 'angstrom'
filterdata = filterfits[1].data
try:
band = synphot.SpectralElement(synphot.models.Empirical1D, points=filterdata.WAVELENGTH,
lookup_table=filterdata.THROUGHPUT, keep_neg=False)
except AttributeError:
raise ValueError("The supplied file, %s, does not appear to be a FITS table "
"with WAVELENGTH and THROUGHPUT columns." % filter_info.filename)
filterfits.close()
return band
def psf_grid(self, num_psfs=16, all_detectors=True, save=False,
outdir=None, outfile=None, overwrite=True, verbose=True,
use_detsampled_psf=False, single_psf_centered=True, **kwargs):
"""
Create a PSF library in the form of a grid of PSFs across the detector
based on the specified instrument, filter, and detector. The output
GriddedPSFModel object will contain a 3D array with axes [i, y, x]
where i is the PSF position on the detector grid and (y,x) is the 2D
PSF.
Parameters
----------
num_psfs : int
The total number of fiducial PSFs to be created and saved in the files.
This number must be a square number. Default is 16.
E.g. num_psfs = 16 will create a 4x4 grid of fiducial PSFs.
all_detectors : bool
If True, run all detectors for the instrument. If False, run for
the detector set in the instance. Default is True
save : bool
True/False boolean if you want to save your file. Default is False.
outdir : str
If "save" keyword is set to True, your file will be saved in the
specified directory. Default of None will save it in the current
directory
outfile : str
If "save" keyword is set to True, your file will be saved as
{outfile}_det.fits. Default of None will save it as
instr_det_filt_fovp#_samp#_npsf#.fits
overwrite : bool
True/False boolean to overwrite the output file if it already exists.
Default is True.
verbose : bool
True/False boolean to print status updates. Default is True.
use_detsampled_psf : bool
If True, the grid of PSFs returned will be detector sampled (made
by binning down the oversampled PSF). If False, the PSFs will be
oversampled by the factor defined by the
oversample/detector_oversample/fft_oversample keywords. Default is False.
This is rarely needed - if uncertain, leave this alone.
single_psf_centered : bool
If num_psfs is set to 1, this defines where that psf is located.
If True it will be the center of the detector, if False it will
be the location defined in the WebbPSF attribute detector_position
(reminder - detector_position is (x,y)). Default is True
This is also rarely needed.
**kwargs
Any extra arguments to pass the WebbPSF calc_psf() method call.
Returns
-------
gridmodel : photutils GriddedPSFModel object or list of objects
Returns a GriddedPSFModel object or a list of objects if more than one
configuration is specified (1 per instrument, detector, and filter)
User also has the option to save the grid as a fits.HDUlist object.
Use
----
nir = webbpsf.NIRCam()
nir.filter = "F090W"
list_of_grids = nir.psf_grid(all_detectors=True, num_psfs=4)
wfi = webbpsf.WFI()
wfi.filter = "Z087"
wfi.detector = "SCA02"
grid = wfi.psf_grid(all_detectors=False, oversample=5, fov_pixels=101)
"""
# Keywords that could be set before the method call
filt = self.filter
if all_detectors is True:
detectors = "all"
else:
detectors = self.detector
if single_psf_centered is True:
psf_location = (int((self._detector_npixels - 1) / 2), int((self._detector_npixels - 1) / 2)) # center pt
else:
psf_location = self.detector_position[::-1] # (y,x)
# Call CreatePSFLibrary class
inst = gridded_library.CreatePSFLibrary(instrument=self, filter_name=filt, detectors=detectors,
num_psfs=num_psfs, psf_location=psf_location,
use_detsampled_psf=use_detsampled_psf, save=save,
outdir=outdir, filename=outfile, overwrite=overwrite,
verbose=verbose, **kwargs)
gridmodel = inst.create_grid()
return gridmodel
####### JWInstrument classes #####
@utils.combine_docstrings
class JWInstrument(SpaceTelescopeInstrument):
""" Superclass for all JWST instruments
Notable attributes
-------------------
telescope : name of telescope
pupilopd : filename or FITS file object
include_si_wfe : boolean (default: True)
Should SI internal WFE be included in models? Requires
the presence of ``si_zernikes_isim_cv3.fits`` in the
``WEBBPSF_PATH``.
"""
telescope = "JWST"
pupilopd = None
"""Filename *or* fits.HDUList for JWST pupil OPD.
This can be either a full absolute filename, or a relative name in which
case it is assumed to be within the instrument's `data/OPDs/` directory,
or an actual fits.HDUList object corresponding to such a file. If the file
contains a datacube, you may set this to a tuple (filename, slice) to
select a given slice, or else the first slice will be used."""
def __init__(self, *args, **kwargs):
super(JWInstrument, self).__init__(*args, **kwargs)
self.siaf = pysiaf.Siaf(self.name)
opd_path = os.path.join(self._datapath, 'OPD')
self.opd_list = []
for filename in glob.glob(os.path.join(opd_path, 'OPD*.fits*')):
self.opd_list.append(os.path.basename(os.path.abspath(filename)))
for filename in glob.glob(os.path.join(self._WebbPSF_basepath, 'JWST_OTE_OPD*.fits*')):
self.opd_list.append(os.path.basename(os.path.abspath(filename)))
if not len(self.opd_list) > 0:
raise RuntimeError("No pupil OPD files found for {name} in {path}".format(name=self.name, path=opd_path))
self.opd_list.sort()
self.pupilopd = 'JWST_OTE_OPD_cycle1_example_2022-07-30.fits' # Default is now an on-orbit measured example OPD
self.pupil = os.path.abspath(os.path.join(
self._WebbPSF_basepath,
"jwst_pupil_RevW_npix1024.fits.gz"
))
"Filename *or* fits.HDUList for JWST pupil mask. Usually there is no need to change this."
self._aperturename = None
self._detector = None
# where is the source on the detector, in 'Science frame' pixels?
self.detector_position = (1024, 1024)
self.include_si_wfe = True
self.include_ote_field_dependence = True # Note, this will be implicitly ignored if pupilopd=None
"""Should calculations include the Science Instrument internal WFE?"""
self.options['jitter'] = 'gaussian'
self.options['jitter_sigma'] = constants.JWST_TYPICAL_LOS_JITTER_PER_AXIS
# class name to use for SI internal WFE, which can be overridden in subclasses
self._si_wfe_class = optics.WebbFieldDependentAberration
def _get_default_fov(self):
""" Return default FOV in arcseconds """
return 5 # default for all NIR instruments
def get_optical_system(self, fft_oversample=2, detector_oversample=None, fov_arcsec=2, fov_pixels=None, options=None):
# invoke superclass version of this
# then add a few display tweaks
optsys = SpaceTelescopeInstrument.get_optical_system(self,
fft_oversample=fft_oversample,
detector_oversample=detector_oversample,
fov_arcsec=fov_arcsec, fov_pixels=fov_pixels,
options=options)
# If the OTE model in the entrance pupil is a plain FITSOpticalElement, cast it to the linear model class
if not isinstance(optsys.planes[0], opds.OTE_Linear_Model_WSS):
lom_ote = opds.OTE_Linear_Model_WSS()
# FIXME seems like some code is missing here...? But in practice this code path
# never gets executed due to the _get_telescope_pupil_and_aberrations() function doing the right thing.
lom_ote
optsys.planes[0].display_annotate = utils.annotate_ote_pupil_coords
return optsys
def _get_aberrations(self):
""" return OpticalElement modeling wavefront aberrations for a given instrument,
including field dependence based on a lookup table of Zernike coefficients derived from
ISIM cryovac test data.
"""
if not self.include_si_wfe:
return None
optic = self._si_wfe_class(self)
return optic
def get_opd_file_full_path(self, opdfilename=None):
"""Return full path to the named OPD file.
The OPD may be:
- a local or absolute path,
- or relative implicitly within an SI directory, e.g. $WEBBPSF_PATH/NIRCam/OPD
- or relative implicitly within $WEBBPSF_PATH
This function handles filling in the implicit path in the latter cases.
"""
if opdfilename is None:
opdfilename = self.pupilopd
if os.path.exists(opdfilename):
return opdfilename
elif self.name in opdfilename:
return os.path.join(self._datapath, "OPD", opdfilename)
else:
return os.path.join(self._WebbPSF_basepath, opdfilename)
def _get_telescope_pupil_and_aberrations(self):
"""return OpticalElement modeling wavefront aberrations for the telescope.
This is nearly identical to the version of this function in SpaceTelescopeInstrument, differing only at the
very end. Here, we load the selected OPD file from disk into an instance of opds.OTE_Linear_Model_WSS if possible.
It falls back to a plain FITSOpticalElement for nonstandard sizes of input pupil, since the linear model is not
yet generalized to work on arbitrary sizes of pupil other than 1024 pixels.
See also get_aberrations for the SI aberrations.
"""
# ---- set pupil OPD
opd_index = None # default assumption: OPD file is not a datacube
if isinstance(self.pupilopd, str): # simple filename
opd_map = self.get_opd_file_full_path(self.pupilopd)
elif hasattr(self.pupilopd, '__getitem__') and isinstance(self.pupilopd[0], str):
# tuple with filename and slice, for a datacube
opd_map = self.get_opd_file_full_path(self.pupilopd[0])
opd_index = self.pupilopd[1]
elif isinstance(self.pupilopd, (fits.HDUList, poppy.OpticalElement)):
opd_map = self.pupilopd # not a path per se but this works correctly to pass it to poppy
elif self.pupilopd is None:
opd_map = None
else:
raise TypeError("Not sure what to do with a pupilopd of that type:" + str(type(self.pupilopd)))
# ---- set pupil intensity
if self.pupil is None:
raise RuntimeError("The pupil shape must be specified in the "
"instrument class or by setting self.pupil")
if isinstance(self.pupil, poppy.OpticalElement):
# supply to POPPY as-is
pupil_optic = self.pupil
else:
# wrap in an optic and supply to POPPY
if isinstance(self.pupil, str): # simple filename
if os.path.exists(self.pupil):
pupil_transmission = self.pupil
else:
pupil_transmission = os.path.join(
self._WebbPSF_basepath,
self.pupil
)
# Get npix from pupil_transmission
npix = int(pupil_transmission.split('npix')[-1].split('.')[0])
elif isinstance(self.pupil, fits.HDUList):
# POPPY can use self.pupil as-is
pupil_transmission = self.pupil
# Get npix from the shape of the data
npix = self.pupil[0].data.shape[0]
else:
raise TypeError("Not sure what to do with a pupil of "
"that type: {}".format(type(self.pupil)))
# ---- apply pupil intensity and OPD to the optical model
pupil_optic = opds.OTE_Linear_Model_WSS(
name='{} Entrance Pupil'.format(self.telescope),
transmission=pupil_transmission,
opd=opd_map,
opd_index=opd_index,
v2v3=self._tel_coords(), npix=npix,
include_nominal_field_dependence=self.include_ote_field_dependence
)
return pupil_optic
@SpaceTelescopeInstrument.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.
"""
# 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}')
if self.detector not in value:
raise ValueError(f'Aperture name {value} does not match currently selected detector {self.detector}. '
f'Change detector attribute first, then set desired aperture.')
# Only update if new value is different
if self._aperturename != value:
# First, check some info from current settings, wich 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(self._aperturename))
# 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(self._aperturename)
_log.debug(f"Pixelscale updated to {self.pixelscale} based on average X+Y SciScale at SIAF aperture {self._aperturename}")
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.
"""
return self._detector_geom_info.pix2angle(self.detector_position[0], self.detector_position[1])
def _xan_yan_coords(self):
""" Convert from detector pixel coordinates to the XAN, YAN coordinate system
which was used for much of ISIM optical testing. The origin of XAN, YAN is
centered at the master chief ray, which passes through the ISIM focal plane
between the NIRCam A3 and B4 detectors. The sign of YAN is flipped relative to V3.
"""
coords = self._tel_coords()
# XAN is the same as V2, therefore no change to first element
# YAN is opposite direction as V3, and offset by 468 arcseconds
coords[1] = -coords[1] - 468 * units.arcsec
return coords
def set_position_from_aperture_name(self, aperture_name):
""" Set the simulated center point of the array based on a named SIAF aperture.
This will adjust the detector and detector position attributes.
"""
try:
ap = self.siaf[aperture_name]
# setting the detector must happen -before- we set the position
detname = aperture_name.split('_')[0]
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:
# NIRSpec slit apertures need some separate handling, since they don't map directly to detector pixels
ref_in_tel = ap.V2Ref, ap.V3Ref
nrs_full_aperture = self.siaf[aperture_name.split('_')[0]+"_FULL"]
ref_in_sci = nrs_full_aperture.tel_to_sci(*ref_in_tel)
self.detector_position = ref_in_sci
_log.debug("From {} set det. pos. to {} {}".format(aperture_name, detname, self.detector_position))
except KeyError:
raise ValueError("Not a valid aperture name for {}: {}".format(self.name, aperture_name))
def _get_pixelscale_from_apername(self, apername):
"""Simple utility function to look up pixelscale from apername"""
ap = self.siaf[apername]