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):
52 Initialize the overall geometry of the AO system, including pupil and WFS
55 context: (carmaWrap_context) : context
57 p_geom: (Param_geom) : geom settings
59 p_tel: (Param_tel) : telescope settings
61 r0: (float) : atmos r0 @ 0.5 microns
63 ittime: (float) : 1/loop frequency [s]
65 p_wfss: (list of Param_wfs) : wfs settings
68 dm: (list of Param_dm) : dms settings [=None]
71 telescope: (Telescope): Telescope object
74 if p_wfss
is not None:
76 nsensors = len(p_wfss)
78 any_sh = [o.type
for o
in p_wfss].count(scons.WFSType.SH) > 0
80 if (p_wfss[0].dms_seen
is None and dm
is not None):
81 for i
in range(nsensors):
82 if (
not p_wfss[i].open_loop):
83 p_wfss[i].set_dms_seen(np.arange(len(dm), dtype=np.int32))
88 indmax = np.argsort([o.nxsub
for o
in p_wfss
89 if o.type == scons.WFSType.SH])[-1]
91 indmax = np.argsort([o.nxsub
for o
in p_wfss])[-1]
93 print(
"*-----------------------")
94 print(
"Computing geometry of WFS", indmax)
96 init_wfs_geom(p_wfss[indmax], r0, p_tel, p_geom, ittime, verbose=1)
98 for i
in range(nsensors):
100 print(
"*-----------------------")
101 print(
"Computing geometry of WFS", i)
102 init_wfs_geom(p_wfss[i], r0, p_tel, p_geom, ittime, verbose=1)
107 telescope = Telescope(context, p_geom._spupil.shape[0],
108 np.where(p_geom._spupil > 0)[0].size,
109 (p_geom._spupil * p_geom._apodizer).astype(np.float32),
110 p_geom._mpupil.shape[0], p_geom._mpupil)
112 if (p_geom._phase_ab_M1
is not None):
113 telescope.set_phase_ab_M1(p_geom._phase_ab_M1)
114 telescope.set_phase_ab_M1_m(p_geom._phase_ab_M1_m)
119 def init_wfs_geom(p_wfs: conf.Param_wfs, r0: float, p_tel: conf.Param_tel,
120 p_geom: conf.Param_geom, ittime: float, verbose=1):
121 """Compute the geometry of WFSs: valid subaps, positions of the subaps,
122 flux per subap, etc...
125 p_wfs: (Param_wfs) : wfs settings
127 r0: (float) : atmos r0 @ 0.5 microns
129 p_tel: (Param_tel) : telescope settings
131 geom: (Param_geom) : geom settings
133 ittime: (float) : 1/loop frequency [s]
135 verbose: (int) : (optional) display informations if 0
141 if p_wfs.type == scons.WFSType.SH
or p_wfs.type == scons.WFSType.PYRHR
or p_wfs.type == scons.WFSType.PYRLR:
142 pdiam = p_geom.pupdiam // p_wfs.nxsub
143 if ((pdiam * p_wfs.nxsub) % 2):
152 if not p_geom.is_init:
156 if p_geom.pupdiam != 0
and p_geom.pupdiam != p_wfs._pdiam * p_wfs.nxsub:
157 print(
"WARNING: pupdiam set value not correct")
158 p_geom.pupdiam = p_wfs._pdiam * p_wfs.nxsub
159 print(
"pupdiam used: ", p_geom.pupdiam)
160 if p_wfs.type == scons.WFSType.PYRHR
or p_wfs.type == scons.WFSType.PYRLR:
161 geom_init(p_geom, p_tel, padding=p_wfs._nrebin)
162 elif (p_wfs.type == scons.WFSType.SH):
165 raise RuntimeError(
"This WFS can not be used")
167 if (p_wfs.type == scons.WFSType.PYRHR
or p_wfs.type == scons.WFSType.PYRLR):
169 elif (p_wfs.type == scons.WFSType.SH):
170 init_sh_geom(p_wfs, r0, p_tel, p_geom, ittime, verbose=1)
172 raise RuntimeError(
"This WFS can not be used")
175 def init_wfs_size(p_wfs: conf.Param_wfs, r0: float, p_tel: conf.Param_tel, verbose=1):
176 """Compute all the parameters usefull for further WFS image computation (array sizes)
179 p_wfs: (Param_wfs) : wfs settings
181 r0: (float) : atmos r0 @ 0.5 microns
183 p_tel: (Param_tel) : telescope settings
185 verbose: (int) : (optional) display informations if 0
187 Scheme to determine arrays sizes
190 p = k * d/r0 # size of seeing blob
191 n = int(2*d*v/lambda/CONST.RAD2ARCSEC)+1
192 N = fft_goodsize(k*n/v*lambda/r0*CONST.RAD2ARCSEC)
193 u = k * lambda / r0 * CONST.RAD2ARCSEC / N
194 n = v/u - int(v/u) > 0.5 ? int(v/u)+1 : int(v/u)
199 Fs = field stop radius in arcsec
200 N size of big array for FFT in pixels
201 P pupil diameter in pixels
202 D diameter of telescope in m
203 Nssp : number of pyr measurement points in the pupil
205 # Rf is the radius of field stop in pixels
206 Rf = Fs . N . D / lambda / P
207 ideally we choose : Fs = lambda / D . Nssp / 2
209 If we want a good sampling of r0 (to avoid aliasing of speckles that may
210 roll over the FoV), we have to specify a FoV of the FFT support that is
211 larger than seeing by a factor <m> at least equal to 2 (at a very minimum)
212 or 3 (for a better comfort..). This writes as:
214 with m = 2 or 3. This condition is equivalent to put <m> pixels per r0.
216 To get reasonable space between pupil images : N > P.(2 + 3S)
218 N must be a power of 2
220 to ease computation lets assume that Nssp is a divider of P
221 scaling factor between high res pupil images in pyramid model
222 and actual pupil size on camera images would be P / Nssp
225 r0 = r0 * (p_wfs.Lambda * 2)**(6. / 5)
229 print(
"r0 for WFS :",
"%4.3f" % r0,
" m")
232 print(
"seeing for WFS : ",
233 "%4.3f" % (CONST.RAD2ARCSEC * (p_wfs.Lambda * 1.e-6) / r0),
"\"")
235 if (p_wfs._pdiam <= 0):
238 if (p_wfs.type == scons.WFSType.SH):
240 subapdiam = p_tel.diam / float(p_wfs.nxsub)
242 pdiam = int(k * subapdiam / r0)
247 if ((pdiam * p_wfs.nxsub) % 2):
250 nrebin = int(2 * subapdiam * p_wfs.pixsize /
251 (p_wfs.Lambda * 1.e-6) / CONST.RAD2ARCSEC) + 1
252 nrebin = max(2, nrebin)
256 Nfft = util.fft_goodsize(
257 int(pdiam / subapdiam * nrebin / p_wfs.pixsize * CONST.RAD2ARCSEC *
258 (p_wfs.Lambda * 1.e-6)))
260 elif (p_wfs.type == scons.WFSType.PYRHR
or p_wfs.type == scons.WFSType.PYRLR):
262 k = 3
if p_wfs.type == scons.WFSType.PYRHR
else 2
263 pdiam = int(p_tel.diam / r0 * k)
264 while (pdiam % p_wfs.nxsub != 0):
266 pdiam = pdiam // p_wfs.nxsub
274 Nfft = util.fft_goodsize(2 * pdiam)
281 if (p_wfs.type == scons.WFSType.SH):
282 subapdiam = p_tel.diam / float(p_wfs.nxsub)
287 qpixsize = (pdiam * (p_wfs.Lambda * 1.e-6) / subapdiam * CONST.RAD2ARCSEC) / Nfft
290 if (p_wfs.pixsize / qpixsize - int(p_wfs.pixsize / qpixsize) > 0.5):
291 nrebin = int(p_wfs.pixsize / qpixsize) + 1
293 nrebin = int(p_wfs.pixsize / qpixsize)
296 pixsize = nrebin * qpixsize
298 if (pixsize * p_wfs.npix > qpixsize * Nfft):
299 Ntot = util.fft_goodsize(int(pixsize * p_wfs.npix / qpixsize) + 1)
303 if (Ntot % 2 != Nfft % 2):
306 elif (p_wfs.type == scons.WFSType.PYRHR
or p_wfs.type == scons.WFSType.PYRLR):
307 pdiam = pdiam * p_wfs.nxsub
310 Nfft = int(2**np.ceil(np.log2(m * pdiam)))
312 nrebin = pdiam // p_wfs.nxsub
313 while (Nfft % nrebin != 0):
315 pdiam = nrebin * p_wfs.nxsub
316 Nfft = int(2**np.ceil(np.log2(m * pdiam)))
320 (p_wfs.Lambda * 1.e-6) / p_tel.diam * CONST.RAD2ARCSEC) / Nfft
322 Ntot = Nfft // pdiam * p_wfs.nxsub
323 pixsize = qpixsize * nrebin
324 pdiam = pdiam // p_wfs.nxsub
327 p_wfs.pixsize = pixsize
328 p_wfs._qpixsize = qpixsize
331 p_wfs._nrebin = nrebin
332 p_wfs._subapd = p_tel.diam / p_wfs.nxsub
335 if (p_wfs.type == scons.WFSType.SH):
336 print(
"quantum pixsize : ",
"%5.4f" % qpixsize,
"\"")
337 print(
"simulated FoV : ",
"%3.2f" % (Ntot * qpixsize),
"\" x ",
338 "%3.2f" % (Ntot * qpixsize),
"\"")
339 print(
"actual pixsize : ",
"%5.4f" % pixsize)
340 print(
"actual FoV : ",
"%3.2f" % (pixsize * p_wfs.npix),
"\" x ",
341 "%3.2f" % (pixsize * p_wfs.npix),
"\"")
342 print(
"number of phase points : ", p_wfs._pdiam)
343 print(
"size of fft support : ", Nfft)
344 print(
"size of HR spot support : ", Ntot)
346 elif (p_wfs.type == scons.WFSType.PYRHR
or p_wfs.type == scons.WFSType.PYRLR):
347 print(
"quantum pixsize in pyr image : ",
"%5.4f" % qpixsize,
"\"")
348 print(
"simulated FoV : ",
"%3.2f" % (Nfft * qpixsize),
"\" x ",
349 "%3.2f" % (Nfft * qpixsize),
"\"")
350 print(
"number of phase points : ", p_wfs._pdiam * p_wfs.nxsub)
351 print(
"size of fft support : ", Nfft)
354 def compute_nphotons(wfs_type, ittime, optthroughput, diam, cobs=0, nxsub=0, zerop=0,
355 gsmag=0, lgsreturnperwatt=0, laserpower=0, verbose=1):
356 ''' Determines the number of photons TBC
359 wfs_type: (scons.WFSType) : wfs type: SH or PYRHR.
361 ittime: (float) : 1/loop frequency [s].
363 optthroughput: (float) : wfs global throughput.
365 diam: (float) : telescope diameter.
367 cobs: (float) : (optional for SH) telescope central obstruction.
369 nxsub: (int) : (optional for PYRHR) linear number of subaps.
371 zerop: (float) : (optional for LGS) detector zero point expressed in ph/m**2/s in the bandwidth of the WFS.
373 gsmag: (float) : (optional for LGS) magnitude of guide star.
375 lgsreturnperwatt: (float) : (optional for NGS) return per watt factor (high season : 10 ph/cm2/s/W).
377 laserpower: (float) : (optional for NGS) laser power in W.
379 verbose: (bool) : (optional) display informations if True.
381 for PYRHR WFS: nphotons = compute_nphotons(scons.WFSType.PYRHR, ittime,
382 optthroughput, diam, cobs=?, zerop=?, gsmag=?)
383 for NGS SH WFS: nphotons = compute_nphotons(scons.WFSType.SH, ittime,
384 optthroughput, diam, nxsub=?, zerop=?, gsmag=?)
385 for LGS SH WFS: nphotons = compute_nphotons(scons.WFSType.SH, ittime,
386 optthroughput, diam, nxsub=?,
387 lgsreturnperwatt=?, laserpower=?)
391 if (wfs_type == scons.WFSType.PYRHR
or wfs_type == scons.WFSType.PYRLR):
392 surface = np.pi / 4. * (1 - cobs**2.) * diam**2.
393 elif (wfs_type == scons.WFSType.SH):
395 if (laserpower == 0):
398 surface = (diam / nxsub)**2.
402 nphotons = lgsreturnperwatt * laserpower * \
403 optthroughput * (diam / nxsub) ** 2. * 1e4 * ittime
408 print(
"nphotons : ", nphotons)
411 raise RuntimeError(
"WFS unknown")
413 nphotons = zerop * 10. ** (-0.4 * gsmag) * ittime * \
414 optthroughput * surface
417 print(
"nphotons : ", nphotons)
421 def init_pyrhr_geom(p_wfs: conf.Param_wfs, r0: float, p_tel: conf.Param_tel,
422 p_geom: conf.Param_geom, ittime: float, verbose: bool =
True):
423 """Compute the geometry of PYRHR WFSs: valid subaps, positions of the subaps,
424 flux per subap, etc...
427 p_wfs: (Param_wfs) : wfs settings
429 r0: (float) : atmos r0 @ 0.5 microns
431 p_tel: (Param_tel) : telescope settings
433 geom: (Param_geom) : geom settings
435 ittime: (float) : 1/loop frequency [s]
437 verbose: (bool) : (optional) display informations if True
441 p_wfs.npix = p_wfs._pdiam
442 p_wfs._Ntot = p_geom._n
445 fsradius_pixels = int(p_wfs.fssize / p_wfs._qpixsize / 2.)
446 if (p_wfs.fstop == scons.FieldStopType.ROUND):
447 focmask = util.dist(p_wfs._Nfft, xc=p_wfs._Nfft / 2. - 0.5,
448 yc=p_wfs._Nfft / 2. - 0.5) < (fsradius_pixels)
449 elif (p_wfs.fstop == scons.FieldStopType.SQUARE):
450 X = np.indices((p_wfs._Nfft, p_wfs._Nfft)) + 1
451 x = X[1] - (p_wfs._Nfft + 1.) / 2.
452 y = X[0] - (p_wfs._Nfft + 1.) / 2.
453 focmask = (np.abs(x) <= (fsradius_pixels)) * \
454 (np.abs(y) <= (fsradius_pixels))
456 msg =
"PYRHR wfs fstop must be FieldStopType.[ROUND|SQUARE]"
457 raise ValueError(msg)
459 pyr_focmask = focmask * 1.0
460 p_wfs._submask = np.fft.fftshift(pyr_focmask)
463 pyrsize = p_wfs._Nfft
465 rpup = p_geom.pupdiam / 2.0
466 dpup = p_geom.pupdiam
467 nrebin = p_wfs._nrebin
468 fracsub = p_wfs.fracsub
469 if p_wfs.pyr_pup_sep == -1:
470 pup_sep = int(p_wfs.nxsub)
472 pup_sep = p_wfs.pyr_pup_sep
475 y = np.tile(np.arange(pyrsize) - pyrsize / 2, (pyrsize, 1))
478 Pangle = pup_sep * nrebin
479 if p_wfs.nPupils == 0:
482 centers = Pangle / np.sin(
483 np.pi / p_wfs.nPupils) * np.c_[np.cos(
484 (2 * np.arange(p_wfs.nPupils) + 1) * np.pi / p_wfs.nPupils),
485 np.sin((2 * np.arange(p_wfs.nPupils) + 1) *
486 np.pi / p_wfs.nPupils)]
488 if p_wfs.nPupils == 4:
489 centers = np.round(centers[[1, 0, 2, 3], :]).astype(np.int32)
491 if p_wfs.pyr_misalignments
is not None:
492 mis = np.asarray(p_wfs.pyr_misalignments) * nrebin
494 mis = np.zeros((p_wfs.nPupils, 2), dtype=np.float32)
496 pyr = 2 * np.pi / pyrsize * np.min(
497 np.asarray([(c[0] + m[0]) * x + (c[1] + m[1]) * y
498 for c, m
in zip(centers, mis)]), axis=0)
499 p_wfs._halfxy = np.fft.fftshift(pyr.T)
503 pup = np.zeros((pyrsize, pyrsize))
504 pup[pyrsize // 2 - p_geom._n // 2:pyrsize // 2 + p_geom._n // 2,
505 pyrsize // 2 - p_geom._n // 2:pyrsize // 2 + p_geom._n // 2] = \
509 for qIdx
in range(p_wfs.nPupils):
510 quadOnCenter = np.roll(pup, tuple((centers[qIdx] + mis[qIdx]).astype(np.int32)),
512 mskRebin = util.rebin(quadOnCenter.copy(),
513 [pyrsize // nrebin, pyrsize // nrebin])
515 stackedSubap = np.roll(mskRebin >= fracsub,
516 tuple((-centers[qIdx] / nrebin).astype(np.int32)),
518 mskRebTot = (mskRebin >= fracsub).astype(np.int32)
520 stackedSubap += np.roll(mskRebin >= fracsub,
521 tuple((-centers[qIdx] / nrebin).astype(np.int32)),
523 mskRebTot += (mskRebin >= fracsub).astype(np.int32) * (qIdx + 1)
525 validRow, validCol = [], []
528 if p_wfs.nPupils == 4:
529 centers = centers[[2, 1, 3, 0], :]
531 for qIdx
in range(p_wfs.nPupils):
533 np.roll(stackedSubap, tuple((centers[qIdx] / nrebin).astype(np.int32)),
535 validRow += [tmpWh[0].astype(np.int32)]
536 validCol += [tmpWh[1].astype(np.int32)]
537 nvalid = validRow[0].size
538 validRow = np.asarray(validRow).flatten()
539 validCol = np.asarray(validCol).flatten()
540 stack, index = np.unique(np.c_[validRow, validCol], axis=0, return_index=
True)
542 p_wfs._nvalid = nvalid
543 p_wfs._validsubsx = validRow[np.sort(index)]
544 p_wfs._validsubsy = validCol[np.sort(index)]
545 p_wfs._hrmap = mskRebTot.astype(np.int32)
547 if (p_wfs.pyr_pos
is None):
548 pixsize = (np.pi * p_wfs._qpixsize) / (3600 * 180)
552 scale_fact = 2 * np.pi / pyrsize * \
553 (p_wfs.Lambda * 1e-6 / p_tel.diam) / pixsize * p_wfs.pyr_ampl
555 np.sin((np.arange(p_wfs.pyr_npts)) * 2. * np.pi / p_wfs.pyr_npts)
557 np.cos((np.arange(p_wfs.pyr_npts)) * 2. * np.pi / p_wfs.pyr_npts)
558 p_wfs.set_pyr_scale_pos(scale_fact)
562 print(
"Using user-defined positions [arcsec] for the pyramid modulation")
563 cx = p_wfs.pyr_pos[:, 0] / p_wfs._qpixsize
564 cy = p_wfs.pyr_pos[:, 1] / p_wfs._qpixsize
567 p_wfs._pyr_cx = cx.copy()
568 p_wfs._pyr_cy = cy.copy()
574 p_wfs._nphotons =
compute_nphotons(scons.WFSType.PYRHR, ittime, p_wfs.optthroughput,
575 p_tel.diam, cobs=p_tel.cobs, zerop=p_wfs.zerop,
576 gsmag=p_wfs.gsmag, verbose=verbose)
580 y = np.tile(np.arange(pyrsize) - pyrsize // 2, (pyrsize, 1))
584 sincar = np.fft.fftshift(np.sinc(x) * np.sinc(y))
587 p_wfs._sincar = sincar.astype(np.float32)
590 a = pyrsize // nrebin
591 b = p_geom._n // nrebin
592 pupvalid = stackedSubap[a // 2 - b // 2:a // 2 + b // 2, a // 2 - b // 2:a // 2 +
594 p_wfs._isvalid = pupvalid.T.astype(np.int32)
596 validsubsx = np.where(pupvalid)[0].astype(np.int32)
597 validsubsy = np.where(pupvalid)[1].astype(np.int32)
599 istart = np.arange(p_wfs.nxsub + 2) * p_wfs.npix
600 jstart = np.copy(istart)
603 fluxPerSub = np.zeros((p_wfs.nxsub + 2, p_wfs.nxsub + 2), dtype=np.float32)
605 for i
in range(p_wfs.nxsub + 2):
607 for j
in range(p_wfs.nxsub + 2):
609 fluxPerSub[i, j] = np.sum(
610 p_geom._mpupil[indi:indi + p_wfs.npix, indj:indj + p_wfs.npix])
612 fluxPerSub = fluxPerSub / p_wfs.nxsub**2.
614 p_wfs._fluxPerSub = fluxPerSub.copy()
616 phasemap = np.zeros((p_wfs.npix * p_wfs.npix, p_wfs._nvalid), dtype=np.int32)
617 X = np.indices((p_geom._n, p_geom._n))
618 tmp = X[1] + X[0] * p_geom._n
620 pyrtmp = np.zeros((p_geom._n, p_geom._n), dtype=np.int32)
622 ttprojmat = np.zeros([p_wfs.npix * p_wfs.npix,p_wfs._nvalid*2], dtype=np.float32)
623 for i
in range(len(validsubsx)):
624 indi = istart[validsubsy[i]]
625 indj = jstart[validsubsx[i]]
626 phasemap[:, i] = tmp[indi:indi + p_wfs.npix, indj:indj + p_wfs.npix].flatten(
"C")
627 pyrtmp[indi:indi + p_wfs.npix, indj:indj + p_wfs.npix] = i
629 Sx = np.zeros([p_wfs.npix,p_wfs.npix],dtype=np.float32)
630 Sy = np.zeros([p_wfs.npix,p_wfs.npix],dtype=np.float32)
632 subapmask = p_geom._mpupil.T[indi:indi + p_wfs._pdiam,
633 indj:indj + p_wfs._pdiam]
634 for ii
in range(p_wfs.npix):
635 for jj
in range(p_wfs.npix):
636 Sx[ii,jj] += subapmask[ii,jj-1]
if jj>0
else 0
637 Sx[ii,jj] -= subapmask[ii,jj+1]
if jj+1 < p_wfs.npix
else 0
638 Sx[ii,jj] *= subapmask[ii,jj]
639 Sy[ii,jj] += subapmask[ii-1,jj]
if ii>0
else 0
640 Sy[ii,jj] -= subapmask[ii+1,jj]
if ii+1 < p_wfs.npix
else 0
641 Sy[ii,jj] *= subapmask[ii,jj]
642 Sx_den = -Sx.cumsum(axis=1).sum()
643 Sy_den = -Sy.cumsum(axis=0).sum()
644 if Sx_den == 0
or Sy_den == 0:
648 ttprojmat[:,i] = Sx.flatten()
649 ttprojmat[:,i+p_wfs._nvalid] = Sy.flatten()
651 p_wfs._phasemap = phasemap
652 p_wfs._ttprojmat = ttprojmat * p_geom.get_pupdiam() / p_tel.get_diam() \
653 * CONST.RAD2ARCSEC * 1e-6
655 p_wfs._pyr_offsets = pyrtmp
658 def init_sh_geom(p_wfs: conf.Param_wfs, r0: float, p_tel: conf.Param_tel,
659 p_geom: conf.Param_geom, ittime: float, verbose: bool =
True):
660 """Compute the geometry of SH WFSs: valid subaps, positions of the subaps,
661 flux per subap, etc...
664 p_wfs: (Param_wfs) : wfs settings
666 r0: (float) : atmos r0 @ 0.5 microns
668 p_tel: (Param_tel) : telescope settings
670 geom: (Param_geom) : geom settings
672 ittime: (float) : 1/loop frequency [s]
674 verbose: (bool) : (optional) display informations if True
679 istart = ((np.linspace(0, p_geom.pupdiam, p_wfs.nxsub + 1))[:-1]).astype(np.int64)
682 jstart = np.copy(istart)
685 fluxPerSub = np.zeros((p_wfs.nxsub, p_wfs.nxsub), dtype=np.float32)
687 for i
in range(p_wfs.nxsub):
689 for j
in range(p_wfs.nxsub):
691 fluxPerSub[i, j] = np.sum(
692 p_geom._mpupil[indi:indi + p_wfs._pdiam, indj:indj + p_wfs._pdiam])
695 fluxPerSub = fluxPerSub / p_wfs._pdiam**2.
697 pupvalid = (fluxPerSub >= p_wfs.fracsub) * 1
698 p_wfs._isvalid = pupvalid.astype(np.int32)
699 p_wfs._nvalid = int(np.sum(p_wfs._isvalid))
700 p_wfs._fluxPerSub = fluxPerSub.copy()
701 validx = np.where(p_wfs._isvalid.T)[1].astype(np.int32)
702 validy = np.where(p_wfs._isvalid.T)[0].astype(np.int32)
703 p_wfs._validsubsx = validx
704 p_wfs._validsubsy = validy
705 p_wfs._validpuppixx = validx * p_wfs._pdiam + 2
706 p_wfs._validpuppixy = validy * p_wfs._pdiam + 2
709 phasemap = np.zeros((p_wfs._pdiam * p_wfs._pdiam, p_wfs._nvalid), dtype=np.int32)
711 X = np.indices((p_geom._n, p_geom._n))
712 tmp = X[1] + X[0] * p_geom._n
716 ttprojmat = np.zeros([p_wfs._pdiam**2,p_wfs._nvalid*2], dtype=np.float32)
718 indi = istart[p_wfs._validsubsy[i]]
719 indj = jstart[p_wfs._validsubsx[i]]
720 phasemap[:, i] = tmp[indi:indi + p_wfs._pdiam, indj:indj +
721 p_wfs._pdiam].flatten()
723 Sx = np.zeros([p_wfs._pdiam,p_wfs._pdiam],dtype=np.float32)
724 Sy = np.zeros([p_wfs._pdiam,p_wfs._pdiam],dtype=np.float32)
726 subapmask = p_geom._mpupil.T[indi:indi + p_wfs._pdiam,
727 indj:indj + p_wfs._pdiam]
728 for ii
in range(p_wfs._pdiam):
729 for jj
in range(p_wfs._pdiam):
730 Sx[ii,jj] += subapmask[ii,jj-1]
if jj>0
else 0
731 Sx[ii,jj] -= subapmask[ii,jj+1]
if jj+1 < p_wfs._pdiam
else 0
732 Sx[ii,jj] *= subapmask[ii,jj]
733 Sy[ii,jj] += subapmask[ii-1,jj]
if ii>0
else 0
734 Sy[ii,jj] -= subapmask[ii+1,jj]
if ii+1 < p_wfs._pdiam
else 0
735 Sy[ii,jj] *= subapmask[ii,jj]
736 Sx_den = -Sx.cumsum(axis=1).sum()
737 Sy_den = -Sy.cumsum(axis=0).sum()
738 if Sx_den == 0
or Sy_den == 0:
742 ttprojmat[:,i] = Sx.flatten()
743 ttprojmat[:,i+p_wfs._nvalid] = Sy.flatten()
745 p_wfs._phasemap = phasemap
746 p_wfs._validsubsx *= p_wfs.npix
747 p_wfs._validsubsy *= p_wfs.npix
748 p_wfs._ttprojmat = ttprojmat * p_geom.get_pupdiam() / p_tel.get_diam() \
749 * CONST.RAD2ARCSEC * 1e-6
752 halfxy = np.linspace(0, 2 * np.pi, p_wfs._Nfft + 1)[0:p_wfs._pdiam] / 2.
753 halfxy = np.tile(halfxy, (p_wfs._pdiam, 1))
758 if (p_wfs.npix % 2 == 1
and p_wfs._nrebin % 2 == 1):
760 halfxy = np.zeros((p_wfs._pdiam, p_wfs._pdiam), dtype=np.float32)
761 p_wfs._halfxy = halfxy.astype(np.float32)
763 p_wfs._halfxy = halfxy.astype(np.float32)
766 if (p_wfs._Ntot != p_wfs._Nfft):
767 indi = int((p_wfs._Ntot - p_wfs._Nfft) / 2.)
768 indj = int(indi + p_wfs._Nfft)
769 X = np.indices((p_wfs._Nfft, p_wfs._Nfft)) + 1
773 tmp = np.zeros((p_wfs._Ntot, p_wfs._Ntot))
774 tmp[indi:indj, indi:indj] = np.roll(x + (y - 1) * p_wfs._Nfft, p_wfs._Nfft // 2,
776 tmp[indi:indj, indi:indj] = np.roll(tmp[indi:indj, indi:indj], p_wfs._Nfft // 2,
779 tmp = np.roll(tmp, p_wfs._Ntot // 2, axis=0)
780 tmp = np.roll(tmp, p_wfs._Ntot // 2, axis=1)
782 tmp = np.where(tmp.flatten())[0]
784 p_wfs._hrmap = np.copy(tmp.astype(np.int32))
785 p_wfs._hrmap = np.reshape(p_wfs._hrmap, (p_wfs._hrmap.shape[0], 1))
789 tmp = np.zeros((2, 2))
790 p_wfs._hrmap = np.copy(tmp.astype(np.int32))
793 if (p_wfs._nrebin * p_wfs.npix % 2 != p_wfs._Ntot % 2):
795 indi = int((p_wfs._Ntot - p_wfs._nrebin * p_wfs.npix) / 2.) + 1
797 indi = int((p_wfs._Ntot - p_wfs._nrebin * p_wfs.npix) / 2.) + 0
798 indj = int(indi + p_wfs._nrebin * p_wfs.npix)
800 X = np.indices((p_wfs._nrebin * p_wfs.npix, p_wfs._nrebin * p_wfs.npix))
801 x = (X[1] / p_wfs._nrebin).astype(np.int64)
802 y = (X[0] / p_wfs._nrebin).astype(np.int64)
805 binindices = np.zeros((p_wfs._Ntot, p_wfs._Ntot))
806 binindices[indi:indj, indi:indj] = x + y * p_wfs.npix + 1
808 binmap = np.zeros((p_wfs._nrebin * p_wfs._nrebin, p_wfs.npix * p_wfs.npix))
810 X = np.indices((p_wfs._Ntot, p_wfs._Ntot))
811 tmp = X[1] + X[0] * p_wfs._Ntot
813 if (p_wfs.gsalt <= 0):
814 binindices = np.roll(binindices, binindices.shape[0] // 2, axis=0)
815 binindices = np.roll(binindices, binindices.shape[1] // 2, axis=1)
817 for i
in range(p_wfs.npix * p_wfs.npix):
818 binmap[:, i] = tmp[np.where(binindices == i + 1)]
820 p_wfs._binmap = np.copy(binmap.astype(np.int32))
822 dr0 = p_tel.diam / r0 * \
823 (0.5 / p_wfs.Lambda) ** 1.2 / \
824 np.cos(p_geom.zenithangle * CONST.DEG2RAD) ** 0.6
825 fwhmseeing = p_wfs.Lambda / \
826 (p_tel.diam / np.sqrt(p_wfs.nxsub ** 2. + (dr0 / 1.5) ** 2.)) / 4.848
827 fwhmseeing = min(fwhmseeing, 2 * p_wfs.pixsize)
828 kernelfwhm = np.sqrt(fwhmseeing**2. + p_wfs.kernel**2.)
830 tmp = util.makegaussian(p_wfs._Ntot, kernelfwhm / p_wfs._qpixsize, p_wfs._Ntot // 2,
831 p_wfs._Ntot // 2).astype(np.float32)
833 tmp = np.roll(tmp, tmp.shape[0] // 2, axis=0)
834 tmp = np.roll(tmp, tmp.shape[1] // 2, axis=1)
837 tmp = tmp / np.sum(tmp)
838 tmp = np.fft.fft2(tmp).astype(np.complex64) / (p_wfs._Ntot * p_wfs._Ntot)
839 p_wfs._ftkernel = np.copy(tmp).astype(np.complex64)
845 if (p_wfs.gsalt == 0):
846 if (p_wfs.zerop == 0):
854 p_wfs._nphotons =
compute_nphotons(scons.WFSType.SH, ittime, p_wfs.optthroughput,
855 p_tel.diam, nxsub=p_wfs.nxsub,
856 zerop=p_wfs.zerop, gsmag=p_wfs.gsmag,
870 p_wfs._nphotons =
compute_nphotons(scons.WFSType.SH, ittime, p_wfs.optthroughput,
871 p_tel.diam, nxsub=p_wfs.nxsub,
872 lgsreturnperwatt=p_wfs.lgsreturnperwatt,
873 laserpower=p_wfs.laserpower, verbose=verbose)
875 if(p_wfs.fssize != 0):
876 fftsize = util.fft_goodsize(p_geom._mpupil.shape[0])
877 fspixsize = (p_geom.pupdiam *
878 (p_wfs.Lambda * 1.e-6) / p_tel.diam * CONST.RAD2ARCSEC) / fftsize
880 fsradius_pixels = int(p_wfs.fssize / fspixsize / 2.)
881 if (p_wfs.fstop == scons.FieldStopType.ROUND):
882 focmask = util.dist(fftsize, xc=fftsize / 2. - 0.5,
883 yc=fftsize / 2. - 0.5) < (fsradius_pixels)
884 elif (p_wfs.fstop == scons.FieldStopType.SQUARE):
885 X = np.indices((fftsize, fftsize)) + 1
886 x = X[1] - (fftsize + 1.) / 2.
887 y = X[0] - (fftsize + 1.) / 2.
888 focmask = (np.abs(x) <= (fsradius_pixels)) * \
889 (np.abs(y) <= (fsradius_pixels))
891 msg =
"wfs fstop must be FieldStopType.[ROUND|SQUARE]"
892 raise ValueError(msg)
894 pyr_focmask = focmask * 1.0
895 p_wfs._submask = np.fft.fftshift(pyr_focmask)
897 def geom_init(p_geom: conf.Param_geom, p_tel: conf.Param_tel, padding=2):
899 Initialize the system geometry
902 p_geom: (Param_geom) : geometry settings
903 p_tel: (Param_tel) : telescope settings
904 padding: (optional) : padding factor for PYRHR geometry
908 p_geom.ssize = int(2**np.ceil(np.log2(p_geom.pupdiam) + 1))
910 p_geom.cent = p_geom.ssize / 2 - 0.5
912 p_geom._p1 = int(np.ceil(p_geom.cent - p_geom.pupdiam / 2.))
913 p_geom._p2 = int(np.floor(p_geom.cent + p_geom.pupdiam / 2.))
915 p_geom.pupdiam = p_geom._p2 - p_geom._p1 + 1
917 p_geom._n = p_geom.pupdiam + 2 * padding
918 p_geom._n1 = p_geom._p1 - padding
919 p_geom._n2 = p_geom._p2 + padding
921 cent = p_geom.pupdiam / 2. - 0.5
924 p_geom._spupil = mkP.make_pupil(p_geom.pupdiam, p_geom.pupdiam, p_tel, cent,
925 cent).astype(np.float32)
928 p_geom._ipupil = util.pad_array(p_geom._spupil, p_geom.ssize).astype(np.float32)
931 p_geom._mpupil = util.pad_array(p_geom._spupil, p_geom._n).astype(np.float32)
933 if (p_tel.std_piston
or p_tel.std_tt):
934 p_geom._phase_ab_M1 = mkP.make_phase_ab(p_geom.pupdiam, p_geom.pupdiam, p_tel,
935 p_geom._spupil, cent,
936 cent).astype(np.float32)
937 p_geom._phase_ab_M1_m = util.pad_array(p_geom._phase_ab_M1,
938 p_geom._n).astype(np.float32)
942 if p_geom.apodFile is None or p_geom.apodFile == '':
943 apod_filename = shesha_savepath + \
944 "apodizer/SP_HARMONI_I4_C6_N1024.npy"
945 p_geom._apodizer = makeA.make_apodizer(
946 p_geom.pupdiam, p_geom.pupdiam,
947 apod_filename.encode(), 180. / 12.).astype(np.float32)
950 p_geom._apodizer = np.ones(p_geom._spupil.shape, dtype=np.int32)
951 p_geom._pixsize = p_tel.diam / p_geom.pupdiam
952 p_geom.is_init =
True
955 def geom_init_generic(p_geom, pupdiam, t_spiders=0.01, spiders_type="six", xc=0, yc=0,
957 """Initialize the system geometry
960 pupdiam: (long) : linear size of total pupil
962 t_spiders: (float) : secondary supports ratio.
964 spiders_type: (str) : secondary supports type: "four" or "six".
972 cobs: (float) : central obstruction ratio.
976 p_geom.ssize = int(2**np.ceil(np.log2(pupdiam) + 1))
978 p_geom.cent = p_geom.ssize / 2 + 0.5
980 pupdiam = int(pupdiam)
981 p_geom._p1 = int(np.ceil(p_geom.cent - pupdiam / 2.))
982 p_geom._p2 = int(np.floor(p_geom.cent + pupdiam / 2.))
983 p_geom.pupdiam = p_geom._p2 - p_geom._p1 + 1
984 p_geom._n = p_geom.pupdiam + 4
985 p_geom._n1 = p_geom._p1 - 2
986 p_geom._n2 = p_geom._p2 + 2
989 p_geom._spupil = mkP.make_pupil_generic(pupdiam, pupdiam, t_spiders, spiders_type,
993 p_geom._ipupil =
pad_array(p_geom._spupil, p_geom.ssize).astype(np.float32)
996 p_geom._mpupil =
pad_array(p_geom._spupil, p_geom._n).astype(np.float32)
1001 D1 = (N - S[0]) // 2
1002 D2 = (N - S[1]) // 2
1003 padded = np.zeros((N, N))
1004 padded[D1:D1 + S[0], D2:D2 + S[1]] = A
Parameter classes for COMPASS.
Numerical constants for shesha and config enumerations for safe-typing.
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)
def geom_init(conf.Param_geom p_geom, conf.Param_tel p_tel, padding=2)
Initialize the system geometry.
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.
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.
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,...
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,...
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.
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,...
Pupil creation functions.
Basic utilities function.