COMPASS  5.0.0
End-to-end AO simulation tool using GPU acceleration
geom_init.py
1 
37 
38 import shesha.config as conf
39 import shesha.constants as scons
40 from shesha.constants import CONST
41 
42 import shesha.util.make_pupil as mkP
43 import shesha.util.utilities as util
44 from shesha.sutra_wrap import carmaWrap_context, Telescope
45 from shesha.constants import ApertureType
46 import numpy as np
47 
48 
49 def tel_init(context: carmaWrap_context, p_geom: conf.Param_geom, p_tel: conf.Param_tel,
50  r0=None, ittime=None, p_wfss=None, dm=None):
51  """
52  Initialize the overall geometry of the AO system, including pupil and WFS
53 
54  :parameters:
55  context: (carmaWrap_context) : context
56  p_geom: (Param_geom) : geom settings
57  p_tel: (Param_tel) : telescope settings
58  r0: (float) : atmos r0 @ 0.5 microns
59  ittime: (float) : 1/loop frequency [s]
60  p_wfss: (list of Param_wfs) : wfs settings
61  dm: (list of Param_dm) : (optional) dms settings [=None]
62  :return:
63  telescope: (Telescope): Telescope object
64 
65  """
66  if p_wfss is not None:
67  # WFS geometry
68  nsensors = len(p_wfss)
69 
70  any_sh = [o.type for o in p_wfss].count(scons.WFSType.SH) > 0
71  # dm = None
72  if (p_wfss[0].dms_seen is None and dm is not None):
73  for i in range(nsensors):
74  if (not p_wfss[i].open_loop):
75  p_wfss[i].set_dms_seen(np.arange(len(dm), dtype=np.int32))
76 
77  # first get the wfs with max # of subaps
78  # we'll derive the geometry from the requirements in terms of sampling
79  if (any_sh):
80  indmax = np.argsort([o.nxsub for o in p_wfss
81  if o.type == scons.WFSType.SH])[-1]
82  else:
83  indmax = np.argsort([o.nxsub for o in p_wfss])[-1]
84 
85  print("*-----------------------")
86  print("Computing geometry of WFS", indmax)
87 
88  init_wfs_geom(p_wfss[indmax], r0, p_tel, p_geom, ittime, verbose=1)
89  # #do the same for other wfs
90  for i in range(nsensors):
91  if (i != indmax):
92  print("*-----------------------")
93  print("Computing geometry of WFS", i)
94  init_wfs_geom(p_wfss[i], r0, p_tel, p_geom, ittime, verbose=1)
95 
96  else:
97  geom_init(p_geom, p_tel)
98 
99  telescope = Telescope(context, p_geom._spupil.shape[0],
100  np.where(p_geom._spupil > 0)[0].size,
101  (p_geom._spupil * p_geom._apodizer).astype(np.float32),
102  p_geom._mpupil.shape[0], p_geom._mpupil)
103 
104  if (p_geom._phase_ab_M1 is not None):
105  telescope.set_phase_ab_M1(p_geom._phase_ab_M1)
106  telescope.set_phase_ab_M1_m(p_geom._phase_ab_M1_m)
107 
108  return telescope
109 
110 
111 def init_wfs_geom(p_wfs: conf.Param_wfs, r0: float, p_tel: conf.Param_tel,
112  p_geom: conf.Param_geom, ittime: float, verbose=1):
113  """Compute the geometry of WFSs: valid subaps, positions of the subaps,
114  flux per subap, etc...
115 
116  :parameters:
117  p_wfs: (Param_wfs) : wfs settings
118 
119  r0: (float) : atmos r0 @ 0.5 microns
120 
121  p_tel: (Param_tel) : telescope settings
122 
123  geom: (Param_geom) : geom settings
124 
125  ittime: (float) : 1/loop frequency [s]
126 
127  verbose: (int) : (optional) display informations if 0
128 
129 
130  """
131 
132  if p_geom.pupdiam:
133  if p_wfs.type == scons.WFSType.SH or p_wfs.type == scons.WFSType.PYRHR or p_wfs.type == scons.WFSType.PYRLR:
134  pdiam = p_geom.pupdiam // p_wfs.nxsub
135  if (p_geom.pupdiam % p_wfs.nxsub > 0):
136  pdiam += 1
137 
138  else:
139  pdiam = -1
140 
141  p_wfs._pdiam = pdiam
142 
143  init_wfs_size(p_wfs, r0, p_tel, verbose)
144 
145  if not p_geom.is_init:
146  # this is the wfs with largest # of subaps
147  # the overall geometry is deduced from it
148  #if not p_geom.pupdiam:
149  if p_geom.pupdiam != 0 and p_geom.pupdiam != p_wfs._pdiam * p_wfs.nxsub:
150  print("WARNING: pupdiam set value not correct")
151  p_geom.pupdiam = p_wfs._pdiam * p_wfs.nxsub
152  print("pupdiam used: ", p_geom.pupdiam)
153  if p_wfs.type == scons.WFSType.PYRHR or p_wfs.type == scons.WFSType.PYRLR:
154  geom_init(p_geom, p_tel, padding=p_wfs._nrebin)
155  elif (p_wfs.type == scons.WFSType.SH):
156  geom_init(p_geom, p_tel)
157  else:
158  raise RuntimeError("This WFS can not be used")
159 
160  if (p_wfs.type == scons.WFSType.PYRHR or p_wfs.type == scons.WFSType.PYRLR):
161  init_pyrhr_geom(p_wfs, r0, p_tel, p_geom, ittime, verbose=1)
162  elif (p_wfs.type == scons.WFSType.SH):
163  init_sh_geom(p_wfs, r0, p_tel, p_geom, ittime, verbose=1)
164  else:
165  raise RuntimeError("This WFS can not be used")
166 
167 
168 def init_wfs_size(p_wfs: conf.Param_wfs, r0: float, p_tel: conf.Param_tel, verbose=1):
169  """Compute all the parameters usefull for further WFS image computation (array sizes)
170 
171  :parameters:
172  p_wfs: (Param_wfs) : wfs settings
173 
174  r0: (float) : atmos r0 @ 0.5 microns
175 
176  p_tel: (Param_tel) : telescope settings
177 
178  verbose: (int) : (optional) display informations if 0
179 
180  Scheme to determine arrays sizes
181  sh :
182  k = 6
183  p = k * d/r0
184  n = int(2*d*v/lambda/CONST.RAD2ARCSEC)+1
185  N = fft_goodsize(k*n/v*lambda/r0*CONST.RAD2ARCSEC)
186  u = k * lambda / r0 * CONST.RAD2ARCSEC / N
187  n = v/u - int(v/u) > 0.5 ? int(v/u)+1 : int(v/u)
188  v = n * u
189  Nt = v * Npix
190 
191  pyr :
192  Fs = field stop radius in arcsec
193  N size of big array for FFT in pixels
194  P pupil diameter in pixels
195  D diameter of telescope in m
196  Nssp : number of pyr measurement points in the pupil
197 
198  Rf = Fs . N . D / lambda / P
199  ideally we choose : Fs = lambda / D . Nssp / 2
200 
201  if we want good sampling of r0 (avoid aliasing of speckles)
202  P > D / r0 . m
203  with m = 2 or 3
204 
205  to get reasonable space between pupil images : N > P.(2 + 3S)
206  with S close to 1
207  N must be a power of 2
208 
209  to ease computation lets assume that Nssp is a divider of P
210  scaling factor between high res pupil images in pyramid model
211  and actual pupil size on camera images would be P / Nssp
212 
213  """
214 
215  r0 = r0 * (p_wfs.Lambda * 2)**(6. / 5)
216 
217  if (r0 != 0):
218  if (verbose):
219  print("r0 for WFS :", "%3.2f" % r0, " m")
220  # seeing = CONST.RAD2ARCSEC * (p_wfs.lambda * 1.e-6) / r0
221  if (verbose):
222  print("seeing for WFS : ",
223  "%3.2f" % (CONST.RAD2ARCSEC * (p_wfs.Lambda * 1.e-6) / r0), "\"")
224 
225  if (p_wfs._pdiam <= 0):
226  # this case is usualy for the wfs with max # of subaps
227  # we look for the best compromise between pixsize and fov
228  if (p_wfs.type == scons.WFSType.SH):
229 
230  subapdiam = p_tel.diam / float(p_wfs.nxsub) # diam of subap
231  k = 6
232  pdiam = int(k * subapdiam / r0) # number of phase points per subap
233  if (pdiam < 16):
234  pdiam = 16
235 
236  # Must be even to keep ssp and actuators grids aligned in the pupil
237  if ((pdiam * p_wfs.nxsub) % 2):
238  pdiam += 1
239 
240  nrebin = int(2 * subapdiam * p_wfs.pixsize /
241  (p_wfs.Lambda * 1.e-6) / CONST.RAD2ARCSEC) + 1
242  nrebin = max(2, nrebin)
243  # first atempt on a rebin factor
244 
245  # since we clipped pdiam we have to be carreful in nfft computation
246  Nfft = util.fft_goodsize(
247  int(pdiam / subapdiam * nrebin / p_wfs.pixsize * CONST.RAD2ARCSEC *
248  (p_wfs.Lambda * 1.e-6)))
249 
250  elif (p_wfs.type == scons.WFSType.PYRHR or p_wfs.type == scons.WFSType.PYRLR):
251  # while (pdiam % p_wfs.npix != 0) pdiam+=1;
252  k = 3 if p_wfs.type == scons.WFSType.PYRHR else 2
253  pdiam = int(p_tel.diam / r0 * k)
254  while (pdiam % p_wfs.nxsub != 0):
255  pdiam += 1 # we choose to have a multiple of p_wfs.nxsub
256  pdiam = pdiam // p_wfs.nxsub
257  if (pdiam < 8):
258  pdiam = 8
259 
260  # quantum pixel size
261  else:
262  pdiam = p_wfs._pdiam
263  # this case is for a wfs with fixed # of phase points
264  Nfft = util.fft_goodsize(2 * pdiam)
265  # size of the support in fourier domain
266 
267  # qpixsize = pdiam * \
268  # (p_wfs.Lambda * 1.e-6) / subapdiam * CONST.RAD2ARCSEC / Nfft
269  # # quantum pixel size
270 
271  if (p_wfs.type == scons.WFSType.SH):
272  subapdiam = p_tel.diam / float(p_wfs.nxsub) # diam of subap
273 
274  # size of the support in fourier domain
275 
276  # qpixsize = k * (p_wfs.Lambda*1.e-6) / r0 * CONST.RAD2ARCSEC / Nfft
277  qpixsize = (pdiam * (p_wfs.Lambda * 1.e-6) / subapdiam * CONST.RAD2ARCSEC) / Nfft
278 
279  # # actual rebin factor
280  if (p_wfs.pixsize / qpixsize - int(p_wfs.pixsize / qpixsize) > 0.5):
281  nrebin = int(p_wfs.pixsize / qpixsize) + 1
282  else:
283  nrebin = int(p_wfs.pixsize / qpixsize)
284 
285  # actual pixel size
286  pixsize = nrebin * qpixsize
287 
288  if (pixsize * p_wfs.npix > qpixsize * Nfft):
289  Ntot = util.fft_goodsize(int(pixsize * p_wfs.npix / qpixsize) + 1)
290  else:
291  Ntot = Nfft
292 
293  if (Ntot % 2 != Nfft % 2):
294  Ntot += 1
295 
296  elif (p_wfs.type == scons.WFSType.PYRHR or p_wfs.type == scons.WFSType.PYRLR):
297  pdiam = pdiam * p_wfs.nxsub
298  m = 3
299  # fft_goodsize( m * pdiam)
300  Nfft = int(2**np.ceil(np.log2(m * pdiam)))
301 
302  nrebin = pdiam // p_wfs.nxsub
303  while (Nfft % nrebin != 0):
304  nrebin += 1 # we choose to have a divisor of Nfft
305  pdiam = nrebin * p_wfs.nxsub
306  Nfft = int(2**np.ceil(np.log2(m * pdiam)))
307 
308  qpixsize = (pdiam *
309  (p_wfs.Lambda * 1.e-6) / p_tel.diam * CONST.RAD2ARCSEC) / Nfft
310 
311  Ntot = Nfft // pdiam * p_wfs.nxsub
312  pixsize = qpixsize * nrebin
313  pdiam = pdiam // p_wfs.nxsub
314 
315  p_wfs._pdiam = pdiam
316  p_wfs.pixsize = pixsize
317  p_wfs._qpixsize = qpixsize
318  p_wfs._Nfft = Nfft
319  p_wfs._Ntot = Ntot
320  p_wfs._nrebin = nrebin
321  p_wfs._subapd = p_tel.diam / p_wfs.nxsub
322 
323  if (verbose):
324  if (p_wfs.type == scons.WFSType.SH):
325  print("quantum pixsize : ", "%5.4f" % qpixsize, "\"")
326  print("simulated FoV : ", "%3.2f" % (Ntot * qpixsize), "\" x ",
327  "%3.2f" % (Ntot * qpixsize), "\"")
328  print("actual pixsize : ", "%5.4f" % pixsize)
329  print("actual FoV : ", "%3.2f" % (pixsize * p_wfs.npix), "\" x ",
330  "%3.2f" % (pixsize * p_wfs.npix), "\"")
331  print("number of phase points : ", p_wfs._pdiam)
332  print("size of fft support : ", Nfft)
333  print("size of HR spot support : ", Ntot)
334 
335  elif (p_wfs.type == scons.WFSType.PYRHR or p_wfs.type == scons.WFSType.PYRLR):
336  print("quantum pixsize in pyr image : ", "%5.4f" % qpixsize, "\"")
337  print("simulated FoV : ", "%3.2f" % (Nfft * qpixsize), "\" x ",
338  "%3.2f" % (Nfft * qpixsize), "\"")
339  print("number of phase points : ", p_wfs._pdiam * p_wfs.nxsub)
340  print("size of fft support : ", Nfft)
341 
342 
343 def compute_nphotons(wfs_type, ittime, optthroughput, diam, cobs=0, nxsub=0, zerop=0,
344  gsmag=0, lgsreturnperwatt=0, laserpower=0, verbose=1):
345  ''' Determines the number of photons TBC
346 
347  :parameters:
348  wfs_type: (scons.WFSType) : wfs type: SH or PYRHR.
349 
350  ittime: (float) : 1/loop frequency [s].
351 
352  optthroughput: (float) : wfs global throughput.
353 
354  diam: (float) : telescope diameter.
355 
356  cobs: (float) : (optional for SH) telescope central obstruction.
357 
358  nxsub: (int) : (optional for PYRHR) linear number of subaps.
359 
360  zerop: (float) : (optional for LGS) detector zero point expressed in ph/m**2/s in the bandwidth of the WFS.
361 
362  gsmag: (float) : (optional for LGS) magnitude of guide star.
363 
364  lgsreturnperwatt: (float) : (optional for NGS) return per watt factor (high season : 10 ph/cm2/s/W).
365 
366  laserpower: (float) : (optional for NGS) laser power in W.
367 
368  verbose: (bool) : (optional) display informations if True.
369 
370  for PYRHR WFS: nphotons = compute_nphotons(scons.WFSType.PYRHR, ittime,
371  optthroughput, diam, cobs=?, zerop=?, gsmag=?)
372  for NGS SH WFS: nphotons = compute_nphotons(scons.WFSType.SH, ittime,
373  optthroughput, diam, nxsub=?, zerop=?, gsmag=?)
374  for LGS SH WFS: nphotons = compute_nphotons(scons.WFSType.SH, ittime,
375  optthroughput, diam, nxsub=?,
376  lgsreturnperwatt=?, laserpower=?)
377  '''
378  surface = 0
379  nphotons = 0
380  if (wfs_type == scons.WFSType.PYRHR or wfs_type == scons.WFSType.PYRLR):
381  surface = np.pi / 4. * (1 - cobs**2.) * diam**2.
382  elif (wfs_type == scons.WFSType.SH):
383  # from the guide star
384  if (laserpower == 0):
385  if (zerop == 0):
386  zerop = 1e11
387  surface = (diam / nxsub)**2.
388  # include throughput to WFS for unobstructed
389  # subaperture per iteration
390  else: # we are dealing with a LGS
391  nphotons = lgsreturnperwatt * laserpower * \
392  optthroughput * (diam / nxsub) ** 2. * 1e4 * ittime
393  # detected by WFS
394  # ... for given power include throughput to WFS
395  # for unobstructed subaperture per iteration
396  if (verbose):
397  print("nphotons : ", nphotons)
398  return nphotons
399  else:
400  raise RuntimeError("WFS unknown")
401 
402  nphotons = zerop * 10. ** (-0.4 * gsmag) * ittime * \
403  optthroughput * surface
404 
405  if (verbose):
406  print("nphotons : ", nphotons)
407  return nphotons
408 
409 
410 def init_pyrhr_geom(p_wfs: conf.Param_wfs, r0: float, p_tel: conf.Param_tel,
411  p_geom: conf.Param_geom, ittime: float, verbose: bool = True):
412  """Compute the geometry of PYRHR WFSs: valid subaps, positions of the subaps,
413  flux per subap, etc...
414 
415  :parameters:
416  p_wfs: (Param_wfs) : wfs settings
417 
418  r0: (float) : atmos r0 @ 0.5 microns
419 
420  p_tel: (Param_tel) : telescope settings
421 
422  geom: (Param_geom) : geom settings
423 
424  ittime: (float) : 1/loop frequency [s]
425 
426  verbose: (bool) : (optional) display informations if True
427 
428  """
429 
430  p_wfs.npix = p_wfs._pdiam
431  p_wfs._Ntot = p_geom._n
432 
433  # Creating field stop mask
434  fsradius_pixels = int(p_wfs.fssize / p_wfs._qpixsize / 2.)
435  if (p_wfs.fstop == scons.FieldStopType.ROUND):
436  focmask = util.dist(p_wfs._Nfft, xc=p_wfs._Nfft / 2. - 0.5,
437  yc=p_wfs._Nfft / 2. - 0.5) < (fsradius_pixels)
438  elif (p_wfs.fstop == scons.FieldStopType.SQUARE):
439  X = np.indices((p_wfs._Nfft, p_wfs._Nfft)) + 1 # TODO: +1 ??
440  x = X[1] - (p_wfs._Nfft + 1.) / 2.
441  y = X[0] - (p_wfs._Nfft + 1.) / 2.
442  focmask = (np.abs(x) <= (fsradius_pixels)) * \
443  (np.abs(y) <= (fsradius_pixels))
444  else:
445  msg = "PYRHR wfs fstop must be FieldStopType.[ROUND|SQUARE]"
446  raise ValueError(msg)
447 
448  pyr_focmask = focmask * 1.0 # np.fft.fftshift(focmask*1.0)
449  p_wfs._submask = np.fft.fftshift(pyr_focmask)
450 
451  # Creating pyramid mask
452  pyrsize = p_wfs._Nfft
453  cobs = p_tel.cobs
454  rpup = p_geom.pupdiam / 2.0
455  dpup = p_geom.pupdiam
456  nrebin = p_wfs._nrebin
457  fracsub = p_wfs.fracsub
458  if p_wfs.pyr_pup_sep == -1:
459  pup_sep = int(p_wfs.nxsub)
460  else:
461  pup_sep = p_wfs.pyr_pup_sep
462  # number of pix between two centers two pupil images
463 
464  y = np.tile(np.arange(pyrsize) - pyrsize / 2, (pyrsize, 1))
465  x = y.T
466 
467  Pangle = pup_sep * nrebin # Pyramid angle in HR pixels
468  if p_wfs.nPupils == 0:
469  p_wfs.nPupils = 4
470  # Centers is a nPupils x 2 array describing the position of the quadrant centers
471  centers = Pangle / np.sin(
472  np.pi / p_wfs.nPupils) * np.c_[np.cos(
473  (2 * np.arange(p_wfs.nPupils) + 1) * np.pi / p_wfs.nPupils),
474  np.sin((2 * np.arange(p_wfs.nPupils) + 1) *
475  np.pi / p_wfs.nPupils)]
476  # In case nPupils == 4, we put the centers back in the normal A-B-C-D ordering scheme, for misalignment processing.
477  if p_wfs.nPupils == 4:
478  centers = np.round(centers[[1, 0, 2, 3], :]).astype(np.int32)
479  # Add in the misalignments of quadrant centers, in LR px units
480  if p_wfs.pyr_misalignments is not None:
481  mis = np.asarray(p_wfs.pyr_misalignments) * nrebin
482  else:
483  mis = np.zeros((p_wfs.nPupils, 2), dtype=np.float32)
484  # Pyramid as minimal intersect of 4 tilting planes
485  pyr = 2 * np.pi / pyrsize * np.min(
486  np.asarray([(c[0] + m[0]) * x + (c[1] + m[1]) * y
487  for c, m in zip(centers, mis)]), axis=0)
488  p_wfs._halfxy = np.fft.fftshift(pyr.T)
489 
490  # Valid pixels identification
491  # Generate buffer with pupil at center
492  pup = np.zeros((pyrsize, pyrsize), dtype=np.int32)
493  pup[pyrsize // 2 - p_geom._n // 2:pyrsize // 2 + p_geom._n // 2,
494  pyrsize // 2 - p_geom._n // 2:pyrsize // 2 + p_geom._n // 2] = \
495  p_geom._mpupil
496 
497  # Shift the mask to obtain the geometrically valid pixels
498  for qIdx in range(p_wfs.nPupils):
499  quadOnCenter = np.roll(pup, tuple((centers[qIdx] + mis[qIdx]).astype(np.int32)),
500  (0, 1))
501  mskRebin = util.rebin(quadOnCenter.copy(),
502  [pyrsize // nrebin, pyrsize // nrebin])
503  if qIdx == 0:
504  stackedSubap = np.roll(mskRebin >= fracsub,
505  tuple((-centers[qIdx] / nrebin).astype(np.int32)),
506  (0, 1))
507  mskRebTot = (mskRebin >= fracsub).astype(np.int32)
508  else:
509  stackedSubap += np.roll(mskRebin >= fracsub,
510  tuple((-centers[qIdx] / nrebin).astype(np.int32)),
511  (0, 1))
512  mskRebTot += (mskRebin >= fracsub).astype(np.int32) * (qIdx + 1)
513 
514  validRow, validCol = [], []
515  # If n == 4, we need to tweak -again- the order to feed
516  # compass XY controllers in the appropriate order
517  if p_wfs.nPupils == 4:
518  centers = centers[[2, 1, 3, 0], :]
519  # mis = mis[[2, 1, 3, 0], :] # We don't need mis anymore - but if so keep the order of centers
520  for qIdx in range(p_wfs.nPupils):
521  tmpWh = np.where(
522  np.roll(stackedSubap, tuple((centers[qIdx] / nrebin).astype(np.int32)),
523  (0, 1)))
524  validRow += [tmpWh[0].astype(np.int32)]
525  validCol += [tmpWh[1].astype(np.int32)]
526  nvalid = validRow[0].size
527  validRow = np.asarray(validRow).flatten()
528  validCol = np.asarray(validCol).flatten()
529  stack, index = np.unique(np.c_[validRow, validCol], axis=0, return_index=True)
530 
531  p_wfs._nvalid = nvalid
532  p_wfs._validsubsx = validRow[np.sort(index)]
533  p_wfs._validsubsy = validCol[np.sort(index)]
534  p_wfs._hrmap = mskRebTot.astype(np.int32)
535 
536  if (p_wfs.pyr_pos is None):
537  pixsize = (np.pi * p_wfs._qpixsize) / (3600 * 180)
538  # scale_fact = 2 * np.pi / npup * \
539  # (p_wfs.Lambda / p_tel.diam / 4.848) / pixsize * p_wfs.pyr_ampl
540  # Proposition de Flo
541  scale_fact = 2 * np.pi / pyrsize * \
542  (p_wfs.Lambda * 1e-6 / p_tel.diam) / pixsize * p_wfs.pyr_ampl
543  cx = scale_fact * \
544  np.sin((np.arange(p_wfs.pyr_npts)) * 2. * np.pi / p_wfs.pyr_npts)
545  cy = scale_fact * \
546  np.cos((np.arange(p_wfs.pyr_npts)) * 2. * np.pi / p_wfs.pyr_npts)
547  p_wfs.set_pyr_scale_pos(scale_fact)
548  # mod_npts = p_wfs.pyr_npts #UNUSED
549  else:
550  if (verbose):
551  print("Using user-defined positions [arcsec] for the pyramid modulation")
552  cx = p_wfs.pyr_pos[:, 0] / p_wfs._qpixsize
553  cy = p_wfs.pyr_pos[:, 1] / p_wfs._qpixsize
554  # mod_npts=cx.shape[0] #UNUSED
555 
556  p_wfs._pyr_cx = cx.copy()
557  p_wfs._pyr_cy = cy.copy()
558 
559  # telSurf = np.pi / 4. * (1 - p_tel.cobs**2.) * p_tel.diam**2.
560  # p_wfs._nphotons = p_wfs.zerop * \
561  # 10. ** (-0.4 * p_wfs.gsmag) * ittime * \
562  # p_wfs.optthroughput * telSurf
563  p_wfs._nphotons = compute_nphotons(scons.WFSType.PYRHR, ittime, p_wfs.optthroughput,
564  p_tel.diam, cobs=p_tel.cobs, zerop=p_wfs.zerop,
565  gsmag=p_wfs.gsmag, verbose=verbose)
566 
567  # spatial filtering by the pixel extent:
568  # *2/2 intended. min should be 0.40 = sinc(0.5)^2.
569  y = np.tile(np.arange(pyrsize) - pyrsize // 2, (pyrsize, 1))
570  x = y.T
571  x = x * 1. / pyrsize
572  y = y * 1. / pyrsize
573  sincar = np.fft.fftshift(np.sinc(x) * np.sinc(y))
574 
575  # sincar = np.roll(np.pi*x*np.pi*y,x.shape[1],axis=1)
576  p_wfs._sincar = sincar.astype(np.float32)
577 
578  #pup = p_geom._mpupil
579  a = pyrsize // nrebin
580  b = p_geom._n // nrebin
581  pupvalid = stackedSubap[a // 2 - b // 2:a // 2 + b // 2, a // 2 - b // 2:a // 2 +
582  b // 2]
583  p_wfs._isvalid = pupvalid.T.astype(np.int32)
584 
585  validsubsx = np.where(pupvalid)[0].astype(np.int32)
586  validsubsy = np.where(pupvalid)[1].astype(np.int32)
587 
588  istart = np.arange(p_wfs.nxsub + 2) * p_wfs.npix
589  jstart = np.copy(istart)
590 
591  # sorting out valid subaps
592  fluxPerSub = np.zeros((p_wfs.nxsub + 2, p_wfs.nxsub + 2), dtype=np.float32)
593 
594  for i in range(p_wfs.nxsub + 2):
595  indi = istart[i]
596  for j in range(p_wfs.nxsub + 2):
597  indj = jstart[j]
598  fluxPerSub[i, j] = np.sum(
599  p_geom._mpupil[indi:indi + p_wfs.npix, indj:indj + p_wfs.npix])
600 
601  fluxPerSub = fluxPerSub / p_wfs.nxsub**2.
602 
603  p_wfs._fluxPerSub = fluxPerSub.copy()
604 
605  phasemap = np.zeros((p_wfs.npix * p_wfs.npix, p_wfs._nvalid), dtype=np.int32)
606  X = np.indices((p_geom._n, p_geom._n)) # we need c-like indice
607  tmp = X[1] + X[0] * p_geom._n
608 
609  pyrtmp = np.zeros((p_geom._n, p_geom._n), dtype=np.int32)
610 
611  for i in range(len(validsubsx)):
612  indi = istart[validsubsy[i]] # +2-1 (yorick->python
613  indj = jstart[validsubsx[i]]
614  phasemap[:, i] = tmp[indi:indi + p_wfs.npix, indj:indj + p_wfs.npix].flatten("C")
615  pyrtmp[indi:indi + p_wfs.npix, indj:indj + p_wfs.npix] = i
616 
617  p_wfs._phasemap = phasemap
618 
619  p_wfs._pyr_offsets = pyrtmp # pshift
620 
621 
622 def init_sh_geom(p_wfs: conf.Param_wfs, r0: float, p_tel: conf.Param_tel,
623  p_geom: conf.Param_geom, ittime: float, verbose: bool = True):
624  """Compute the geometry of SH WFSs: valid subaps, positions of the subaps,
625  flux per subap, etc...
626 
627  :parameters:
628  p_wfs: (Param_wfs) : wfs settings
629 
630  r0: (float) : atmos r0 @ 0.5 microns
631 
632  p_tel: (Param_tel) : telescope settings
633 
634  geom: (Param_geom) : geom settings
635 
636  ittime: (float) : 1/loop frequency [s]
637 
638  verbose: (bool) : (optional) display informations if True
639 
640  """
641  p_wfs.nPupils = 1
642  # this is the i,j index of lower left pixel of subap in _spupil
643  istart = ((np.linspace(0, p_geom.pupdiam, p_wfs.nxsub + 1))[:-1]).astype(np.int64)
644  # Translation in _mupil useful for raytracing
645  istart += 2
646  jstart = np.copy(istart)
647 
648  # sorting out valid subaps
649  fluxPerSub = np.zeros((p_wfs.nxsub, p_wfs.nxsub), dtype=np.float32)
650 
651  for i in range(p_wfs.nxsub):
652  indi = istart[i] # +2-1 (yorick->python)
653  for j in range(p_wfs.nxsub):
654  indj = jstart[j] # +2-1 (yorick->python)
655  fluxPerSub[i, j] = np.sum(
656  p_geom._mpupil[indi:indi + p_wfs._pdiam, indj:indj + p_wfs._pdiam])
657  # fluxPerSub[i,j] = np.where(p_geom._mpupil[indi:indi+pdiam,indj:indj+pdiam] > 0)[0].size
658 
659  fluxPerSub = fluxPerSub / p_wfs._pdiam**2.
660 
661  pupvalid = (fluxPerSub >= p_wfs.fracsub) * 1
662  p_wfs._isvalid = pupvalid.astype(np.int32)
663  p_wfs._nvalid = int(np.sum(p_wfs._isvalid))
664  p_wfs._fluxPerSub = fluxPerSub.copy()
665  validx = np.where(p_wfs._isvalid.T)[1].astype(np.int32)
666  validy = np.where(p_wfs._isvalid.T)[0].astype(np.int32)
667  p_wfs._validsubsx = validx
668  p_wfs._validsubsy = validy
669  p_wfs._validpuppixx = validx * p_wfs._pdiam + 2
670  p_wfs._validpuppixy = validy * p_wfs._pdiam + 2
671 
672  # this defines how we cut the phase into subaps
673  phasemap = np.zeros((p_wfs._pdiam * p_wfs._pdiam, p_wfs._nvalid), dtype=np.int32)
674 
675  X = np.indices((p_geom._n, p_geom._n))
676  tmp = X[1] + X[0] * p_geom._n
677 
678  n = p_wfs._nvalid
679  # for i in range(p_wfs._nvalid):
680  for i in range(n):
681  indi = istart[p_wfs._validsubsy[i]] # +2-1 (yorick->python)
682  indj = jstart[p_wfs._validsubsx[i]]
683  phasemap[:, i] = tmp[indi:indi + p_wfs._pdiam, indj:indj +
684  p_wfs._pdiam].flatten()
685  p_wfs._phasemap = phasemap
686  p_wfs._validsubsx *= p_wfs.npix
687  p_wfs._validsubsy *= p_wfs.npix
688 
689  # this is a phase shift of 1/2 pix in x and y
690  halfxy = np.linspace(0, 2 * np.pi, p_wfs._Nfft + 1)[0:p_wfs._pdiam] / 2.
691  halfxy = np.tile(halfxy, (p_wfs._pdiam, 1))
692  halfxy += halfxy.T
693  # p_wfs._halfxy = <float*>(halfxy*0.).data #dont work: half*0 is temp
694  # python obj
695 
696  if (p_wfs.npix % 2 == 1 and p_wfs._nrebin % 2 == 1):
697  # p_wfs._halfxy = <float*>(halfxy*0.)
698  halfxy = np.zeros((p_wfs._pdiam, p_wfs._pdiam), dtype=np.float32)
699  p_wfs._halfxy = halfxy.astype(np.float32)
700  else:
701  p_wfs._halfxy = halfxy.astype(np.float32)
702 
703  # this defines how we create a larger fov if required
704  if (p_wfs._Ntot != p_wfs._Nfft):
705  indi = int((p_wfs._Ntot - p_wfs._Nfft) / 2.) # +1 -1 (yorick>python)
706  indj = int(indi + p_wfs._Nfft)
707  X = np.indices((p_wfs._Nfft, p_wfs._Nfft)) + 1 # TODO: +1 ??
708  x = X[1]
709  y = X[0]
710  # hrpix
711  tmp = np.zeros((p_wfs._Ntot, p_wfs._Ntot))
712  tmp[indi:indj, indi:indj] = np.roll(x + (y - 1) * p_wfs._Nfft, p_wfs._Nfft // 2,
713  axis=0)
714  tmp[indi:indj, indi:indj] = np.roll(tmp[indi:indj, indi:indj], p_wfs._Nfft // 2,
715  axis=1)
716  # hrmap=roll(hrpix)
717  tmp = np.roll(tmp, p_wfs._Ntot // 2, axis=0)
718  tmp = np.roll(tmp, p_wfs._Ntot // 2, axis=1)
719 
720  tmp = np.where(tmp.flatten())[0]
721 
722  p_wfs._hrmap = np.copy(tmp.astype(np.int32))
723  p_wfs._hrmap = np.reshape(p_wfs._hrmap, (p_wfs._hrmap.shape[0], 1))
724  # must be set even if unused
725 
726  else:
727  tmp = np.zeros((2, 2))
728  p_wfs._hrmap = np.copy(tmp.astype(np.int32))
729  # must be set even if unused
730 
731  if (p_wfs._nrebin * p_wfs.npix % 2 != p_wfs._Ntot % 2):
732  # +2-1 yorick>python
733  indi = int((p_wfs._Ntot - p_wfs._nrebin * p_wfs.npix) / 2.) + 1
734  else:
735  indi = int((p_wfs._Ntot - p_wfs._nrebin * p_wfs.npix) / 2.) + 0
736  indj = int(indi + p_wfs._nrebin * p_wfs.npix)
737 
738  X = np.indices((p_wfs._nrebin * p_wfs.npix, p_wfs._nrebin * p_wfs.npix))
739  x = (X[1] / p_wfs._nrebin).astype(np.int64)
740  y = (X[0] / p_wfs._nrebin).astype(np.int64)
741 
742  # binindices
743  binindices = np.zeros((p_wfs._Ntot, p_wfs._Ntot))
744  binindices[indi:indj, indi:indj] = x + y * p_wfs.npix + 1
745 
746  binmap = np.zeros((p_wfs._nrebin * p_wfs._nrebin, p_wfs.npix * p_wfs.npix))
747 
748  X = np.indices((p_wfs._Ntot, p_wfs._Ntot))
749  tmp = X[1] + X[0] * p_wfs._Ntot
750 
751  if (p_wfs.gsalt <= 0):
752  binindices = np.roll(binindices, binindices.shape[0] // 2, axis=0)
753  binindices = np.roll(binindices, binindices.shape[1] // 2, axis=1)
754 
755  for i in range(p_wfs.npix * p_wfs.npix):
756  binmap[:, i] = tmp[np.where(binindices == i + 1)]
757  # binmap=np.reshape(binmap.flatten("F"),(binmap.shape[0],binmap.shape[1]))
758  p_wfs._binmap = np.copy(binmap.astype(np.int32))
759 
760  dr0 = p_tel.diam / r0 * \
761  (0.5 / p_wfs.Lambda) ** 1.2 / \
762  np.cos(p_geom.zenithangle * CONST.DEG2RAD) ** 0.6
763  fwhmseeing = p_wfs.Lambda / \
764  (p_tel.diam / np.sqrt(p_wfs.nxsub ** 2. + (dr0 / 1.5) ** 2.)) / 4.848
765  kernelfwhm = np.sqrt(fwhmseeing**2. + p_wfs.kernel**2.)
766 
767  tmp = util.makegaussian(p_wfs._Ntot, kernelfwhm / p_wfs._qpixsize, p_wfs._Ntot // 2,
768  p_wfs._Ntot // 2).astype(np.float32)
769 
770  tmp = np.roll(tmp, tmp.shape[0] // 2, axis=0)
771  tmp = np.roll(tmp, tmp.shape[1] // 2, axis=1)
772 
773  tmp[0, 0] = 1. # this insures that even with fwhm=0, the kernel is a dirac
774  tmp = tmp / np.sum(tmp)
775  tmp = np.fft.fft2(tmp).astype(np.complex64) / (p_wfs._Ntot * p_wfs._Ntot)
776  p_wfs._ftkernel = np.copy(tmp).astype(np.complex64)
777 
778  # dealing with photometry
779  # telSurf = np.pi / 4. * (1 - p_tel.cobs**2.) * p_tel.diam**2.
780 
781  # from the guide star
782  if (p_wfs.gsalt == 0):
783  if (p_wfs.zerop == 0):
784  p_wfs.zerop = 1e11
785  # p_wfs._nphotons = p_wfs.zerop * 10 ** (-0.4 * p_wfs.gsmag) * \
786  # p_wfs.optthroughput * \
787  # (p_tel.diam / p_wfs.nxsub) ** 2. * ittime
788  # include throughput to WFS
789  # for unobstructed subaperture
790  # per iteration
791  p_wfs._nphotons = compute_nphotons(scons.WFSType.SH, ittime, p_wfs.optthroughput,
792  p_tel.diam, nxsub=p_wfs.nxsub,
793  zerop=p_wfs.zerop, gsmag=p_wfs.gsmag,
794  verbose=verbose)
795 
796  else: # we are dealing with a LGS
797  # p_wfs._nphotons = p_wfs.lgsreturnperwatt * \
798  # p_wfs.laserpower * \
799  # p_wfs.optthroughput * \
800  # (p_tel.diam / p_wfs.nxsub) ** 2. * 1e4 * \
801  # ittime
802  # detected by WFS
803  # ... for given power
804  # include throughput to WFS
805  # for unobstructed subaperture
806  # per iteration
807  p_wfs._nphotons = compute_nphotons(scons.WFSType.SH, ittime, p_wfs.optthroughput,
808  p_tel.diam, nxsub=p_wfs.nxsub,
809  lgsreturnperwatt=p_wfs.lgsreturnperwatt,
810  laserpower=p_wfs.laserpower, verbose=verbose)
811 
812 
813 def geom_init(p_geom: conf.Param_geom, p_tel: conf.Param_tel, padding=2):
814  """
815  Initialize the system geometry
816 
817  :parameters:
818  p_geom: (Param_geom) : geometry settings
819  p_tel: (Param_tel) : telescope settings
820  padding: (optional) : padding factor for PYRHR geometry
821  """
822 
823  # First power of 2 greater than pupdiam
824  p_geom.ssize = int(2**np.ceil(np.log2(p_geom.pupdiam) + 1))
825  # Using images centered on 1/2 pixels
826  p_geom.cent = p_geom.ssize / 2 - 0.5
827 
828  p_geom._p1 = int(np.ceil(p_geom.cent - p_geom.pupdiam / 2.))
829  p_geom._p2 = int(np.floor(p_geom.cent + p_geom.pupdiam / 2.))
830 
831  p_geom.pupdiam = p_geom._p2 - p_geom._p1 + 1
832 
833  p_geom._n = p_geom.pupdiam + 2 * padding
834  p_geom._n1 = p_geom._p1 - padding
835  p_geom._n2 = p_geom._p2 + padding
836 
837  cent = p_geom.pupdiam / 2. - 0.5
838 
839  # Useful pupil
840  p_geom._spupil = mkP.make_pupil(p_geom.pupdiam, p_geom.pupdiam, p_tel, cent,
841  cent).astype(np.float32)
842 
843  # large pupil (used for image formation)
844  p_geom._ipupil = util.pad_array(p_geom._spupil, p_geom.ssize).astype(np.float32)
845 
846  # useful pupil + 4 pixels
847  p_geom._mpupil = util.pad_array(p_geom._spupil, p_geom._n).astype(np.float32)
848 
849  if (p_tel.std_piston and p_tel.std_tt):
850  p_geom._phase_ab_M1 = mkP.make_phase_ab(p_geom.pupdiam, p_geom.pupdiam, p_tel,
851  p_geom._spupil, cent,
852  cent).astype(np.float32)
853  p_geom._phase_ab_M1_m = util.pad_array(p_geom._phase_ab_M1,
854  p_geom._n).astype(np.float32)
855  #TODO: apodizer
856  """
857  if p_geom.apod:
858  if p_geom.apodFile is None or p_geom.apodFile == '':
859  apod_filename = shesha_savepath + \
860  "apodizer/SP_HARMONI_I4_C6_N1024.npy"
861  p_geom._apodizer = makeA.make_apodizer(
862  p_geom.pupdiam, p_geom.pupdiam,
863  apod_filename.encode(), 180. / 12.).astype(np.float32)
864  else:
865  """
866  p_geom._apodizer = np.ones(p_geom._spupil.shape, dtype=np.int32)
867  p_geom._pixsize = p_tel.diam / p_geom.pupdiam
868  p_geom.is_init = True
869 
870 
871 def geom_init_generic(p_geom, pupdiam, t_spiders=0.01, spiders_type="six", xc=0, yc=0,
872  real=0, cobs=0):
873  """Initialize the system geometry
874 
875  :parameters:
876  pupdiam: (long) : linear size of total pupil
877 
878  t_spiders: (float) : secondary supports ratio.
879 
880  spiders_type: (str) : secondary supports type: "four" or "six".
881 
882  xc: (int)
883 
884  yc: (int)
885 
886  real: (int)
887 
888  cobs: (float) : central obstruction ratio.
889  """
890  # Initialize the system pupil
891  # first poxer of 2 greater than pupdiam
892  p_geom.ssize = int(2**np.ceil(np.log2(pupdiam) + 1))
893  # using images centered on 1/2 pixels
894  p_geom.cent = p_geom.ssize / 2 + 0.5
895  # valid pupil geometry
896  pupdiam = int(pupdiam)
897  p_geom._p1 = int(np.ceil(p_geom.cent - pupdiam / 2.))
898  p_geom._p2 = int(np.floor(p_geom.cent + pupdiam / 2.))
899  p_geom.pupdiam = p_geom._p2 - p_geom._p1 + 1
900  p_geom._n = p_geom.pupdiam + 4
901  p_geom._n1 = p_geom._p1 - 2
902  p_geom._n2 = p_geom._p2 + 2
903 
904  # useful pupil
905  p_geom._spupil = mkP.make_pupil_generic(pupdiam, pupdiam, t_spiders, spiders_type,
906  xc, yc, real, cobs)
907 
908  # large pupil (used for image formation)
909  p_geom._ipupil = pad_array(p_geom._spupil, p_geom.ssize).astype(np.float32)
910 
911  # useful pupil + 4 pixels
912  p_geom._mpupil = pad_array(p_geom._spupil, p_geom._n).astype(np.float32)
913 
914 
915 def pad_array(A, N):
916  S = A.shape
917  D1 = (N - S[0]) // 2
918  D2 = (N - S[1]) // 2
919  padded = np.zeros((N, N))
920  padded[D1:D1 + S[0], D2:D2 + S[1]] = A
921  return padded
shesha.sutra_wrap
Definition: sutra_wrap.py:1
shesha.util.utilities
Basic utilities function.
Definition: utilities.py:1
shesha.constants
Numerical constants for shesha and config enumerations for safe-typing.
Definition: constants.py:1
shesha.util.make_pupil
Pupil creation functions.
Definition: make_pupil.py:1
shesha.init.geom_init.init_pyrhr_geom
def init_pyrhr_geom(conf.Param_wfs p_wfs, float r0, conf.Param_tel p_tel, conf.Param_geom p_geom, float ittime, bool verbose=True)
Compute the geometry of PYRHR WFSs: valid subaps, positions of the subaps, flux per subap,...
Definition: geom_init.py:427
shesha.config
Parameter classes for COMPASS.
Definition: shesha/shesha/config/__init__.py:1
shesha.init.geom_init.geom_init
def geom_init(conf.Param_geom p_geom, conf.Param_tel p_tel, padding=2)
Initialize the system geometry.
Definition: geom_init.py:821
shesha.init.geom_init.compute_nphotons
def compute_nphotons(wfs_type, ittime, optthroughput, diam, cobs=0, nxsub=0, zerop=0, gsmag=0, lgsreturnperwatt=0, laserpower=0, verbose=1)
Determines the number of photons TBC.
Definition: geom_init.py:376
shesha.init.geom_init.init_sh_geom
def init_sh_geom(conf.Param_wfs p_wfs, float r0, conf.Param_tel p_tel, conf.Param_geom p_geom, float ittime, bool verbose=True)
Compute the geometry of SH WFSs: valid subaps, positions of the subaps, flux per subap,...
Definition: geom_init.py:639
shesha.init.geom_init.init_wfs_size
def init_wfs_size(conf.Param_wfs p_wfs, float r0, conf.Param_tel p_tel, verbose=1)
Compute all the parameters usefull for further WFS image computation (array sizes)
Definition: geom_init.py:213
shesha.init.geom_init.pad_array
def pad_array(A, N)
Definition: geom_init.py:915
shesha.init.geom_init.init_wfs_geom
def init_wfs_geom(conf.Param_wfs p_wfs, float r0, conf.Param_tel p_tel, conf.Param_geom p_geom, float ittime, verbose=1)
Compute the geometry of WFSs: valid subaps, positions of the subaps, flux per subap,...
Definition: geom_init.py:129
shesha.init.geom_init.tel_init
def tel_init(carmaWrap_context context, conf.Param_geom p_geom, conf.Param_tel p_tel, r0=None, ittime=None, p_wfss=None, dm=None)
Initialize the overall geometry of the AO system, including pupil and WFS.
Definition: geom_init.py:64
shesha.init.geom_init.geom_init_generic
def geom_init_generic(p_geom, pupdiam, t_spiders=0.01, spiders_type="six", xc=0, yc=0, real=0, cobs=0)
Initialize the system geometry.
Definition: geom_init.py:888