-
Notifications
You must be signed in to change notification settings - Fork 58
/
opds.py
3422 lines (2759 loc) · 144 KB
/
opds.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
###############################################################################
#
# JWST Linear Optical Models for Mirror Motions
#
# This file implements two independent linear optical models for
# how the JWST OTE wavefront will change in response to commanded mirror
# motions. There are two models for historical reasons; both give fairly
# consistent answers within the uncertainties set by the inherent limits of
# linearizing a nonlinear optical problem. One optical model is derived from
# the OTE influence functions in Elliott 2011, JWST-STScI-02356, "Sensitivity
# of wavefront error to segment motions in the JWST primary mirror."
# The second is derived from OTE influence functions consistent with those used
# in the JWST Wavefront Analysis Software by Ball Aerospace.
# All code and algorithms in this file are by Marshall Perrin
#
# This information is *not* considered sensitive under ITAR or EAR:
# See GSFC Export Control Checklist "JWST Performance and Calibration
# Technical Data Assessment (Jan. 2010)" which describes technical data
# on the JWST program that is not considered ITAR controlled, in particular
# items 1.1.2: "models and measures of wavefront error on both the telescope
# and in the instrument" and item 4.7: "the effects of the OTE's optical
# elements on a PSF."
# See also the STScI Mission Office June 2015 determination that the
# Elliott paper JWST-STScI-02356 is not ITAR sensitive, in part because
# its contents are derivable from the publicly available OTE optical
# prescription, as published for instance in Lightsey et al. 2012 Opt. Eng.
#
###############################################################################
import functools
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import scipy.special as sp
import scipy
import astropy.table
import astropy.io.fits as fits
import astropy.units as u
import logging
from collections import OrderedDict
import warnings
from packaging.version import Version
import copy
import poppy
import poppy.zernike as zernike
import pysiaf
from . import constants
from . import surs
from . import utils
import webbpsf
_log = logging.getLogger('webbpsf')
__doc__ = """
Linear Models for JWST segment motions
JWST_OTE_LOM_Elliot : OTE Linear optical model based on
the influence functions in Elliott 2011, JWST-STScI-02356.
JWST_OTE_LOM_WSS: OTE linear optical model based on the
influence functions used in the JWST WFSC Software System
These give very similar results, though not precisely consistent
because of small differences in the parameter ranges used when
deriving the linear approximations.
"""
__location__ = os.path.dirname(os.path.abspath(__file__)) + os.sep
################################################################################
class OPD(poppy.FITSOpticalElement):
""" Base class for JWST Optical Path Difference (OPD) files
Provides methods for analyzing and displaying OPDS.
If you want to manipulate one using the linear optical model, see the OPDbender class.
This class is implemented as a child of poppy.FITSOpticalElement; as such it
follows all the same behavior as that, including having all OPDs internally
in units of meters.
Key class attributes include
.OPD the OPD data (as a fits.PrimaryHDU)
.pupilmask a mask for the pupil (as a fits.ImageHDU)
.pixelscale Scale in meters/pixel
.header reference to the OPD data's header
"""
def __init__(self, name='unnamed OPD', opd=None, opd_index=0, transmission=None,
segment_mask_file=None, npix=1024, **kwargs):
"""
Parameters
----------
opd : string, path to FITS file.
FITS file to load an OPD from. The OPD must be specified in microns. This FITS file must be
compatible with the format expected by poppy.FITSOpticalElement.
transmission: str
FITS file for pupil mask, with throughput from 0-1. If not explicitly provided, will be inferred from
wherever is nonzero in the OPD file
npix: int
Number of pixels per side for the OPD. Can be 1024, 2048, 4096 with current data files.
segment_mask_file : string, or None
if None, will infer based on npix.
ext : int, optional
FITS extension to load OPD from
opd_index : int, optional
slice of a datacube to load OPD from, if the selected extension contains a datacube.
"""
self.npix = npix
if opd is None and transmission is None:
_log.debug('Neither a pupil mask nor OPD were specified. Using the default JWST pupil.')
transmission = os.path.join(utils.get_webbpsf_data_path(), f"jwst_pupil_RevW_npix{self.npix}.fits.gz")
super(OPD, self).__init__(name=name,
opd=opd, transmission=transmission,
opd_index=opd_index, transmission_index=0,
planetype=poppy.poppy_core.PlaneType.pupil, **kwargs)
# Check that the shape of the OPD that has been passed, matches the npix parameters
if self.opd is not None:
if not self.opd.shape == (self.npix, self.npix):
raise ValueError(f"npix value of {self.npix} does not match shape of OPD provided: {self.opd.shape}")
if self.opd_header is None:
self.opd_header = self.amplitude_header.copy()
self.segnames = np.asarray([a[0:2] for a in constants.SEGNAMES_WSS_ORDER])
if segment_mask_file is None:
segment_mask_file = f'JWpupil_segments_RevW_npix{self.npix}.fits.gz'
full_seg_mask_file = os.path.join(utils.get_webbpsf_data_path(), segment_mask_file)
if not os.path.exists(full_seg_mask_file):
# try without .gz
full_seg_mask_file = os.path.join(utils.get_webbpsf_data_path(), f"JWpupil_segments_RevW_npix{npix}.fits")
self._segment_masks = fits.getdata(full_seg_mask_file)
if not self._segment_masks.shape[0] == self.npix:
raise ValueError(f"The shape of the segment mask file {self._segment_masks.shape} does not match the shape expect: ({self.npix}, {self.npix})")
self._segment_masks_version = fits.getheader(full_seg_mask_file)['VERSION']
# Where are the centers of each segment? From OTE design geometry
self._seg_centers_m = {seg[0:2]: np.asarray(cen)
for seg, cen in constants.JWST_PRIMARY_SEGMENT_CENTERS}
# convert the center of each segment to pixels for the current array sampling:
self._seg_centers_pixels = {seg[0:2]: self.shape[0] / 2 + np.asarray(cen) / self.pixelscale.value
for seg, cen in constants.JWST_PRIMARY_SEGMENT_CENTERS}
# And what are the angles of the local control coordinate systems?
self._rotations = {}
self._control_xaxis_rotations = {}
self._control_xaxis_rot_base = {'A': 180, # Rotations of local control coord
'B': 0, # X axes, CCW relative to V2 axis
'C': 60} # for A,B,C1
for i in range(18):
seg = self.segnames[i]
self._rotations[seg] = (int(seg[1]) - 1) * 60 # seg_rotangles[i]
self._control_xaxis_rotations[seg] = (self._control_xaxis_rot_base[seg[0]] +
-1 * self._rotations[seg])
self._seg_tilt_angles = {'A': -4.7644, # Tilts of local normal vector
'B': 9.4210, # relative to the V1 axis, around
'C': -8.1919} # the local X axis, or Y axis for Cs
self.name = name
self.header = self.opd_header # convenience reference
def copy(self):
""" Make a copy of a wavefront object """
from copy import deepcopy
return deepcopy(self)
# We can add and subtract OPD objects using the following.
def __add__(self, x):
""" Operator add """
output = self.copy()
if isinstance(x, OPD):
output.opd += x.data
output.name += " + " + x.name
else:
output.opd += x
output.name += " + " + str(x)
return output
def __radd__(self, x):
return self.__add__(x)
def __sub__(self, x):
""" Operator subtract """
output = self.copy()
if isinstance(x, OPD):
output.opd -= x.opd
output.name += " - " + x.name
else:
output.opd -= x
output.name += " - " + str(x)
return output
def as_fits(self, include_pupil=True):
""" Return the OPD as a fits.HDUList object
Parameters
-----------
include_pupil : bool
Include the pupil mask as a FITS extension?
"""
output = fits.HDUList([fits.PrimaryHDU(self.opd, self.opd_header)])
output[0].header['EXTNAME'] = 'OPD'
output[0].header['BUNIT'] = 'meter' # Rescaled to meters in poppy_core
if include_pupil:
self.amplitude_header['EXTNAME'] = 'PUPIL'
output.append(fits.ImageHDU(self.amplitude, self.amplitude_header))
return output
def writeto(self, outname, overwrite=True, **kwargs):
""" Write OPD to a FITS file on disk """
self.as_fits(**kwargs).writeto(outname, overwrite=overwrite)
# ---- display and analysis
def powerspectrum(self, max_cycles=50, sampling=5, vmax=100, iterate=False):
"""
Compute the spatial power spectrum of an aberrated wavefront.
Produces nice plots on screen.
Returns an array [low, mid, high] giving the RMS spatial frequencies in the
different JWST-defined spatial frequency bins:
low: <=5 cycles/aperture
mid: 5 < cycles <= 30
high: 30 < cycles
"""
# import SFT
def tvcircle(radius=1, xcen=0, ycen=0, center=None, **kwargs):
"""
draw a circle on an image.
radius
xcen
ycen
center= tuple in (Y,X) order.
"""
if center is not None:
xcen = center[1]
ycen = center[0]
t = np.arange(0, np.pi * 2.0, 0.01)
t = t.reshape((len(t), 1))
x = radius * np.cos(t) + xcen
y = radius * np.sin(t) + ycen
plt.plot(x, y, **kwargs)
cmap = copy.copy(matplotlib.cm.get_cmap(poppy.conf.cmap_diverging))
cmap.set_bad('0.3')
plt.clf()
plt.subplot(231)
self.display(title='full wavefront', clear=False, colorbar=False, vmax=vmax)
ps_pixel_size = 1. / sampling # how many cycles per pixel
trans = SFT.SFT3(self.data, max_cycles * 2, max_cycles * 2 * sampling)
abstrans = np.abs(trans)
extent = [-max_cycles, max_cycles, -max_cycles, max_cycles]
plt.subplot(233)
plt.imshow(abstrans, extent=extent, origin='lower')
plt.title("Power Spectrum of the phase")
plt.ylabel("cycles/aperture")
tvcircle(radius=5, color='k', linewidth=1) # , ls='--')
tvcircle(radius=30, color='k', linewidth=1) # 2, ls='--')
plt.gca().set_xbound(-max_cycles, max_cycles)
plt.gca().set_ybound(-max_cycles, max_cycles)
y, x = np.indices(abstrans.shape)
y -= abstrans.shape[0] / 2.
x -= abstrans.shape[1] / 2.
r = np.sqrt(x ** 2 + y ** 2) * ps_pixel_size
mask = np.ones_like(self.data)
mask[np.where(self.amplitude == 0)] = np.nan
wgood = np.where(self.amplitude != 0)
components = []
for i, label in enumerate(['low', 'mid', 'high']):
plt.subplot(2, 3, i + 4)
if label == 'low':
condition = r <= 5
elif label == 'mid':
condition = (r > 5) & (r <= 30)
else:
condition = (r > 30)
filtered = trans * condition
inverse = SFT.SFT3(filtered, max_cycles * 2, self.opd.shape[0], inverse=True)
inverse = inverse[::-1, ::-1] # I thought SFT did this but apparently this is necessary to get the high freqs right...
plt.imshow(inverse.real * mask, vmin=(-vmax) / 1000., vmax=vmax / 1000,
cmap=cmap, origin='lower') # vmax is in nm, but WFE is in microns, so convert
plt.title(label + " spatial frequencies")
rms = (np.sqrt((inverse.real[wgood] ** 2).mean()) * 1000)
components.append(rms)
plt.xlabel("%.3f nm RMS WFE" % rms)
return np.asarray(components)
def display_opd(self, ax=None, labelsegs=True, vmax=150., colorbar=True, clear=False, title=None, unit='nm',
pupil_orientation='entrance_pupil',
cbpad=None, colorbar_orientation='vertical',
show_axes=False, show_rms=True, show_v2v3=False,
cmap=None):
""" Draw on screen the perturbed OPD
Parameters
-----------
ax : matplotlib.Axes
axes instance to display into.
labelsegs : bool
draw segment name labels on each segment? default True.
show_axes: bool
Draw local control axes per each segment
show_rms : bool
Annotate the RMS wavefront value
show_v2v3:
Draw the observatory V2V3 coordinate axes
pupil_orientation : string
either 'entrance_pupil' or 'exit_pupil', for which orientation we should display the OPD in.
clear : bool
Clear plot window before display? default true
unit : str
Unit for WFE. default is 'nm'
"""
if unit == 'nm':
scalefact = 1e9
elif unit == 'micron':
scalefact = 1e6
elif unit == 'm' or unit == 'meter':
scalefact = 1
else:
raise ValueError("unknown unit keyword")
# The actual OPD data is stored in "entrance pupil" orientation, looking in to JWST,
# by convention and for consistency with the JWST WSS.
# In the case of JWST it turns out the "exit pupil" orientation is in effect a vertical
# flip. This is because of (1) inversion in both X and Y axes from the physical entrance pupil
# at the primary to the real image of that pupil on the FSM which is the OTE exit pupil, plus
# (2) the change in viewing convention from "in front of OTE looking in" to "at the instruments
# looking outward". That results in a flip in the X axis. Thus overall we can just flip the
# Y axis to get to the exit pupil orientation. In this display function we can do that using the
# origin parameter to matplotlib.imshow. In the actual optical propagation we do that with a
# poppy coordinate transform instance.
if pupil_orientation=='entrance_pupil':
origin='lower'
elif pupil_orientation == 'exit_pupil':
origin = 'upper'
else:
raise ValueError("pupil_orientation must be 'entrance_pupil' or 'exit_pupil'.")
if clear:
if ax is not None:
raise RuntimeError("'clear=True' is incompatible with passing in an Axes instance.")
plt.clf()
if cmap is None:
cmap = copy.copy(matplotlib.cm.get_cmap(poppy.conf.cmap_diverging))
cmap.set_bad('0.3')
mask = np.ones_like(self.opd)
mask[np.where(self.amplitude == 0)] = np.nan
try:
pupilscale = self.opd_header['PUPLSCAL']
s = self.opd.shape
extent = [a * pupilscale for a in [-s[0] / 2, s[0] / 2, -s[1] / 2, s[1] / 2]]
except KeyError:
extent = None
if ax is None:
ax = plt.gca()
plot = ax.imshow(self.opd * mask * scalefact, vmin=-vmax, vmax=vmax, cmap=cmap, extent=extent, origin=origin)
_log.debug("Displaying OPD. Vmax is %f, data max is %f " % (vmax, self.opd.max()))
if title is None:
title = self.name
ax.set_title(title)
if show_rms:
ax.set_xlabel("RMS WFE = %.1f nm" % self.rms())
if show_v2v3:
utils.annotate_ote_pupil_coords(None, ax, orientation=pupil_orientation)
if labelsegs:
for seg in self.segnames[0:18]:
self.label_seg(seg, ax=ax, show_axes=show_axes, pupil_orientation=pupil_orientation)
if colorbar:
if cbpad is None:
cbpad = 0.05 if colorbar_orientation == 'vertical' else 0.15
cb = plt.colorbar(plot, ax=ax, pad=cbpad, orientation=colorbar_orientation)
cb.set_label("WFE [%s]" % unit)
else:
cb = None
plt.draw()
return ax, cb
def label_seg(self, segment, ax=None, show_axes=False, color='black', pupil_orientation='entrance_pupil'):
""" Annotate a plot with a text label for a particular segment """
cx, cy = self._seg_centers_m[segment]
# See note about pupil orientations in OPD.display_opd() for explanation
# of the Y axis signs here
if pupil_orientation == 'entrance_pupil':
ysign = 1
elif pupil_orientation == 'exit_pupil':
ysign = -1
else:
raise ValueError("pupil_orientation must be 'entrance_pupil' or 'exit_pupil'.")
offset = 0.2 if show_axes else 0
if ax is None:
ax = plt.gca()
label = ax.text(cx + offset, (cy + offset)*ysign, segment, color=color, horizontalalignment='center', verticalalignment='center')
if show_axes:
ax_arrow_len = .3
# if not ('C' in segment ):
if True:
for i, color, label in zip([0, 1, 2], ['green', 'blue', 'red'], ['x', 'y', 'z']):
vec = np.matrix([0, 0, 0]).transpose() # xyz order
vec[i] = 1
b = self._rot_matrix_local_to_global(segment) * vec
b = np.asarray(b).flatten() # Inelegant but it works
ax.arrow(cx, cy*ysign, ax_arrow_len * b[0], (ax_arrow_len * b[1])*ysign, color=color,
# width=ax,
head_width=.050, head_length=.080) # in units of mm
xoffset = 0.1 if i == 2 else 0
ax.text(cx + ax_arrow_len * b[0] * 1.5 + xoffset, (cy + ax_arrow_len * b[1] * 1.5)*ysign, label,
color=color, fontsize=8,
horizontalalignment='center', verticalalignment='center'
)
ax.get_figure().canvas.draw()
return label
def _rot_matrix_local_to_global(self, segname):
""" Rotation matrix from Local to Global coordinates
Inverse of _rot_matrix_global_to_local
"""
from . import geometry
tilt = self._seg_tilt_angles[segname[0]]
xaxis_rot = self._control_xaxis_rotations[segname]
if 'C' in segname: # Cs are tilted about Y, ABs tilted about X
return geometry.rot_matrix_z(xaxis_rot) * geometry.rot_matrix_y(tilt)
else:
return geometry.rot_matrix_z(xaxis_rot) * geometry.rot_matrix_x(tilt)
def _rot_matrix_global_to_local(self, segname):
""" Rotation matrix from Global to Local coordinates
"""
from . import geometry
tilt = self._seg_tilt_angles[segname[0]]
xaxis_rot = self._control_xaxis_rotations[segname]
if 'C' in segname: # Cs are tilted about Y, ABs tilted about X
return geometry.rot_matrix_y(-tilt) * geometry.rot_matrix_z(-xaxis_rot)
else:
return geometry.rot_matrix_x(-tilt) * geometry.rot_matrix_z(-xaxis_rot)
def zern_seg(self, segment, vmax=150, unit='nm'):
""" Show the Zernike terms applied to a given segment"""
# the actual wavefront is always in units of microns.
if unit == 'nm':
scalefact = 1000.
elif unit == 'micron':
scalefact = 1.0
else:
raise ValueError("unknown unit keyword")
nzerns = 11
title = "Zernike components for " + segment
zerns = np.zeros(nzerns)
if segment + "-decenter" in self.state.keys():
zerns += self._move(segment, type='decenter', vector=self.state[segment + "-decenter"], return_zernikes=True)
title += ", displacement=[%.2e, %.2e, %.2e] um" % tuple(self.state[segment + "-decenter"])
if segment + "-tilt" in self.state.keys():
zerns += self._move(segment, type='tilt', vector=self.state[segment + "-tilt"], return_zernikes=True)
title += ", tilts =[%.5f, %.5f, %.5f] urad" % tuple(self.state[segment + "-tilt"])
_log.info("Zerns: " + str(zerns))
fig = plt.gcf()
fig.clf()
npix = 200
hexap = zernike.hex_aperture(npix)
hexap[np.where(hexap == 0)] = np.nan
cmap = copy.copy(matplotlib.cm.get_cmap('jet'))
cmap.set_bad('0.5', alpha=0.0)
for j in np.arange(nzerns) + 1:
ax = fig.add_subplot(3, 4, j, frameon=False, xticks=[], yticks=[])
# n, m = zernike.noll_indices(j)
Z = zernike.zernike1(j, npix=npix)
ax.imshow(Z * zerns[j - 1] * hexap * scalefact, vmin=-1 * vmax, vmax=vmax, cmap=cmap, origin='lower')
ax.text(npix * 0.95, npix * 0.8, "$Z{:d}$".format(j), fontsize=20, horizontalalignment='right')
ax.text(npix * 0.95, npix * 0.1, "{:.2e}".format(zerns[j - 1]), fontsize=15, horizontalalignment='right')
fig.text(0.5, 0.95, title, horizontalalignment='center', fontsize=15)
fig.text(0.95, 0.15, 'Segment RMS WFE = {:.2f} {} '.format(self.rms(segment) * (scalefact / 1000), unit), fontsize=10,
horizontalalignment='right')
# fig.text(0.95, 0.05, 'OPDs scaled from %.2e - %.2e um' % (-vmax, vmax), horizontalalignment='right', fontsize=10)
fig.text(0.95, 0.05, 'OPDs scaled from {:.2f} - {:.2f} {}'.format(-vmax, vmax, unit), horizontalalignment='right', fontsize=10)
plt.draw()
def rms(self, segment=None):
""" Return RMS WFE in nanometers
Parameters
----------
segment : string
Segment name, to compute RMS for a single segment. Leave unspecified (None) to compute
RMS WFE for the entire aperture. Segments are identified by name: A1-A6, B1-B6, C1-C6
"""
if segment is None:
# RMS for whole aperture
wgood = np.where((self.amplitude != 0) & (np.isfinite(self.opd)))
else:
assert (segment in self.segnames)
iseg = np.where(self.segnames == segment)[0][0] + 1 # segment index from 1 - 18
wseg = np.where((self._segment_masks == iseg) & (np.isfinite(self.opd)))
wgood = wseg
return np.sqrt((self.opd[wgood] ** 2).mean()) * 1e9
def ptv(self, segment=None):
"""return peak to valley WFE in nanometers
Parameters
----------
segment : string
Segment name, to compute RMS for a single segment. Leave unspecified (None) to compute
for the entire aperture.
"""
if segment is None:
# RMS for whole aperture
wgood = np.where(self.amplitude != 0)
else:
assert (segment in self.segnames)
iseg = np.where(self.segnames == segment)[0][0] + 1 # segment index from 1 - 18
wseg = np.where(self._segment_masks == iseg)
wgood = wseg
peak = self.opd[wgood].max()
valley = self.opd[wgood].min()
return (peak - valley) * 1000
def estimated_Strehl(self, wavelength, verbose=True):
""" Compute an estimated Strehl given a wavelength in meters
Parameters
-----------
wavelength : float
in meters
verbose : bool
should I print out an informative message?
"""
# rms = sqrt((-1.0)*alog(strehl))*wavelength/(2*!pi)
rms = self.rms() * 1e-9
strehl = np.exp(-(rms * 2 * np.pi / wavelength) ** 2)
if verbose:
print("For wavelength = {0} meters, estimated Strehl = {1} for {2} nm rms WFE".format(wavelength, strehl, rms * 1e9))
return strehl
################################################################################
class OTE_Linear_Model_Elliott(OPD):
""" Perturb an existing JWST OTE wavefront OPD file, by applying changes in WFE
derived from Erin Elliott's linear optical model. See Elliott 20122, JWST-STScI-02356
Note that this model is strictly accurate only for a single NIRCam field
point, but we can apply it elsewhere and it should not be too far off for small motions.
"""
def __init__(self, opd=None, opd_index=0, transmission=None, rm_ptt=False, rm_piston=False, zero=False):
"""
Parameters
----------
opdfile : str or fits.HDUList
FITS file to load an OPD from. The OPD must be specified in microns.
ext : int, optional
FITS extension to load OPD from
opd_index : int, optional
slice of a datacube to load OPD from, if the selected extension contains a datacube.
pupilfile : str
FITS file for pupil mask, with throughput from 0-1. If not explicitly provided, will be inferred from
wherever is nonzero in the OPD file.
rm_ptt : bool
Remove piston, tip and tilt from all segments if set. Default is False.
rm_piston : bool
Remove piston only from all segments if set. Default is False
zero : bool
If set, reate an OPD which is precisely zero in all locations. Default is False.
"""
OPD.__init__(self, opd=opd, opd_index=opd_index, transmission=transmission)
self._sensitivities = astropy.table.Table.read(os.path.join(__location__, 'otelm', 'seg_sens.txt'), format='ascii', delimiter='\t')
self.state = {}
self.remove_piston_tip_tilt = rm_ptt
self.remove_piston_only = rm_piston
self._opd_original = self.opd.copy()
if zero:
self.zero()
# ---- overall state manipulation
def reset(self):
""" Reset an OPD to the state it was loaded from disk.
i.e. undo all segment moves.
"""
self.opd = self._opd_original.copy()
_log.info("Reset to unperturbed OPD")
def zero(self):
self.opd *= 0
self.state = {}
_log.info("Set OPD to zero WFE!")
def load_state(self, newstate):
_log.info("Loading new state.")
self.zero()
for k in newstate.keys():
segname, motion = k.split('-')
self._move(segname, motion, newstate[k])
def print_state_v1(self):
keys = self.state.keys()
keys.sort()
if keys is None:
print("No perturbations")
else:
print("state = {")
for k in keys:
print(" '%s' : %s," % (k, repr(self.state[k])))
print(" }")
def print_state(self, type='local'):
keys = self.state.keys()
if type == 'local': # control coords
print(
"Segment poses in local Control coordinates: "
"(microns for decenter & piston, microradians for tilts, milliradians for clocking):")
print(" \t %10s %10s %10s %10s %10s %10s" % ("X dec", "Y dec", "Z piston", "X tilt", "Y tilt", "Clocking"))
for segment in self.segnames[0:18]:
if segment + "-tilt" in keys:
tilts = self.state[segment + "-tilt"].tolist()
else:
tilts = [0, 0, 0]
if segment + "-decenter" in keys:
decenters = self.state[segment + "-decenter"].tolist()
else:
decenters = [0, 0, 0]
print("%2s\t %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f" % tuple([segment] + decenters + tilts))
elif type == 'Report':
raise NotImplementedError("Coord conversions need work")
print("Segment positions in Report coordinates: "
"(microns for decenter, microradians for alpha & beta, milliradians for gamma):")
print(" \t %10s %10s %10s %10s %10s %10s" % ("X dec", "Y dec", "Z dec", "alpha", "beta", "gamma"))
for segment in self.segnames:
if segment + "-tilt" in keys:
tilts = self.state[segment + "-tilt"].tolist()
else:
tilts = [0, 0, 0]
if segment + "-decenter" in keys:
decenters = self.state[segment + "-decenter"].tolist()
else:
decenters = [0, 0, 0]
print("%2s\t %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f" % tuple([segment] + decenters + tilts))
elif type == 'Matlab':
print("Segment positions in Matlab coordinates: (millimeters and degrees)")
print(" \t %12s %12s %12s %12s %12s %12s" % ("X dec", "Y dec", "Z dec", "alpha", "beta", "gamma"))
for segment in self.segnames:
if segment + "-tilt" in keys:
tilts = self.state[segment + "-tilt"].tolist()
else:
tilts = [0, 0, 0]
if segment + "-decenter" in keys:
decenters = self.state[segment + "-decenter"].tolist()
else:
decenters = [0, 0, 0]
decenters = (np.array([decenters[1], decenters[0], -1 * decenters[2]]) / 1000).tolist()
tilts[2] *= 1000 # convert millirad to microrad for consistency
tilts = (np.array([tilts[1], tilts[0], -1 * tilts[2]]) * 1e-6 * 180 / np.pi).tolist()
print("%2s\t %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e" % tuple([segment] + decenters + tilts))
# ---- segment manipulation via linear model
def _get_seg_sensitivities(self, segment='A1', type='decenter'):
assert (segment in self.segnames)
refseg = 2 if segment[0] == 'C' else 1
fields = ['%s_%s_%s%d' % (axis, type, segment[0], refseg) for axis in ['X', 'Y', 'Z']]
# divide by 1e6 since sensitivities in units of microns and we want meters
return np.asarray([self._sensitivities[f] for f in fields]) / 1e6
def _record(self, segment='A1', type='decenter', values=np.zeros(3)):
""" update the state structure with the current segment positions """
key = "%s-%s" % (segment, type)
if key in self.state.keys():
self.state[key] += values
else:
self.state[key] = values
def print_(self):
keys = self.state.keys()
keys.sort()
if keys is None:
print("No perturbations")
else:
for k in keys:
print("%s\t%s" % (k, str(self.state[k])))
def _apply_zernikes_to_seg(self, segment, zernike_coeffs, coordsys='local', debug=False):
""" Apply Zernike perturbations to a given segment
Parameters
----------
segment : string
name of segment, A1 through C6
zernike_coeffs : iterable of floats
Zernike coefficients, in units of microns
coordsys : string
Coordinate system to apply the Zernikes in, either "global" or "local".
Local is the BATC-defined "Control" coordinates for each segment.
"""
assert (segment in self.segnames)
iseg = np.where(self.segnames == segment)[0][0] + 1 # segment index from 1 - 18
wseg = np.where(self._segment_masks == iseg)
# determine the X and Y zernike coefficients for each segment
# determine the center of each segment, as the mean of the min and max X and Y values
# FIXME this should just use the BATC-provided center coordinates; see jwst_ote3d.py
Y, X = np.indices(self.opd.shape, dtype=float)
cx = np.mean([X[wseg].min(), X[wseg].max()])
cy = np.mean([Y[wseg].min(), Y[wseg].max()])
Y -= cy
X -= cx
# We have to compute the segment radius in pixels exactly here, to make sure we don't
# have any pixels that are > that radius, based on the discrete pixelization and the
# varying illuminated areas of some of the segments. FIXME this could be done better once
# the segment centers are update as noted above. I think.
# Hmm, and needs slight rounding up so add a ceil() call?
R = np.sqrt(X[wseg] ** 2 + Y[wseg] ** 2)
seg_radius = np.ceil(R.max())
_log.debug("Segment {} is centered at pixel loc ({:.1f}, {:.1f}) with radius {:.1f} pix".format(segment, cx, cy, seg_radius))
# Ah hell, this already does the rotations here??
ang = -self._rotations[segment] # This definitely has to be a positive sign,
# to get the right X, Y for locally-defined zernikes
# 2016-06: No that seems incorrect. Add negative sign;
# aha it depends on which kind of segment! This gets things
# working correctly for the A segments.
if 'B' in segment:
ang += 180 # As are opposite Bs
elif 'C' in segment:
ang += 60 # Cs are 60 deg back
ang = np.deg2rad(ang)
Xr = X * np.cos(ang) + Y * np.sin(ang)
Yr = X * -np.sin(ang) + Y * np.cos(ang)
# These are the BATC "Control" coordinates for each segment
Xrc = Xr[wseg]
Yrc = Yr[wseg]
# But Erin Elliott's linear optical model was computed on a slightly
# different coordinate system, which has a different sign convention,
# and which depends on segment type
if 'A' in segment:
Xrc_elliott = -Xrc
Yrc_elliott = Yrc
else:
Xrc_elliott = Xrc
Yrc_elliott = -Yrc
theta = np.arctan2(Yrc_elliott, Xrc_elliott)
Rw = np.sqrt(
Xrc_elliott ** 2 + Yrc_elliott ** 2) / seg_radius
# This normalizes the Zernike input radial coord to maximum 1 on the segment aperture.
if debug:
plt.subplot(121)
plt.imshow(Xr * (self._segment_masks == iseg), origin='lower')
plt.title("Local X_Control for " + segment)
plt.subplot(122)
plt.imshow(Yr * (self._segment_masks == iseg), origin='lower')
plt.title("Local Y_Control for" + segment)
plt.draw()
if self.remove_piston_tip_tilt:
zernike_coeffs[0:3] = 0
elif self.remove_piston_only:
zernike_coeffs[0] = 0
for i in range(len(zernike_coeffs)):
zern = zernike.zernike1(i + 1, rho=Rw, theta=theta) * zernike_coeffs[i]
self.opd[wseg] += zern
outtxt = "Zs=[" + ", ".join(['%.1e' % z for z in zernike_coeffs]) + "]"
_log.debug(" " + outtxt)
def _move(self, segment, type='tilt', vector=None, display=False, return_zernikes=False):
"""
Internal function that actually modifies the OPD by applying the linear optical model.
Don't use this directly, use tilt() or displace() instead.
The provided vector must be in **control** coordinates.
"""
# print "Segment rotation for %s is %f" % (segment, self._rotations[segment])
# ang = self._rotations[segment] * np.pi/180
# local_coordX = vector[0] * np.cos(ang) + vector[1] * np.sin(ang)
# local_coordY = vector[0] *-np.sin(ang) + vector[1] * np.cos(ang)
# local_vector = np.array([local_coordX, local_coordY, vector[2]])
local_vector = vector
if type == 'tilt':
local_vector[
2] /= 1000 # convert Z tilt to milliradians instead of microradians because that is what the sensitivity tables use
units = 'microradians for tip/tilt, milliradians for clocking'
else:
units = 'microns'
sensitivities = self._get_seg_sensitivities(segment, type=type)
zernike_coeffs_2d = sensitivities * local_vector[:, np.newaxis]
zernike_coeffs = zernike_coeffs_2d.sum(axis=0)
if return_zernikes:
return zernike_coeffs
else:
_log.info("Segment %s requested %s: %s in %s" % (segment, type, str(vector), units))
_log.debug(" local %s: %s" % (type, str(local_vector)))
self._record(segment, type, vector)
self._apply_zernikes_to_seg(segment, zernike_coeffs)
if display:
self.display()
def tilt(self, segment, tiltX=0.0, tiltY=0.0, tiltZ=0.0, unit='urad', display=False, coordsys='elliott'):
""" Tilt/rotate a segment some angle around X, Y, or Z.
Parameters
-----------
segment : str
Segment name, e.g. 'A1'
tiltX, tiltY, tiltZ : floats
Tilt angle, in microradians by default.
unit : str
Unit for displacements. Can be 'urad', 'radian', 'milliradian', 'arcsec', 'arcmin', 'milliarcsec'
display : bool
Display after moving?
coordsys : string
Name of coordinate system the input tilts are with respect to. Can be 'local', 'global', 'control', or 'elliott'
"""
if coordsys == 'control':
# flip signs to match Ball's control coordinates convention relative to Erin's LOM
# don't use *= to avoid editing the value in the calling function
if 'A' in segment:
tiltX = -1 * tiltX
elif 'B' in segment:
tiltY = -1 * tiltY
elif 'C' in segment:
tiltX = -1 * tiltX
tiltZ = -1 * tiltZ
tilts = np.array([tiltX, tiltY, tiltZ]).astype(float)
self.opd_header.add_history('Rotation: %s %s' % (str(tuple(tilts)), unit))
# new sensitivity matrices are in urad for alpha and beta, mrad for gamma.
# first just convert all to urad.
if unit.endswith('s'):
unit = unit[:-1]
unit = unit.lower()
if unit == 'urad':
pass
elif unit == 'milliarcsec':
tilts *= (1e6 * np.pi / (180. * 60 * 60 * 1000))
elif unit == 'arcsec':
tilts *= (1e6 * np.pi / (180. * 60 * 60))
elif unit == 'arcmin':
tilts *= (1e6 * np.pi / (180. * 60))
elif unit == 'radian' or unit == 'rad':
tilts *= 1e6
elif unit == 'milliradian' or unit == 'mrad':
tilts *= 1e3
else:
raise NotImplementedError('unknown unit')
self._move(segment, 'tilt', tilts, display=display)
def displace(self, segment, distX=0.0, distY=0.0, distZ=0.0, unit='micron', display=False, coordsys='elliott'):
""" Move a segment some distance in X, Y, and Z.
Parameters
-----------
segment : str
Segment name, e.g. 'A1'
distX, distY, distZ : floats
Displacement distance, in microns by default.
unit : str
Unit for displacements. Can be 'micron', 'millimeter','nanometer', 'mm', 'nm', 'um'
display : bool
Display after moving?
coordsys : string
Name of coordinate system the input tilts are with respect to. Can be 'control', 'global' or 'elliott'
"""
# Convert the provided coordinates into CONTROL coords
if coordsys == 'elliott':
pass # no conversion needed
elif coordsys == 'control':
# flip signs of displacements as needed, due to differing conventions and
# handedness. Depends on segment type.
if 'A' in segment:
# distY = -distY
distZ = -distZ
else:
# distX = -distX
distZ = -distZ
elif coordsys == 'global':
raise NotImplementedError('global not yet implemented')
else:
raise ValueError("invalid coordsys param")
vector = np.array([distX, distY, distZ])
self.opd_header.add_history('Displacement: %s %s' % (str(tuple(vector)), unit))
if unit.endswith('s'):
unit = unit[:-1]
unit = unit.lower()
if unit == 'micron' or unit == 'um':
pass
elif unit == 'millimeters' or unit == 'millimeter' or unit == 'mm':
vector *= 1000