-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathZernike_Module.py
3521 lines (2799 loc) · 176 KB
/
Zernike_Module.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
"""
First created on Mon Aug 13 10:01:03 2018
Versions
Oct 31, 2018; 0.1 -> 0.11 fixed FRD effect
Nov 1, 2018; 0.11 -> 0.12 added correct edges to the detector; fixed wrong behavior for misaligment
Nov 2, 2018; 0.12 -> 0.13 added lorentzian wings to the illumination of the pupil
Nov 3, 2018; 0.13 -> 0.13b fixed edges of detector when det_vert is not 1
Nov 12, 2018; 0.13b -> 0.13c changed parameter describing hexagonal effect "f" from 0.1 to 0.2
Nov 12, 2018; 0.13c -> 0.14 changed illumination description modifying entrance -> exit pupil illumination
Nov 29, 2018; 0.14 -> 0.14b added fixed scattering slope, deduced from large image in focus
Dec 16, 2018; 0.14b -> 0.14c allparameters_proposal_err from list to array
Dec 18, 2018; 0.14c -> 0.14d strutFrac upper limit to 0.13 in create_parInit
Dec 23, 2018; 0.14d -> 0.15 refactoring so that x_ilum and y_ilum is one
Dec 26, 2018; 0.15 -> 0.15b when in focus, create exactly 10x oversampling
Dec 31, 2018; 0.15b -> 0.16 major rewrite of downsampling algorithm
Jan 8, 2019; 0.16 -> 0.17 added support for zmax=22
Jan 14, 2019; 0.17 -> 0.18 fixed bug with dowsampling algorithm - I was just taking central values
Jan 15, 2019; 0.18 -> 0.19 added simple algorithm to interpolate between 1/10 pixels in the best position
Feb 15, 2019; 0.19 -> 0.20 updated analysis for the new data
Feb 21, 2019; 0.20 -> 0.20b test parameter for showing globalparamers outside their limits
Feb 22, 2019; 0.20 -> 0.21 added support for Zernike higher than 22
Feb 22, 2019; 0.21 -> 0.21b added support for return image along side likelihood
Apr 17, 2019; 0.21b -> 0.21c changed defintion of residuals from (model-data) to (data-model)
Jun 4, 2019; 0.21c -> 0.21d slight cleaning of the code, no functional changes
Jun 26, 2019; 0.21d -> 0.21e included variable ``dataset'', which denots which data we are using in the analysis
Jul 29, 2019; 0.21e -> 0.21f changed the spread of paramters when drawing initial solutions, based on data
Sep 11, 2019; 0.21f -> 0.21g globalparameters_flat_6<1 to globalparameters_flat_6<=1
Oct 10, 2019: 0.21g -> 0.21h scattered_light_kernel saving option
Oct 31, 2019: 0.21h -> 0.22 (re)introduced small amount of apodization (PIPE2D-463)
Oct 31, 2019: 0.22 -> 0.22b introduced verbosity
Nov 07, 2019: 0.22b -> 0.22c nan values can pass through find_single_realization_min_cut
Nov 08, 2019: 0.22c -> 0.22d changes to resizing and centering
Nov 13, 2019: 0.22d -> 0.23 major changes to centering - chief ray in the center of oversampled image
Nov 15, 2019: 0.23 -> 0.24 change likelihood definition
Dec 16, 2019: 0.24 -> 0.24a added iluminaton with z4,z11,z22=0
Jan 14, 2020: 0.24a -> 0.24b added verbosity in find_single_realization_min_cut function
Jan 31, 2020: 0.24b -> 0.25 added support for data contaning spots from two wavelengths
Feb 11, 2020: 0.25 -> 0.26 proper bilinear interpolation of the spots
@author: Neven Caplar
@contact: [email protected]
@web: www.ncaplar.com
"""
########################################
#standard library imports
from __future__ import absolute_import, division, print_function
import os
import time
import sys
import math
import socket
os.environ["MKL_NUM_THREADS"] = "1"
os.environ["NUMEXPR_NUM_THREADS"] = "1"
os.environ["OMP_NUM_THREADS"] = "1"
import numpy as np
np.set_printoptions(suppress=True)
np.seterr(divide='ignore', invalid='ignore')
#print(np.__config__)
from multiprocessing import current_process
from functools import lru_cache
#import pyfftw
#import pandas as pd
########################################
# Related third party imports
# none at the moment
########################################
# Local application/library specific imports
# galsim
import galsim
galsim.GSParams.maximum_fft_size=12000
# astropy
import astropy
import astropy.convolution
from astropy.convolution import Gaussian2DKernel
# scipy and skimage
import scipy.misc
import skimage.transform
import scipy.optimize as optimize
from scipy.ndimage.filters import gaussian_filter
#lmfit
import lmfit
#matplotlib
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
# needed for resizing routines
from typing import Tuple, Iterable
########################################
__all__ = ['PupilFactory', 'Pupil','ZernikeFitter_PFS','LN_PFS_single','LNP_PFS',\
'find_centroid_of_flux','create_parInit',\
'Zernike_Analysis','PFSPupilFactory','custom_fftconvolve','stepK','maxK',\
'sky_scale','sky_size','create_x','remove_pupil_parameters_from_all_parameters',\
'resize','_interval_overlap']
__version__ = "0.26"
############################################################
# name your directory where you want to have files!
PSF_DIRECTORY='/Users/nevencaplar/Documents/PFS/'
############################################################
TESTING_FOLDER=PSF_DIRECTORY+'Testing/'
TESTING_PUPIL_IMAGES_FOLDER=TESTING_FOLDER+'Pupil_Images/'
TESTING_WAVEFRONT_IMAGES_FOLDER=TESTING_FOLDER+'Wavefront_Images/'
TESTING_FINAL_IMAGES_FOLDER=TESTING_FOLDER+'Final_Images/'
class Pupil(object):
"""!Pupil obscuration function.
"""
def __init__(self, illuminated, size, scale):
"""!Construct a Pupil
@param[in] illuminated 2D numpy array indicating which parts of
the pupil plane are illuminated.
@param[in] size Size of pupil plane array in meters. Note
that this may be larger than the actual
diameter of the illuminated pupil to
accommodate zero-padding.
@param[in] scale Sampling interval of pupil plane array in
meters.
"""
self.illuminated = illuminated
self.size = size
self.scale = scale
class PupilFactory(object):
"""!Pupil obscuration function factory for use with Fourier optics.
"""
def __init__(self, pupilSize, npix,input_angle,hscFrac,strutFrac,slitFrac,slitFrac_dy,x_fiber,y_fiber,effective_ilum_radius,frd_sigma,frd_lorentz_factor,det_vert,verbosity=None):
"""!Construct a PupilFactory.
@params[in] others
@param[in] npix Constructed Pupils will be npix x npix.
"""
self.verbosity=verbosity
if self.verbosity==1:
print('Entering PupilFactory class')
self.pupilSize = pupilSize
self.npix = npix
self.input_angle=input_angle
self.hscFrac=hscFrac
self.strutFrac=strutFrac
self.pupilScale = pupilSize/npix
self.slitFrac=slitFrac
self.slitFrac_dy=slitFrac_dy
self.effective_ilum_radius=effective_ilum_radius
self.frd_sigma=frd_sigma
self.frd_lorentz_factor=frd_lorentz_factor
self.det_vert=det_vert
u = (np.arange(npix, dtype=np.float64) - (npix - 1)/2) * self.pupilScale
self.u, self.v = np.meshgrid(u, u)
@staticmethod
def _pointLineDistance(p0, p1, p2):
"""Compute the right-angle distance between the points given by `p0`
and the line that passes through `p1` and `p2`.
@param[in] p0 2-tuple of numpy arrays (x,y coords)
@param[in] p1 2-tuple of scalars (x,y coords)
@param[in] p2 2-tuple of scalars (x,y coords)
@returns numpy array of distances; shape congruent to p0[0]
"""
x0, y0 = p0
x1, y1 = p1
x2, y2 = p2
dy21 = y2 - y1
dx21 = x2 - x1
return np.abs(dy21*x0 - dx21*y0 + x2*y1 - y2*x1)/np.hypot(dy21, dx21)
def _fullPupil(self):
"""Make a fully-illuminated Pupil.
@returns Pupil
"""
illuminated = np.ones(self.u.shape, dtype=np.float32)
return Pupil(illuminated, self.pupilSize, self.pupilScale)
def _cutCircleInterior(self, pupil, p0, r):
"""Cut out the interior of a circular region from a Pupil.
@param[in,out] pupil Pupil to modify in place
@param[in] p0 2-tuple indicating region center
@param[in] r Circular region radius
"""
r2 = (self.u - p0[0])**2 + (self.v - p0[1])**2
pupil.illuminated[r2 < r**2] = False
def _cutCircleExterior(self, pupil, p0, r):
"""Cut out the exterior of a circular region from a Pupil.
@param[in,out] pupil Pupil to modify in place
@param[in] p0 2-tuple indicating region center
@param[in] r Circular region radius
"""
r2 = (self.u - p0[0])**2 + (self.v - p0[1])**2
pupil.illuminated[r2 > r**2] = False
def _cutEllipseExterior(self, pupil, p0, r, b, thetarot):
"""Cut out the exterior of a circular region from a Pupil.
@param[in,out] pupil Pupil to modify in place
@param[in] p0 2-tuple indicating region center
@param[in] r Ellipse region radius = major axis
@param[in] b Ellipse region radius = minor axis
@param[in] thetarot Ellipse region rotation
"""
r2 = (self.u - p0[0])**2 + (self.v - p0[1])**2
theta=np.arctan(self.u/self.v)+thetarot
pupil.illuminated[r2 > r**2*b**2/(b**2*(np.cos(theta))**2+r**2*(np.sin(theta))**2)] = False
def _cutSquare(self,pupil, p0, r,angle,det_vert):
"""Cut out the interior of a circular region from a Pupil.
@param[in,out] pupil Pupil to modify in place
@param[in] p0 2-tuple indicating region center
@param[in] r half lenght of the length of square side
@param[in] angle angle that the camera is rotated
"""
pupil_illuminated_only1=np.ones_like(pupil.illuminated)
time_start_single_square=time.time()
###########################################################
# Central square
if det_vert is None:
det_vert=1
x21 = -r/2*det_vert*1
x22 = +r/2*det_vert*1
y21 = -r/2*1
y22 = +r/2*1
angleRad = angle
pupil_illuminated_only1[np.logical_and(((self.u-p0[0])*np.cos(-angle)+(self.v-p0[1])*np.sin(-angleRad)<x22) & \
((self.u-p0[0])*np.cos(-angleRad)+(self.v-p0[1])*np.sin(-angleRad)>x21),\
((self.v-p0[1])*np.cos(-angleRad)-(self.u-p0[0])*np.sin(-angleRad)<y22) & \
((self.v-p0[1])*np.cos(-angleRad)-(self.u-p0[0])*np.sin(-angleRad)>y21))] = False
f=0.2
###########################################################
# Lower right corner
x21 = -r/2
x22 = +r/2
y21 = -r/2*det_vert
y22 = +r/2*det_vert
angleRad21=-np.pi/4
triangle21=[[p0[0]+x22,p0[1]+y21],[p0[0]+x22,p0[1]+y21-y21*f],[p0[0]+x22-x22*f,p0[1]+y21]]
p21=triangle21[0]
y22=(triangle21[1][1]-triangle21[0][1])/np.sqrt(2)
y21=0
x21=(triangle21[2][0]-triangle21[0][0])/np.sqrt(2)
x22=-(triangle21[2][0]-triangle21[0][0])/np.sqrt(2)
pupil_illuminated_only1[np.logical_and(((self.u-p21[0])*np.cos(-angleRad21)+(self.v-p21[1])*np.sin(-angleRad21)<x22) & \
((self.u-p21[0])*np.cos(-angleRad21)+(self.v-p21[1])*np.sin(-angleRad21)>x21),\
((self.v-p21[1])*np.cos(-angleRad21)-(self.u-p21[0])*np.sin(-angleRad21)<y22) & \
((self.v-p21[1])*np.cos(-angleRad21)-(self.u-p21[0])*np.sin(-angleRad21)>y21)) ] = True
###########################################################
# Upper left corner
x21 = -r/2*1
x22 = +r/2*1
y21 = -r/2*det_vert
y22 = +r/2*det_vert
angleRad12=-np.pi/4
triangle12=[[p0[0]+x21,p0[1]+y22],[p0[0]+x21,p0[1]+y22-y22*f],[p0[0]+x21-x21*f,p0[1]+y22]]
p21=triangle12[0]
y22=0
y21=(triangle12[1][1]-triangle12[0][1])/np.sqrt(2)
x21=-(triangle12[2][0]-triangle12[0][0])/np.sqrt(2)
x22=+(triangle12[2][0]-triangle12[0][0])/np.sqrt(2)
pupil_illuminated_only1[np.logical_and(((self.u-p21[0])*np.cos(-angleRad12)+(self.v-p21[1])*np.sin(-angleRad12)<x22) & \
((self.u-p21[0])*np.cos(-angleRad12)+(self.v-p21[1])*np.sin(-angleRad12)>x21),\
((self.v-p21[1])*np.cos(-angleRad12)-(self.u-p21[0])*np.sin(-angleRad12)<y22) & \
((self.v-p21[1])*np.cos(-angleRad12)-(self.u-p21[0])*np.sin(-angleRad12)>y21)) ] = True
###########################################################
# Upper right corner
x21 = -r/2*1
x22 = +r/2*1
y21 = -r/2*det_vert
y22 = +r/2*det_vert
angleRad12=np.pi/4
triangle22=[[p0[0]+x22,p0[1]+y22],[p0[0]+x22,p0[1]+y22-y22*f],[p0[0]+x22-x22*f,p0[1]+y22]]
p21=triangle22[0]
y22=-0
y21=+(triangle22[1][1]-triangle22[0][1])/np.sqrt(2)
x21=+(triangle22[2][0]-triangle22[0][0])/np.sqrt(2)
x22=-(triangle22[2][0]-triangle22[0][0])/np.sqrt(2)
pupil_illuminated_only1[np.logical_and(((self.u-p21[0])*np.cos(-angleRad12)+(self.v-p21[1])*np.sin(-angleRad12)<x22) & \
((self.u-p21[0])*np.cos(-angleRad12)+(self.v-p21[1])*np.sin(-angleRad12)>x21),\
((self.v-p21[1])*np.cos(-angleRad12)-(self.u-p21[0])*np.sin(-angleRad12)<y22) & \
((self.v-p21[1])*np.cos(-angleRad12)-(self.u-p21[0])*np.sin(-angleRad12)>y21)) ] = True
###########################################################
# Lower right corner
x21 = -r/2*1
x22 = +r/2*1
y21 = -r/2*det_vert
y22 = +r/2*det_vert
angleRad12=np.pi/4
triangle11=[[p0[0]+x21,p0[1]+y21],[p0[0]+x21,p0[1]+y21-y21*f],[p0[0]+x21-x21*f,p0[1]+y21]]
p21=triangle11[0]
y22=-(triangle22[1][1]-triangle22[0][1])/np.sqrt(2)
y21=0
x21=+(triangle22[2][0]-triangle22[0][0])/np.sqrt(2)
x22=-(triangle22[2][0]-triangle22[0][0])/np.sqrt(2)
pupil_illuminated_only1[np.logical_and(((self.u-p21[0])*np.cos(-angleRad12)+(self.v-p21[1])*np.sin(-angleRad12)<x22) & \
((self.u-p21[0])*np.cos(-angleRad12)+(self.v-p21[1])*np.sin(-angleRad12)>x21),\
((self.v-p21[1])*np.cos(-angleRad12)-(self.u-p21[0])*np.sin(-angleRad12)<y22) & \
((self.v-p21[1])*np.cos(-angleRad12)-(self.u-p21[0])*np.sin(-angleRad12)>y21)) ] = True
pupil.illuminated=pupil.illuminated*pupil_illuminated_only1
time_end_single_square=time.time()
if self.verbosity==1:
print('Time for cutting out the square is '+str(time_end_single_square-time_start_single_square))
def _cutRay(self, pupil, p0, angle, thickness,angleunit=None):
"""Cut out a ray from a Pupil.
@param[in,out] pupil Pupil to modify in place
@param[in] p0 2-tuple indicating ray starting point
@param[in] angle Ray angle measured CCW from +x.
@param[in] thickness Thickness of cutout
"""
if angleunit is None:
angleRad = angle.asRadians()
else:
angleRad = angle
# the 1 is arbitrary, just need something to define another point on
# the line
p1 = (p0[0] + 1, p0[1] + np.tan(angleRad))
d = PupilFactory._pointLineDistance((self.u, self.v), p0, p1)
pupil.illuminated[(d < 0.5*thickness) &
((self.u - p0[0])*np.cos(angleRad) +
(self.v - p0[1])*np.sin(angleRad) >= 0)] = False
def _addRay(self, pupil, p0, angle, thickness,angleunit=None):
"""Add a ray from a Pupil.
@param[in,out] pupil Pupil to modify in place
@param[in] p0 2-tuple indicating ray starting point
@param[in] angle Ray angle measured CCW from +x.
@param[in] thickness Thickness of cutout
"""
if angleunit is None:
angleRad = angle.asRadians()
else:
angleRad = angle
# the 1 is arbitrary, just need something to define another point on
# the line
p1 = (p0[0] + 1, p0[1] + np.tan(angleRad))
d = PupilFactory._pointLineDistance((self.u, self.v), p0, p1)
pupil.illuminated[(d < 0.5*thickness) &
((self.u - p0[0])*np.cos(angleRad) +
(self.v - p0[1])*np.sin(angleRad) >= 0)] = True
class PFSPupilFactory(PupilFactory):
"""!Pupil obscuration function factory for PFS
"""
def __init__(self, pupilSize, npix,input_angle,hscFrac,strutFrac,slitFrac,slitFrac_dy,x_fiber,y_fiber,effective_ilum_radius,frd_sigma,frd_lorentz_factor,det_vert,slitHolder_frac_dx,verbosity=None):
"""!Construct a PupilFactory.
@param[in] visitInfo VisitInfo object for a particular exposure.
@param[in] pupilSize Size in meters of constructed Pupils.
@param[in] npix Constructed Pupils will be npix x npix.
"""
self.verbosity=verbosity
if self.verbosity==1:
print('Entering PFSPupilFactory class')
PupilFactory.__init__(self, pupilSize,npix,input_angle,hscFrac,strutFrac,slitFrac,slitFrac_dy,x_fiber,y_fiber,effective_ilum_radius,frd_sigma,frd_lorentz_factor,det_vert,verbosity=self.verbosity)
self.x_fiber=x_fiber
self.y_fiber=y_fiber
self.slitHolder_frac_dx=slitHolder_frac_dx
self._spiderStartPos=[np.array([ 0., 0.]), np.array([ 0., 0.]), np.array([ 0., 0.])]
self._spiderAngles=[1.57-1.57,3.66519-1.57,5.75959-1.57]
self.effective_ilum_radius=effective_ilum_radius
def _horizonRotAngle(self,resultunit=None):
"""!Compute rotation angle of camera with respect to horizontal
coordinates from self.visitInfo.
@returns horizon rotation angle.
"""
if resultunit is None:
print('error')
#parAng = Angle(self.input_angle)
#return parAng.wrap()
else:
return 0
def getPupil(self, point):
"""!Calculate a Pupil at a given point in the focal plane.
@param point Point2D indicating focal plane coordinates.
@returns Pupil
"""
if self.verbosity==1:
print('Entering getPupil')
# called subaruRadius as it was taken from the code fitting pupil for HSC on Subaru
subaruRadius = (self.pupilSize/2)*1
hscFrac = self.hscFrac # linear fraction
hscRadius = hscFrac * subaruRadius
slitFrac = self.slitFrac # linear fraction
subaruSlit = slitFrac*subaruRadius
strutFrac = self.strutFrac # linear fraction
subaruStrutThick = strutFrac*subaruRadius
# y-position of the slit
slitFrac_dy = self.slitFrac_dy
# relic from the HSC code
# See DM-8589 for more detailed description of following parameters
# d(lensCenter)/d(theta) in meters per degree
# lensRate = 0.0276 * 3600 / 128.9 * subaruRadius
# d(cameraCenter)/d(theta) in meters per degree
hscRate = 2.62 / 1000 * subaruRadius
hscPlateScale = 380
thetaX = point[0] * hscPlateScale
thetaY = point[1] * hscPlateScale
pupil = self._fullPupil()
camX = thetaX * hscRate
camY = thetaY * hscRate
# creating FRD effects
single_element=np.linspace(-1,1,len(pupil.illuminated), endpoint=True)
u_manual=np.tile(single_element,(len(single_element),1))
v_manual=np.transpose(u_manual)
center_distance=np.sqrt((u_manual-self.x_fiber*hscRate*hscPlateScale*12)**2+(v_manual-self.y_fiber*hscRate*hscPlateScale*12)**2)
frd_sigma=self.frd_sigma
sigma=2*frd_sigma
pupil_frd=(1/2*(scipy.special.erf((-center_distance+self.effective_ilum_radius)/sigma)+scipy.special.erf((center_distance+self.effective_ilum_radius)/sigma)))
pupil_lorentz=(np.arctan(2*(self.effective_ilum_radius-center_distance)/(4*sigma))+np.arctan(2*(self.effective_ilum_radius+center_distance)/(4*sigma)))/(2*np.arctan((2*self.effective_ilum_radius)/(4*sigma)))
pupil.illuminated= (pupil_frd+1*self.frd_lorentz_factor*pupil_lorentz)/(1+self.frd_lorentz_factor)
# Cout out the acceptance angle of the camera
self._cutCircleExterior(pupil, (0.0, 0.0), subaruRadius)
# Cut out detector shadow
self._cutSquare(pupil, (camX, camY), hscRadius,self.input_angle,self.det_vert)
#No vignetting of this kind for the spectroscopic camera
#self._cutCircleExterior(pupil, (lensX, lensY), lensRadius)
# Cut out spider shadow
for pos, angle in zip(self._spiderStartPos, self._spiderAngles):
x = pos[0] + camX
y = pos[1] + camY
self._cutRay(pupil, (x, y), angle, subaruStrutThick,'rad')
# cut out slit shadow
self._cutRay(pupil, (2,slitFrac_dy/18),-np.pi,subaruSlit,'rad')
# cut out slit holder shadow
#also subaruSlit/3 not fitted, just put roughly correct number
self._cutRay(pupil, (self.slitHolder_frac_dx/18,1),-np.pi/2,subaruSlit/3,'rad')
if self.verbosity==1:
print('Finished with getPupil')
return pupil
class ZernikeFitter_PFS(object):
"""!Class to create donut images in PFS
The model is constructed using FFT, and consists of the convolution of
an OpticalPSF and an input fiber image. The OpticalPSF part includes the
specification of an arbitrary number of zernike wavefront aberrations.
This code uses lmfit to initalize the parameters.
"""
def __init__(self,image=None,image_var=None,pixelScale=None,wavelength=None,
jacobian=None,diam_sic=None,npix=None,pupilExplicit=None,
wf_full_Image=None,radiometricEffectArray_Image=None,
ilum_Image=None,dithering=None,save=None,
pupil_parameters=None,use_pupil_parameters=None,use_optPSF=None,
zmaxInit=None,extraZernike=None,simulation_00=None,verbosity=None,
double_sources=None,double_sources_positions_ratios=None,*args):
"""
@param image image to analyze
@param image_var variance image
@param pixelScale pixel scale in arcseconds
@param wavelength wavelenght
@param jacobian An optional 2x2 Jacobian distortion matrix to apply
to the forward model. Note that this is relative to
the pixelScale above. Default is the identity matrix.
@param diam_sic diameter of the exit pupil
@param npix number of pixels describing the pupil
@param pupilExplicit if you want to pass explicit image of the exit pupil instead of creating one within the class
"""
# if you do not pass the image that you wish to compare, the model will default to creating 41x41 image
if image is None:
image=np.ones((41,41))
self.image = image
else:
self.image = image
# if you do not pass the variance image that you wish to compare, the default variance image has value of '1' everywhere
if image_var is None:
image_var=np.ones((41,41))
self.image_var=image_var
else:
self.image_var = image_var
#flux = number of counts in the image
flux = float(np.sum(image))
self.flux=flux
if jacobian is None:
jacobian = np.eye(2, dtype=np.float64)
else:
self.jacobian = jacobian
# if you do not pass the value for wavelength it will default to 794 nm, which is roughly in the middle of the red detector
if wavelength is None:
wavelength=794
self.wavelength=wavelength
else:
self.wavelength=wavelength
# This is size of the pixel in arcsec for PFS red arm in focus
# calculated with http://www.wilmslowastro.com/software/formulae.htm
# pixel size in microns/focal length in mm x 206.3
# pixel size =15 microns, focal length = 149.2 mm (138 aperature x1.1 f number)
if pixelScale is None:
pixelScale=20.76
self.pixelScale=pixelScale
else:
self.pixelScale=pixelScale
#print('self.pixelScale: '+str(self.pixelScale))
# Exit pupil size in focus, default is 139.5237e-3 meters (taken from Zemax)
if diam_sic is None:
diam_sic=139.5327e-3
self.diam_sic=diam_sic
else:
self.diam_sic=diam_sic
# when creating pupils it will have size of npix pixels
if npix is None:
npix=1024
self.npix=npix
else:
self.npix=npix
# puilExplicit can be used to pass explicitly the image of the pupil instead of creating it from the supplied parameters
if pupilExplicit is None:
pupilExplicit==False
self.pupilExplicit=pupilExplicit
else:
self.pupilExplicit=pupilExplicit
# radiometricEffectArray_Image can be used to pass explicitly the illumination of the exit pupil instead of creating it from the supplied parameters
if radiometricEffectArray_Image is None:
radiometricEffectArray_Image==False
self.radiometricEffectArray_Image=radiometricEffectArray_Image
else:
self.radiometricEffectArray_Image=radiometricEffectArray_Image
# when creating pupils it will have size of npix pixels
if dithering is None:
dithering=1
self.dithering=dithering
self.pixelScale=self.pixelScale/dithering
else:
self.dithering=dithering
self.pixelScale=self.pixelScale/dithering
if save in (None,0):
save=None
self.save=save
else:
save=1
self.save=save
if pupil_parameters is None:
self.pupil_parameters=pupil_parameters
else:
self.pupil_parameters=pupil_parameters
if use_pupil_parameters is None:
self.use_pupil_parameters=use_pupil_parameters
else:
self.use_pupil_parameters=use_pupil_parameters
self.args = args
if use_optPSF is None:
self.use_optPSF=use_optPSF
else:
self.use_optPSF=use_optPSF
self.zmax=zmaxInit
self.simulation_00=simulation_00
self.extraZernike=extraZernike
self.verbosity=verbosity
self.double_sources=double_sources
self.double_sources_positions_ratios=double_sources_positions_ratios
def initParams(self,z4Init=None, dxInit=None,dyInit=None,hscFracInit=None,strutFracInit=None,
focalPlanePositionInit=None,fiber_rInit=None,
slitFracInit=None,slitFrac_dy_Init=None,apodizationInit=None,radiometricEffectInit=None,
trace_valueInit=None,serial_trace_valueInit=None,pixel_effectInit=None,backgroundInit=None,
x_ilumInit=None,y_ilumInit=None,radiometricExponentInit=None,
x_fiberInit=None,y_fiberInit=None,effective_ilum_radiusInit=None,
grating_linesInit=None,scattering_radiusInit=None,scattering_slopeInit=None,scattering_amplitudeInit=None,fluxInit=None,frd_sigmaInit=None,frd_lorentz_factorInit=None,
det_vertInit=None,slitHolder_frac_dxInit=None):
"""Initialize lmfit Parameters object.
@param zmax Total number of Zernike aberrations used
@param z4Init Initial Z4 aberration value in waves (that is 2*np.pi*wavelengths)
# pupil parameters
@param hscFracInit Value determining how much of the exit pupil obscured by the central obscuration(detector)
@param strutFracInit Value determining how much of the exit pupil is obscured by a single strut
@param focalPlanePositionInit 2-tuple for position of the central obscuration(detector) in the focal plane
@param slitFracInit Value determining how much of the exit pupil is obscured by slit
@param slitFrac_dy_Init Value determining what is the vertical position of the slit in the exit pupil
#non-uniform illumination
@param radiometricEffectInit parameter describing non-uniform illumination of the pupil (1-params['radiometricEffect']**2*r**2)**(params['radiometricExponent'])
@param radiometricExponentInit parameter describing non-uniform illumination of the pupil (1-params['radiometricEffect']**2*r**2)**(params['radiometricExponent'])
@param x_ilumInit x-position of the center of illumination of the exit pupil
@param y_ilumInit y-position of the center of illumination of the exit pupil
# illumination due to fiber, parameters
@param x_fiberInit position of the fiber misaligment in the x direction
@param y_fiberInit position of the fiber misaligment in the y direction
@param effective_ilum_radiusInit fraction of the maximal radius of the illumination of the exit pupil
@param frd_sigma sigma of Gaussian convolving only outer edge, mimicking FRD
@param frd_lorentz_factor strength of the lorentzian factor describing wings of the pupil illumination
# further pupil parameters
@param det_vert multiplicative factor determining vertical size of the detector obscuration
@param slitHolder_frac_dx dx position of slit holder
# convolving (postprocessing) parameters
@param grating_lines number of effective lines in the grating
@param scattering_slopeInit slope of scattering
@param scattering_amplitudeInit amplitude of scattering compared to optical PSF
@param pixel_effectInit sigma describing charge diffusion effect [in units of 15 microns]
@param fiber_rInit radius of perfect tophat fiber, as seen on the detector [in units of 15 microns]
@param fluxInit total flux in generated image compared to input image (needs to be 1 or very close to 1)
#not used anymore
@param dxInit (not used in this version of the code - parameter determing position of PSF on detector)
@param dyInit (not used in this version of the code - parameter determing position of PSF on detector )
@param apodizationInit (not used in this iteration of the code) by how much pixels to convolve the pupil image to apodize the strucutre - !
@param trace_valueInit (not used in this iteration of the code) inital value for adding vertical component to the data
@param serial_trace_valueInit (not used in this iteration of the code)inital value for adding horizontal component to the data
"""
if self.verbosity==1:
print(' ')
print('Initializing ZernikeFitter_PFS')
print('Verbosity parameter is: '+str(self.verbosity))
print('Highest Zernike polynomial is (zmax): '+str(self.zmax))
params = lmfit.Parameters()
z_array=[]
if z4Init is None:
params.add('z4', 0.0)
else:
params.add('z4', z4Init)
for i in range(5, self.zmax+1):
params.add('z{}'.format(i), 0.0)
if dxInit is None:
params.add('dx', 0.0)
else:
params.add('dx', dxInit)
if dyInit is None:
params.add('dy', 0.0)
else:
params.add('dy', dyInit)
if hscFracInit is None:
params.add('hscFrac', 0)
else:
params.add('hscFrac',hscFracInit)
if strutFracInit is None:
params.add('strutFrac', 0)
else:
params.add('strutFrac', strutFracInit)
if focalPlanePositionInit is None:
params.add('dxFocal', 0.0)
params.add('dyFocal', 0.0)
else:
params.add('dxFocal', focalPlanePositionInit[0])
params.add('dyFocal', focalPlanePositionInit[1])
if slitFracInit is None:
params.add('slitFrac', 0)
else:
params.add('slitFrac', slitFracInit)
if slitFrac_dy_Init is None:
params.add('slitFrac_dy', 0)
else:
params.add('slitFrac_dy', slitFrac_dy_Init)
if fiber_rInit is None:
params.add('fiber_r', 1.8)
else:
params.add('fiber_r', fiber_rInit)
if radiometricEffectInit is None:
params.add('radiometricEffect', 0)
else:
params.add('radiometricEffect', radiometricEffectInit)
if trace_valueInit is None:
params.add('trace_value', 0)
else:
params.add('trace_value', trace_valueInit)
if serial_trace_valueInit is None:
params.add('serial_trace_value', 0)
else:
params.add('serial_trace_value', serial_trace_valueInit)
if pixel_effectInit is None:
params.add('pixel_effect', 1)
else:
params.add('pixel_effect', pixel_effectInit)
if backgroundInit is None:
params.add('background', 0)
else:
params.add('background', backgroundInit)
if fluxInit is None:
params.add('flux', 1)
else:
params.add('flux', fluxInit)
if x_ilumInit is None:
params.add('x_ilum', 1)
else:
params.add('x_ilum', x_ilumInit)
if y_ilumInit is None:
params.add('y_ilum', 1)
else:
params.add('y_ilum', y_ilumInit)
if radiometricExponentInit is None:
params.add('radiometricExponent', 0.25)
else:
params.add('radiometricExponent', radiometricExponentInit)
if x_ilumInit is None:
params.add('x_fiber', 1)
else:
params.add('x_fiber', x_fiberInit)
if effective_ilum_radiusInit is None:
params.add('effective_ilum_radius', 1)
else:
params.add('effective_ilum_radius', effective_ilum_radiusInit)
if y_fiberInit is None:
params.add('y_fiber', 0)
else:
params.add('y_fiber', y_fiberInit)
if grating_linesInit is None:
params.add('grating_lines', 100000)
else:
params.add('grating_lines', grating_linesInit)
if scattering_slopeInit is None:
params.add('scattering_slope', 2)
else:
params.add('scattering_slope', scattering_slopeInit)
if scattering_amplitudeInit is None:
params.add('scattering_amplitude', 10**-2)
else:
params.add('scattering_amplitude', scattering_amplitudeInit)
if frd_sigmaInit is None:
params.add('frd_sigma', 0.02)
else:
params.add('frd_sigma', frd_sigmaInit)
if frd_lorentz_factorInit is None:
params.add('frd_lorentz_factor', 0.5)
else:
params.add('frd_lorentz_factor', frd_lorentz_factorInit)
if det_vertInit is None:
params.add('det_vert', 1)
else:
params.add('det_vert', det_vertInit)
if slitHolder_frac_dxInit is None:
params.add('slitHolder_frac_dx', 0)
else:
params.add('slitHolder_frac_dx', slitHolder_frac_dxInit)
self.params = params
self.optPsf=None
self.z_array=z_array
def constructModelImage_PFS_naturalResolution(self,params=None,shape=None,pixelScale=None,jacobian=None,use_optPSF=None,extraZernike=None,hashvalue=None):
"""Construct model image from parameters
@param params lmfit.Parameters object or python dictionary with
param values to use, or None to use self.params
@param pixelScale pixel scale in arcseconds to use for model image,
or None to use self.pixelScale.
@param shape (nx, ny) shape for model image, or None to use
the shape of self.maskedImage
@returns numpy array image with the same flux as the input image
"""
if self.verbosity==1:
print(' ')
print('Entering constructModelImage_PFS_naturalResolution')
if params is None:
params = self.params
if shape is None:
shape = self.image.shape
if pixelScale is None:
pixelScale = self.pixelScale
if jacobian is None:
jacobian = np.eye(2, dtype=np.float64)
else:
jacobian = self.jacobian
try:
v = params.valuesdict()
except AttributeError:
v = params
if hashvalue is None:
hashvalue = 0
use_optPSF=self.use_optPSF
if extraZernike is None:
extraZernike=None
self.extraZernike=extraZernike
else:
extraZernike=list(extraZernike)
self.extraZernike=extraZernike
# This give image in nyquist resolution
# if not explicitly stated to the full procedure
if use_optPSF is None:
optPsf=self._getOptPsf_naturalResolution(v)
else:
#if first iteration still generate image
if self.optPsf is None:
optPsf=self._getOptPsf_naturalResolution(v)
self.optPsf=optPsf
else:
optPsf=self.optPsf
optPsf_cut_fiber_convolved_downsampled=self.optPsf_postprocessing(optPsf)
if self.save==1:
if socket.gethostname()=='IapetusUSA':
np.save(TESTING_FINAL_IMAGES_FOLDER+'optPsf',optPsf)
np.save(TESTING_FINAL_IMAGES_FOLDER+'optPsf_cut_fiber_convolved_downsampled',optPsf_cut_fiber_convolved_downsampled)
else:
pass
if self.verbosity==1:
print('Finished with constructModelImage_PFS_naturalResolution')
print(' ')
return optPsf_cut_fiber_convolved_downsampled
def optPsf_postprocessing(self,optPsf):
if self.verbosity==1:
print(' ')
print('Entering optPsf_postprocessing')
params = self.params
shape = self.image.shape
double_sources=self.double_sources
v = params.valuesdict()
# how much is my generated image oversampled compared to final image
oversampling_original=(self.pixelScale)/self.scale_ModelImage_PFS_naturalResolution
if self.verbosity==1:
print('optPsf.shape'+str(optPsf.shape))
print('oversampling_original:' +str(oversampling_original))
# from the large image cut the central portion (1.4 times larger than the size of actual image)
size_of_central_cut=int(oversampling_original*self.image.shape[0]*1.4)
assert size_of_central_cut<optPsf.shape[0]
if self.verbosity==1:
print('size_of_central_cut: '+str(size_of_central_cut))
#cut part which you need image
optPsf_cut=Psf_position.cut_Centroid_of_natural_resolution_image(image=optPsf,size_natural_resolution=size_of_central_cut+1,oversampling=1,dx=0,dy=0)
if self.verbosity==1:
print('optPsf_cut.shape'+str(optPsf_cut.shape))
# reduce oversampling by the factor of 4 to make things easier
oversampling=np.round(oversampling_original/4)
if self.verbosity==1:
print('oversampling:' +str(oversampling))
size_of_optPsf_cut_downsampled=np.round(size_of_central_cut/(oversampling_original/oversampling))
if self.verbosity==1:
print('optPsf_cut.shape[0]'+str(optPsf_cut.shape[0]))
print('size_of_optPsf_cut_downsampled: '+str(size_of_optPsf_cut_downsampled))
#optPsf_cut_downsampled_p1=resize(optPsf_cut,(size_of_optPsf_cut_downsampled,size_of_optPsf_cut_downsampled))
#optPsf_cut_downsampled_m1=resize(optPsf_cut[1:-1,1:-1],(size_of_optPsf_cut_downsampled,size_of_optPsf_cut_downsampled))
#optPsf_cut_downsampled=(optPsf_cut_downsampled_p1+optPsf_cut_downsampled_m1)/2
optPsf_cut_downsampled=skimage.transform.resize(optPsf_cut,(size_of_optPsf_cut_downsampled+1,size_of_optPsf_cut_downsampled+1),mode='constant')
if self.verbosity==1:
print('optPsf_cut_downsampled.shape: '+str(optPsf_cut_downsampled.shape))
# ensure it is shape is even nubmer, needed for fft convolutions later
#if optPsf_cut_downsampled.shape[0] % 2 ==0:
# pass
#else:
# optPsf_cut_downsampled=optPsf_cut_downsampled[:optPsf_cut_downsampled.shape[0]-1,:optPsf_cut_downsampled.shape[0]-1]
# gives middle point of the image to used for calculations of scattered light
mid_point_of_optPsf_cut_downsampled=int(optPsf_cut_downsampled.shape[0]/2)
# gives the size of one pixel in optPsf_downsampled in microns
size_of_pixels_in_optPsf_cut_downsampled=(15/self.dithering)/oversampling