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
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]
63 telescope: (Telescope): Telescope object
66 if p_wfss
is not None:
68 nsensors = len(p_wfss)
70 any_sh = [o.type
for o
in p_wfss].count(scons.WFSType.SH) > 0
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))
80 indmax = np.argsort([o.nxsub
for o
in p_wfss
81 if o.type == scons.WFSType.SH])[-1]
83 indmax = np.argsort([o.nxsub
for o
in p_wfss])[-1]
85 print(
"*-----------------------")
86 print(
"Computing geometry of WFS", indmax)
88 init_wfs_geom(p_wfss[indmax], r0, p_tel, p_geom, ittime, verbose=1)
90 for i
in range(nsensors):
92 print(
"*-----------------------")
93 print(
"Computing geometry of WFS", i)
94 init_wfs_geom(p_wfss[i], r0, p_tel, p_geom, ittime, verbose=1)
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)
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)
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...
117 p_wfs: (Param_wfs) : wfs settings
119 r0: (float) : atmos r0 @ 0.5 microns
121 p_tel: (Param_tel) : telescope settings
123 geom: (Param_geom) : geom settings
125 ittime: (float) : 1/loop frequency [s]
127 verbose: (int) : (optional) display informations if 0
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):
145 if not p_geom.is_init:
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):
158 raise RuntimeError(
"This WFS can not be used")
160 if (p_wfs.type == scons.WFSType.PYRHR
or p_wfs.type == scons.WFSType.PYRLR):
162 elif (p_wfs.type == scons.WFSType.SH):
163 init_sh_geom(p_wfs, r0, p_tel, p_geom, ittime, verbose=1)
165 raise RuntimeError(
"This WFS can not be used")
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)
172 p_wfs: (Param_wfs) : wfs settings
174 r0: (float) : atmos r0 @ 0.5 microns
176 p_tel: (Param_tel) : telescope settings
178 verbose: (int) : (optional) display informations if 0
180 Scheme to determine arrays sizes
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)
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
198 Rf = Fs . N . D / lambda / P
199 ideally we choose : Fs = lambda / D . Nssp / 2
201 if we want good sampling of r0 (avoid aliasing of speckles)
205 to get reasonable space between pupil images : N > P.(2 + 3S)
207 N must be a power of 2
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
215 r0 = r0 * (p_wfs.Lambda * 2)**(6. / 5)
219 print(
"r0 for WFS :",
"%3.2f" % r0,
" m")
222 print(
"seeing for WFS : ",
223 "%3.2f" % (CONST.RAD2ARCSEC * (p_wfs.Lambda * 1.e-6) / r0),
"\"")
225 if (p_wfs._pdiam <= 0):
228 if (p_wfs.type == scons.WFSType.SH):
230 subapdiam = p_tel.diam / float(p_wfs.nxsub)
232 pdiam = int(k * subapdiam / r0)
237 if ((pdiam * p_wfs.nxsub) % 2):
240 nrebin = int(2 * subapdiam * p_wfs.pixsize /
241 (p_wfs.Lambda * 1.e-6) / CONST.RAD2ARCSEC) + 1
242 nrebin = max(2, nrebin)
246 Nfft = util.fft_goodsize(
247 int(pdiam / subapdiam * nrebin / p_wfs.pixsize * CONST.RAD2ARCSEC *
248 (p_wfs.Lambda * 1.e-6)))
250 elif (p_wfs.type == scons.WFSType.PYRHR
or p_wfs.type == scons.WFSType.PYRLR):
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):
256 pdiam = pdiam // p_wfs.nxsub
264 Nfft = util.fft_goodsize(2 * pdiam)
271 if (p_wfs.type == scons.WFSType.SH):
272 subapdiam = p_tel.diam / float(p_wfs.nxsub)
277 qpixsize = (pdiam * (p_wfs.Lambda * 1.e-6) / subapdiam * CONST.RAD2ARCSEC) / Nfft
280 if (p_wfs.pixsize / qpixsize - int(p_wfs.pixsize / qpixsize) > 0.5):
281 nrebin = int(p_wfs.pixsize / qpixsize) + 1
283 nrebin = int(p_wfs.pixsize / qpixsize)
286 pixsize = nrebin * qpixsize
288 if (pixsize * p_wfs.npix > qpixsize * Nfft):
289 Ntot = util.fft_goodsize(int(pixsize * p_wfs.npix / qpixsize) + 1)
293 if (Ntot % 2 != Nfft % 2):
296 elif (p_wfs.type == scons.WFSType.PYRHR
or p_wfs.type == scons.WFSType.PYRLR):
297 pdiam = pdiam * p_wfs.nxsub
300 Nfft = int(2**np.ceil(np.log2(m * pdiam)))
302 nrebin = pdiam // p_wfs.nxsub
303 while (Nfft % nrebin != 0):
305 pdiam = nrebin * p_wfs.nxsub
306 Nfft = int(2**np.ceil(np.log2(m * pdiam)))
309 (p_wfs.Lambda * 1.e-6) / p_tel.diam * CONST.RAD2ARCSEC) / Nfft
311 Ntot = Nfft // pdiam * p_wfs.nxsub
312 pixsize = qpixsize * nrebin
313 pdiam = pdiam // p_wfs.nxsub
316 p_wfs.pixsize = pixsize
317 p_wfs._qpixsize = qpixsize
320 p_wfs._nrebin = nrebin
321 p_wfs._subapd = p_tel.diam / p_wfs.nxsub
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)
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)
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
348 wfs_type: (scons.WFSType) : wfs type: SH or PYRHR.
350 ittime: (float) : 1/loop frequency [s].
352 optthroughput: (float) : wfs global throughput.
354 diam: (float) : telescope diameter.
356 cobs: (float) : (optional for SH) telescope central obstruction.
358 nxsub: (int) : (optional for PYRHR) linear number of subaps.
360 zerop: (float) : (optional for LGS) detector zero point expressed in ph/m**2/s in the bandwidth of the WFS.
362 gsmag: (float) : (optional for LGS) magnitude of guide star.
364 lgsreturnperwatt: (float) : (optional for NGS) return per watt factor (high season : 10 ph/cm2/s/W).
366 laserpower: (float) : (optional for NGS) laser power in W.
368 verbose: (bool) : (optional) display informations if True.
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=?)
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):
384 if (laserpower == 0):
387 surface = (diam / nxsub)**2.
391 nphotons = lgsreturnperwatt * laserpower * \
392 optthroughput * (diam / nxsub) ** 2. * 1e4 * ittime
397 print(
"nphotons : ", nphotons)
400 raise RuntimeError(
"WFS unknown")
402 nphotons = zerop * 10. ** (-0.4 * gsmag) * ittime * \
403 optthroughput * surface
406 print(
"nphotons : ", nphotons)
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...
416 p_wfs: (Param_wfs) : wfs settings
418 r0: (float) : atmos r0 @ 0.5 microns
420 p_tel: (Param_tel) : telescope settings
422 geom: (Param_geom) : geom settings
424 ittime: (float) : 1/loop frequency [s]
426 verbose: (bool) : (optional) display informations if True
430 p_wfs.npix = p_wfs._pdiam
431 p_wfs._Ntot = p_geom._n
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
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))
445 msg =
"PYRHR wfs fstop must be FieldStopType.[ROUND|SQUARE]"
446 raise ValueError(msg)
448 pyr_focmask = focmask * 1.0
449 p_wfs._submask = np.fft.fftshift(pyr_focmask)
452 pyrsize = p_wfs._Nfft
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)
461 pup_sep = p_wfs.pyr_pup_sep
464 y = np.tile(np.arange(pyrsize) - pyrsize / 2, (pyrsize, 1))
467 Pangle = pup_sep * nrebin
468 if p_wfs.nPupils == 0:
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)]
477 if p_wfs.nPupils == 4:
478 centers = np.round(centers[[1, 0, 2, 3], :]).astype(np.int32)
480 if p_wfs.pyr_misalignments
is not None:
481 mis = np.asarray(p_wfs.pyr_misalignments) * nrebin
483 mis = np.zeros((p_wfs.nPupils, 2), dtype=np.float32)
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)
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] = \
498 for qIdx
in range(p_wfs.nPupils):
499 quadOnCenter = np.roll(pup, tuple((centers[qIdx] + mis[qIdx]).astype(np.int32)),
501 mskRebin = util.rebin(quadOnCenter.copy(),
502 [pyrsize // nrebin, pyrsize // nrebin])
504 stackedSubap = np.roll(mskRebin >= fracsub,
505 tuple((-centers[qIdx] / nrebin).astype(np.int32)),
507 mskRebTot = (mskRebin >= fracsub).astype(np.int32)
509 stackedSubap += np.roll(mskRebin >= fracsub,
510 tuple((-centers[qIdx] / nrebin).astype(np.int32)),
512 mskRebTot += (mskRebin >= fracsub).astype(np.int32) * (qIdx + 1)
514 validRow, validCol = [], []
517 if p_wfs.nPupils == 4:
518 centers = centers[[2, 1, 3, 0], :]
520 for qIdx
in range(p_wfs.nPupils):
522 np.roll(stackedSubap, tuple((centers[qIdx] / nrebin).astype(np.int32)),
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)
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)
536 if (p_wfs.pyr_pos
is None):
537 pixsize = (np.pi * p_wfs._qpixsize) / (3600 * 180)
541 scale_fact = 2 * np.pi / pyrsize * \
542 (p_wfs.Lambda * 1e-6 / p_tel.diam) / pixsize * p_wfs.pyr_ampl
544 np.sin((np.arange(p_wfs.pyr_npts)) * 2. * np.pi / p_wfs.pyr_npts)
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)
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
556 p_wfs._pyr_cx = cx.copy()
557 p_wfs._pyr_cy = cy.copy()
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)
569 y = np.tile(np.arange(pyrsize) - pyrsize // 2, (pyrsize, 1))
573 sincar = np.fft.fftshift(np.sinc(x) * np.sinc(y))
576 p_wfs._sincar = sincar.astype(np.float32)
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 +
583 p_wfs._isvalid = pupvalid.T.astype(np.int32)
585 validsubsx = np.where(pupvalid)[0].astype(np.int32)
586 validsubsy = np.where(pupvalid)[1].astype(np.int32)
588 istart = np.arange(p_wfs.nxsub + 2) * p_wfs.npix
589 jstart = np.copy(istart)
592 fluxPerSub = np.zeros((p_wfs.nxsub + 2, p_wfs.nxsub + 2), dtype=np.float32)
594 for i
in range(p_wfs.nxsub + 2):
596 for j
in range(p_wfs.nxsub + 2):
598 fluxPerSub[i, j] = np.sum(
599 p_geom._mpupil[indi:indi + p_wfs.npix, indj:indj + p_wfs.npix])
601 fluxPerSub = fluxPerSub / p_wfs.nxsub**2.
603 p_wfs._fluxPerSub = fluxPerSub.copy()
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))
607 tmp = X[1] + X[0] * p_geom._n
609 pyrtmp = np.zeros((p_geom._n, p_geom._n), dtype=np.int32)
611 for i
in range(len(validsubsx)):
612 indi = istart[validsubsy[i]]
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
617 p_wfs._phasemap = phasemap
619 p_wfs._pyr_offsets = pyrtmp
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...
628 p_wfs: (Param_wfs) : wfs settings
630 r0: (float) : atmos r0 @ 0.5 microns
632 p_tel: (Param_tel) : telescope settings
634 geom: (Param_geom) : geom settings
636 ittime: (float) : 1/loop frequency [s]
638 verbose: (bool) : (optional) display informations if True
643 istart = ((np.linspace(0, p_geom.pupdiam, p_wfs.nxsub + 1))[:-1]).astype(np.int64)
646 jstart = np.copy(istart)
649 fluxPerSub = np.zeros((p_wfs.nxsub, p_wfs.nxsub), dtype=np.float32)
651 for i
in range(p_wfs.nxsub):
653 for j
in range(p_wfs.nxsub):
655 fluxPerSub[i, j] = np.sum(
656 p_geom._mpupil[indi:indi + p_wfs._pdiam, indj:indj + p_wfs._pdiam])
659 fluxPerSub = fluxPerSub / p_wfs._pdiam**2.
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
673 phasemap = np.zeros((p_wfs._pdiam * p_wfs._pdiam, p_wfs._nvalid), dtype=np.int32)
675 X = np.indices((p_geom._n, p_geom._n))
676 tmp = X[1] + X[0] * p_geom._n
681 indi = istart[p_wfs._validsubsy[i]]
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
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))
696 if (p_wfs.npix % 2 == 1
and p_wfs._nrebin % 2 == 1):
698 halfxy = np.zeros((p_wfs._pdiam, p_wfs._pdiam), dtype=np.float32)
699 p_wfs._halfxy = halfxy.astype(np.float32)
701 p_wfs._halfxy = halfxy.astype(np.float32)
704 if (p_wfs._Ntot != p_wfs._Nfft):
705 indi = int((p_wfs._Ntot - p_wfs._Nfft) / 2.)
706 indj = int(indi + p_wfs._Nfft)
707 X = np.indices((p_wfs._Nfft, p_wfs._Nfft)) + 1
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,
714 tmp[indi:indj, indi:indj] = np.roll(tmp[indi:indj, indi:indj], p_wfs._Nfft // 2,
717 tmp = np.roll(tmp, p_wfs._Ntot // 2, axis=0)
718 tmp = np.roll(tmp, p_wfs._Ntot // 2, axis=1)
720 tmp = np.where(tmp.flatten())[0]
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))
727 tmp = np.zeros((2, 2))
728 p_wfs._hrmap = np.copy(tmp.astype(np.int32))
731 if (p_wfs._nrebin * p_wfs.npix % 2 != p_wfs._Ntot % 2):
733 indi = int((p_wfs._Ntot - p_wfs._nrebin * p_wfs.npix) / 2.) + 1
735 indi = int((p_wfs._Ntot - p_wfs._nrebin * p_wfs.npix) / 2.) + 0
736 indj = int(indi + p_wfs._nrebin * p_wfs.npix)
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)
743 binindices = np.zeros((p_wfs._Ntot, p_wfs._Ntot))
744 binindices[indi:indj, indi:indj] = x + y * p_wfs.npix + 1
746 binmap = np.zeros((p_wfs._nrebin * p_wfs._nrebin, p_wfs.npix * p_wfs.npix))
748 X = np.indices((p_wfs._Ntot, p_wfs._Ntot))
749 tmp = X[1] + X[0] * p_wfs._Ntot
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)
755 for i
in range(p_wfs.npix * p_wfs.npix):
756 binmap[:, i] = tmp[np.where(binindices == i + 1)]
758 p_wfs._binmap = np.copy(binmap.astype(np.int32))
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.)
767 tmp = util.makegaussian(p_wfs._Ntot, kernelfwhm / p_wfs._qpixsize, p_wfs._Ntot // 2,
768 p_wfs._Ntot // 2).astype(np.float32)
770 tmp = np.roll(tmp, tmp.shape[0] // 2, axis=0)
771 tmp = np.roll(tmp, tmp.shape[1] // 2, axis=1)
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)
782 if (p_wfs.gsalt == 0):
783 if (p_wfs.zerop == 0):
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,
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)
813 def geom_init(p_geom: conf.Param_geom, p_tel: conf.Param_tel, padding=2):
815 Initialize the system geometry
818 p_geom: (Param_geom) : geometry settings
819 p_tel: (Param_tel) : telescope settings
820 padding: (optional) : padding factor for PYRHR geometry
824 p_geom.ssize = int(2**np.ceil(np.log2(p_geom.pupdiam) + 1))
826 p_geom.cent = p_geom.ssize / 2 - 0.5
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.))
831 p_geom.pupdiam = p_geom._p2 - p_geom._p1 + 1
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
837 cent = p_geom.pupdiam / 2. - 0.5
840 p_geom._spupil = mkP.make_pupil(p_geom.pupdiam, p_geom.pupdiam, p_tel, cent,
841 cent).astype(np.float32)
844 p_geom._ipupil = util.pad_array(p_geom._spupil, p_geom.ssize).astype(np.float32)
847 p_geom._mpupil = util.pad_array(p_geom._spupil, p_geom._n).astype(np.float32)
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)
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)
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
871 def geom_init_generic(p_geom, pupdiam, t_spiders=0.01, spiders_type="six", xc=0, yc=0,
873 """Initialize the system geometry
876 pupdiam: (long) : linear size of total pupil
878 t_spiders: (float) : secondary supports ratio.
880 spiders_type: (str) : secondary supports type: "four" or "six".
888 cobs: (float) : central obstruction ratio.
892 p_geom.ssize = int(2**np.ceil(np.log2(pupdiam) + 1))
894 p_geom.cent = p_geom.ssize / 2 + 0.5
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
905 p_geom._spupil = mkP.make_pupil_generic(pupdiam, pupdiam, t_spiders, spiders_type,
909 p_geom._ipupil =
pad_array(p_geom._spupil, p_geom.ssize).astype(np.float32)
912 p_geom._mpupil =
pad_array(p_geom._spupil, p_geom._n).astype(np.float32)
919 padded = np.zeros((N, N))
920 padded[D1:D1 + S[0], D2:D2 + S[1]] = A