43 from shesha.util import dm_util, influ_util, kl_util
49 from scipy
import interpolate
52 from typing
import List
54 from rich.progress
import track
58 shesha_dm = os.environ[
'SHESHA_DM_ROOT']
59 except KeyError
as err:
60 shesha_dm = os.environ[
'SHESHA_ROOT'] +
"/data/dm-data"
63 def dm_init(context: carmaWrap_context, p_dms: List[conf.Param_dm],
64 p_tel: conf.Param_tel, p_geom: conf.Param_geom,
65 p_wfss: List[conf.Param_wfs] =
None) -> Dms:
66 """Create and initialize a Dms object on the gpu
69 context: (carmaWrap_context): context
70 p_dms: (list of Param_dms) : dms settings
71 p_tel: (Param_tel) : telescope settings
72 p_geom: (Param_geom) : geom settings
73 p_wfss: (list of Param_wfs) : wfs settings
75 Dms: (Dms): Dms object
78 if (p_wfss
is not None):
81 for i
in range(len(p_wfss)):
82 xpos_wfs.append(p_wfss[i].xpos)
83 ypos_wfs.append(p_wfss[i].ypos)
89 types_dm = [p_dm.type
for p_dm
in p_dms]
90 if scons.DmType.TT
in types_dm:
91 first_TT = types_dm.index(scons.DmType.TT)
92 if np.any(np.array(types_dm[first_TT:]) != scons.DmType.TT):
93 raise RuntimeError(
"TT must be defined at the end of the dms parameters")
95 for i
in range(len(p_dms)):
96 max_extent = _dm_init(context, dms, p_dms[i], xpos_wfs, ypos_wfs, p_geom,
97 p_tel.diam, p_tel.cobs, p_tel.pupangle, max_extent)
102 def _dm_init(context: carmaWrap_context, dms: Dms, p_dm: conf.Param_dm, xpos_wfs: list,
103 ypos_wfs: list, p_geom: conf.Param_geom, diam: float, cobs: float,
104 pupAngle: float, max_extent: int):
105 """ inits a Dms object on the gpu
108 context: (carmaWrap_context): context
109 dms: (Dms) : dm object
111 p_dm: (Param_dms) : dm settings
113 xpos_wfs: (list) : list of wfs xpos
115 ypos_wfs: (list) : list of wfs ypos
117 p_geom: (Param_geom) : geom settings
119 diam: (float) : diameter of telescope
121 cobs: (float) : cobs of telescope
123 pupAngle: (float) : rotation/clocking angle of the pupil in degrees
125 max_extent: (int) : maximum dimension of all dms
128 max_extent: (int) : new maximum dimension of all dms
132 if (p_dm.pupoffset
is not None):
133 p_dm._puppixoffset = p_dm.pupoffset / diam * p_geom.pupdiam
135 patchDiam = dm_util.dim_dm_patch(p_geom.pupdiam, diam, p_dm.type, p_dm.alt, xpos_wfs,
138 if (p_dm.type == scons.DmType.PZT):
139 if p_dm.file_influ_fits ==
None:
140 if p_dm._pitch
is None:
141 p_dm._pitch = patchDiam / float(p_dm.nact - 1)
142 print(f
"DM pitch = {p_dm._pitch:8.5f} pix = {p_dm._pitch*diam/p_geom.pupdiam:8.5f} m",
145 extent = p_dm._pitch * (p_dm.nact + p_dm.pzt_extent)
146 p_dm._n1, p_dm._n2 = dm_util.dim_dm_support(p_geom.cent, extent,
155 max_extent = max(max_extent, p_dm._n2 - p_dm._n1 + 1)
157 dim = max(p_dm._n2 - p_dm._n1 + 1, p_geom._mpupil.shape[0])
158 p_dm._dim_screen = dim
159 ninflupos = p_dm._influpos.size
160 n_npts = p_dm._ninflu.size
161 dms.add_dm(context, p_dm.type, p_dm.alt, dim, p_dm._ntotact, p_dm._influsize,
162 ninflupos, n_npts, p_dm.push4imat, 0, context.active_device)
164 dms.d_dms[-1].pzt_loadarrays(p_dm._influ, p_dm._influpos.astype(np.int32),
165 p_dm._ninflu, p_dm._influstart, p_dm._i1, p_dm._j1)
167 elif (p_dm.type == scons.DmType.TT):
169 if (p_dm.alt == 0)
and (max_extent != 0):
170 extent = int(max_extent * 1.05)
171 if (extent % 2 != 0):
174 extent = p_geom.pupdiam + 16
175 p_dm._n1, p_dm._n2 = dm_util.dim_dm_support(p_geom.cent, extent, p_geom.ssize)
177 max_extent = max(max_extent, p_dm._n2 - p_dm._n1 + 1)
179 dim = p_dm._n2 - p_dm._n1 + 1
181 p_dm._dim_screen = dim
182 dms.add_dm(context, p_dm.type, p_dm.alt, dim, 2, dim, 1, 1, p_dm.push4imat, 0,
183 context.active_device)
184 dms.d_dms[-1].tt_loadarrays(p_dm._influ)
186 elif (p_dm.type == scons.DmType.KL):
188 extent = p_geom.pupdiam + 16
189 p_dm._n1, p_dm._n2 = dm_util.dim_dm_support(p_geom.cent, extent, p_geom.ssize)
191 max_extent = max(max_extent, p_dm._n2 - p_dm._n1 + 1)
193 dim = p_dm._n2 - p_dm._n1 + 1
198 p_dm._dim_screen = dim
200 dms.add_dm(context, p_dm.type, p_dm.alt, dim, p_dm.nkl, p_dm._ncp, p_dm._nr,
201 p_dm._npp, p_dm.push4imat, p_dm._ord.max(), context.active_device)
203 dms.d_dms[-1].kl_loadarrays(p_dm._rabas, p_dm._azbas, p_dm._ord, p_dm._cr,
208 raise TypeError(
"This type of DM doesn't exist ")
216 def _dm_init_factorized(context: carmaWrap_context, dms: Dms, p_dm: conf.Param_dm,
217 xpos_wfs: list, ypos_wfs: list, p_geom: conf.Param_geom,
218 diam: float, cobs: float, pupAngle: float, max_extent: int):
219 """ inits a Dms object on the gpu
223 context: (carmaWrap_context): context
225 dms: (Dms) : dm object
227 p_dm: (Param_dms) : dm settings
229 xpos_wfs: (list) : list of wfs xpos
231 ypos_wfs: (list) : list of wfs ypos
233 p_geom: (Param_geom) : geom settings
235 diam: (float) : diameter of telescope
237 cobs: (float) : cobs of telescope
239 pupAngle: (float) : rotation/clocking angle of the pupil in degrees
241 max_extent: (int) : maximum dimension of all dms
244 max_extent: (int) : new maximum dimension of all dms
248 if (p_dm.pupoffset
is not None):
249 p_dm._puppixoffset = p_dm.pupoffset / diam * p_geom.pupdiam
251 patchDiam = dm_util.dim_dm_patch(p_geom.pupdiam, diam, p_dm.type, p_dm.alt, xpos_wfs,
254 if (p_dm.type == scons.DmType.PZT)
and p_dm.file_influ_fits
is not None:
257 if (p_dm.type == scons.DmType.PZT):
258 p_dm._pitch = patchDiam / float(p_dm.nact - 1)
260 extent = p_dm._pitch * (p_dm.nact + p_dm.pzt_extent)
265 elif (p_dm.type == scons.DmType.TT):
266 if (p_dm.alt == 0)
and (max_extent != 0):
267 extent = int(max_extent * 1.05)
268 if (extent % 2 != 0):
271 extent = p_geom.pupdiam + 16
273 elif (p_dm.type == scons.DmType.KL):
274 extent = p_geom.pupdiam + 16
276 raise TypeError(
"This type of DM doesn't exist ")
282 p_dm._n1, p_dm._n2 = dm_util.dim_dm_support(p_geom.cent, extent, p_geom.ssize)
285 max_extent = max(max_extent, p_dm._n2 - p_dm._n1 + 1)
287 dim = max(p_dm._n2 - p_dm._n1 + 1, p_geom._mpupil.shape[0])
289 if (p_dm.type == scons.DmType.PZT):
290 ninflupos = p_dm._influpos.size
291 n_npts = p_dm._ninflu.size
292 dms.add_dm(context, p_dm.type, p_dm.alt, dim, p_dm._ntotact, p_dm._influsize,
293 ninflupos, n_npts, p_dm.push4imat, 0, context.active_device)
295 dms.d_dms[-1].pzt_loadarrays(p_dm._influ, p_dm._influpos.astype(np.int32),
296 p_dm._ninflu, p_dm._influstart, p_dm._i1, p_dm._j1)
297 elif (p_dm.type == scons.DmType.TT):
299 dms.add_dm(context, p_dm.type, p_dm.alt, dim, 2, dim, 1, 1, p_dm.push4imat, 0,
300 context.active_device)
301 dms.d_dms[-1].tt_loadarrays(p_dm._influ)
302 elif (p_dm.type == scons.DmType.KL):
306 dms.add_dm(context, p_dm.type, p_dm.alt, dim, p_dm.nkl, p_dm._ncp, p_dm._nr,
307 p_dm._npp, p_dm.push4imat, p_dm._ord.max(), context.active_device)
308 dms.d_dms[-1].kl_loadarrays(p_dm._rabas, p_dm._azbas, p_dm._ord, p_dm._cr,
314 def dm_init_standalone(context: carmaWrap_context, p_dms: list, p_geom: conf.Param_geom,
315 diam=1., cobs=0., pupAngle=0., wfs_xpos=[0], wfs_ypos=[0]):
316 """Create and initialize a Dms object on the gpu
319 p_dms: (list of Param_dms) : dms settings
321 p_geom: (Param_geom) : geom settings
323 diam: (float) : diameter of telescope (default 1.)
325 cobs: (float) : cobs of telescope (default 0.)
327 pupAngle: (float) : pupil rotation angle (degrees, default 0.)
329 wfs_xpos: (array) : guide star x position on sky (in arcsec).
331 wfs_ypos: (array) : guide star y position on sky (in arcsec).
335 if (len(p_dms) != 0):
337 for i
in range(len(p_dms)):
338 _dm_init(context, dms, p_dms[i], wfs_xpos, wfs_ypos, p_geom, diam, cobs,
339 pupAngle, max_extent)
343 def make_pzt_dm(p_dm: conf.Param_dm, p_geom: conf.Param_geom, cobs: float,
345 """Compute the actuators positions and the influence functions for a pzt DM.
346 NOTE: if the DM is in altitude, central obstruction is forced to 0
349 p_dm: (Param_dm) : dm parameters
351 p_geom: (Param_geom) : geometry parameters
353 cobs: (float) : telescope central obstruction
355 pupAngle: (float) : rotation/clocking angle of the pupil in degrees
358 influ: (np.ndarray(dims=3, dtype=np.float64)) : cube of the IF for each actuator
362 if (p_dm.influ_type == scons.InfluType.PETAL):
367 coupling = p_dm.coupling
371 if (p_dm.influ_type == scons.InfluType.RADIALSCHWARTZ):
372 smallsize = influ_util.makeRadialSchwartz(pitch, coupling)
373 elif (p_dm.influ_type == scons.InfluType.SQUARESCHWARTZ):
374 smallsize = influ_util.makeSquareSchwartz(pitch, coupling)
375 elif (p_dm.influ_type == scons.InfluType.BLACKNUTT):
376 smallsize = influ_util.makeBlacknutt(pitch, coupling)
377 elif (p_dm.influ_type == scons.InfluType.GAUSSIAN):
378 smallsize = influ_util.makeGaussian(pitch, coupling)
379 elif (p_dm.influ_type == scons.InfluType.BESSEL):
380 smallsize = influ_util.makeBessel(pitch, coupling, p_dm.type_pattern)
381 elif (p_dm.influ_type == scons.InfluType.DEFAULT):
382 smallsize = influ_util.makeRigaut(pitch, coupling)
384 print(
"ERROR influtype not recognized ")
385 p_dm._influsize = smallsize
388 if p_dm.type_pattern
is None:
389 p_dm.type_pattern = scons.PatternType.SQUARE
391 if p_dm.type_pattern == scons.PatternType.HEXA:
392 print(
"Pattern type : hexa")
393 xypos = dm_util.createHexaPattern(pitch, p_geom.pupdiam * 1.1)
394 p_dm.keep_all_actu =
True
395 elif p_dm.type_pattern == scons.PatternType.HEXAM4:
396 print(
"Pattern type : hexaM4")
397 p_dm.keep_all_actu =
True
398 xypos = dm_util.createDoubleHexaPattern(pitch, p_geom.pupdiam * 1.1, pupAngle)
399 if p_dm.margin_out
is not None:
400 print(f
'p_dm.margin_out={p_dm.margin_out} is being '
401 'used for pupil-based actuator filtering')
402 pup_side = p_geom._ipupil.shape[0]
403 cub_off = dm_util.filterActuWithPupil(xypos + pup_side // 2 - 0.5,
405 p_dm.margin_out * p_dm.get_pitch())
406 xypos = cub_off - pup_side // 2 + 0.5
407 p_dm.set_ntotact(xypos.shape[1])
408 elif p_dm.type_pattern == scons.PatternType.SQUARE:
409 print(
"Pattern type : square")
410 xypos = dm_util.createSquarePattern(pitch, p_dm.nact + 4)
412 raise ValueError(
"This pattern does not exist for pzt dm")
414 if p_dm.keep_all_actu:
415 inbigcirc = np.arange(xypos.shape[1])
419 inbigcirc = dm_util.select_actuators(xypos[0, :], xypos[1, :], p_dm.nact,
420 p_dm._pitch, cobs, p_dm.margin_in,
421 p_dm.margin_out, p_dm._ntotact)
422 p_dm._ntotact = inbigcirc.size
428 xypos = xypos[:, inbigcirc]
429 ntotact = xypos.shape[1]
432 i1t = (xpos - smallsize / 2 - 0.5 - p_dm._n1).astype(np.int32)
433 j1t = (ypos - smallsize / 2 - 0.5 - p_dm._n1).astype(np.int32)
441 influ = np.zeros((smallsize, smallsize, ntotact), dtype=np.float32)
444 print(
"Computing Influence Function type : ", p_dm.influ_type)
445 for i
in track(range(ntotact)):
448 x = np.tile(np.arange(i1, i1 + smallsize, dtype=np.float32),
456 y = np.tile(np.arange(j1, j1 + smallsize, dtype=np.float32),
462 if (p_dm.influ_type == scons.InfluType.RADIALSCHWARTZ):
463 influ[:, :, i] = influ_util.makeRadialSchwartz(pitch, coupling, x=x, y=y)
464 elif (p_dm.influ_type == scons.InfluType.SQUARESCHWARTZ):
465 influ[:, :, i] = influ_util.makeSquareSchwartz(pitch, coupling, x=x, y=y)
466 elif (p_dm.influ_type == scons.InfluType.BLACKNUTT):
467 influ[:, :, i] = influ_util.makeBlacknutt(pitch, coupling, x=x, y=y)
468 elif (p_dm.influ_type == scons.InfluType.GAUSSIAN):
469 influ[:, :, i] = influ_util.makeGaussian(pitch, coupling, x=x, y=y)
470 elif (p_dm.influ_type == scons.InfluType.BESSEL):
471 influ[:, :, i] = influ_util.makeBessel(pitch, coupling, x=x, y=y,
472 patternType=p_dm.type_pattern)
473 elif (p_dm.influ_type == scons.InfluType.DEFAULT):
474 influ[:, :, i] = influ_util.makeRigaut(pitch, coupling, x=x, y=y)
476 print(
"ERROR influtype not recognized (defaut or gaussian or bessel)")
478 if (p_dm._puppixoffset
is not None):
479 xpos += p_dm._puppixoffset[0]
480 ypos += p_dm._puppixoffset[1]
482 if p_dm.segmented_mirror:
485 s = (p_geom._ipupil.shape[0] - p_geom._mpupil.shape[0]) // 2
487 from skimage.morphology
import label
489 for i
in track(range(ntotact)):
491 i1, j1 = i1t[i] + s - smallsize // 2, j1t[i] + s - smallsize // 2
492 pupilSnapshot = p_geom._ipupil[i1:i1 + smallsize, j1:j1 + smallsize]
493 if np.all(pupilSnapshot):
495 labels, num = label(pupilSnapshot, background=0, return_num=
True)
499 maxPerArea = np.array([
500 (influ[:, :, i] * (labels == k).astype(np.float32)).max()
501 for k
in range(1, num + 1)
503 influ[:, :, i] *= (labels == np.argmax(maxPerArea) + 1).astype(np.float32)
504 print(f
'{k} cross-spider influence functions trimmed.')
506 influ = influ * float(p_dm.unitpervolt / np.max(influ))
512 dim = max(p_geom._mpupil.shape[0], p_dm._n2 - p_dm._n1 + 1)
513 off = (dim - p_dm._influsize) // 2
517 def init_custom_dm(p_dm: conf.Param_dm, p_geom: conf.Param_geom, diam: float):
518 """Read Fits for influence pzt fonction and form
521 p_dm: (Param_dm) : dm settings
523 p_geom: (Param_geom) : geom settings
525 diam: (float) : tel diameter
527 Conversion. There are sereval coordinate systems.
528 Some are coming from the input fits file, others from compass.
529 Those systems differ by scales and offsets.
530 Coord systems from the input fits file:
531 - they all have the same scale: the coordinates are expressed in
533 - one system is the single common frame where fits data are described [i]
534 - one is local to the minimap [l]
535 Coord systems from compass:
536 - they all have the same scale (different from fits one) expressed in
538 - one system is attached to ipupil (i=image, largest support) [i]
539 - one system is attached to mpupil (m=medium, medium support) [m]
540 - one system is local to minimap [l)]
542 Variables will be named using
543 f, c: define either fits or compass
546 from astropy.io
import fits
as pfits
549 file_name = p_dm.file_influ_fits
550 if(
not os.path.isfile(file_name)):
551 file_name = shesha_dm +
"/" + p_dm.file_influ_fits
553 hdul = pfits.open(file_name)
554 print(
"Read influence function from fits file : ", file_name)
556 dm_fits_version = hdul[0].header[
'VERSION']
559 f_xC = hdul[0].header[
'XCENTER']
560 f_yC = hdul[0].header[
'YCENTER']
561 f_pixsize = hdul[0].header[
'PIXSIZE']
562 if(dm_fits_version < 1.2 ):
563 fi_i1 , fi_j1 = hdul[1].data
564 f_influ = hdul[2].data
565 f_xpos, f_ypos = hdul[3].data
567 fi_i1, fi_j1 = hdul[
'I1_J1'].data
568 f_influ = hdul[
'INFLU'].data
569 f_xpos, f_ypos = hdul[
'XPOS_YPOS'].data
572 cases = [ p_dm.diam_dm
is not None, p_dm._pitch
is not None,
573 p_dm.diam_dm_proj
is not None ]
576 if cases == [
False,
False,
False]:
577 f_pupm = hdul[0].header[
'PUPM']
578 scale = diam / f_pupm
579 print(
'Custom DM: stretching DM to fit PUPM (%f) to compass (%f)' % (f_pupm, diam))
580 elif cases == [
True,
False,
False]:
581 scale = diam / p_dm.diam_dm
582 elif cases == [
False,
True,
False]:
583 f_pitchm = hdul[0].header[
'PITCHM']
584 scale = p_dm._pitch / f_pitchm
585 elif cases == [
False,
False,
True]:
586 f_pupm = hdul[0].header[
'PUPM']
587 scale = p_dm.diam_dm_proj / f_pupm
589 err_msg =
'''Invalid rescaling parameters
590 To set the scale of the custom dm, MAXIMUM ONE of the following parameters should be set:
591 - p_dm.set_pitch(val) : if the pitch in the tel pupil plane is known
592 - p_dm.set_diam_dm(val) : if the pupil diameter in the dm plane is known
593 - p_dm.set_diam_dm_proj(val) : if the dm pupil diameter projected in the tel pupil plane is known
594 If none of the above is set, the dm will be scaled so that the PUPM parameter in the fits file matches the tel pupil.
596 raise RuntimeError(err_msg)
599 print(
"Custom dm scaling factor to pupil plane :", scale)
604 scaleToCompass = np.float32(f_pixsize) / np.float32(p_geom._pixsize)
609 offsetXToCompass = p_geom.cent - f_xC * scaleToCompass
610 offsetYToCompass = p_geom.cent - f_yC * scaleToCompass
613 iSize, jSize, ntotact = f_influ.shape
617 ess1 = np.ceil((fi_i1+iSize) * scaleToCompass + offsetXToCompass) \
618 - np.floor(fi_i1 * scaleToCompass + offsetXToCompass)
621 ess2 = np.ceil((fi_j1+jSize) * scaleToCompass + offsetYToCompass) \
622 - np.floor(fi_j1 * scaleToCompass + offsetYToCompass)
624 smallsize = np.maximum(ess1, ess2).astype(int)
627 p_dm._ntotact = ntotact
628 p_dm._influsize = np.int64(smallsize)
629 p_dm._i1 = np.zeros(ntotact, dtype=np.int32)
630 p_dm._j1 = np.zeros(ntotact, dtype=np.int32)
631 p_dm._xpos = np.zeros(ntotact, dtype=np.float32)
632 p_dm._ypos = np.zeros(ntotact, dtype=np.float32)
633 p_dm._influ = np.zeros((smallsize, smallsize, ntotact), dtype=np.float32)
636 for i
in range(ntotact):
638 fi_x = np.arange(iSize) + fi_i1[i]
639 fi_y = np.arange(jSize) + fi_j1[i]
642 ci_x = fi_x * scaleToCompass + offsetXToCompass
643 ci_y = fi_y * scaleToCompass + offsetYToCompass
646 ci_i1 = fi_i1[i] * scaleToCompass + offsetXToCompass
647 ci_j1 = fi_j1[i] * scaleToCompass + offsetYToCompass
648 ci_i1 = np.floor(ci_i1).astype(np.int32)
649 ci_j1 = np.floor(ci_j1).astype(np.int32)
650 ci_xpix = ci_i1 + np.arange(smallsize)
651 ci_ypix = ci_j1 + np.arange(smallsize)
656 f = interpolate.interp2d(ci_y, ci_x, f_influ[:, :, i], kind=
'cubic')
657 temp = f(ci_ypix, ci_xpix) * p_dm.unitpervolt
659 p_dm._influ[:, :, i] = temp
665 tmp = p_geom.ssize - (np.max(p_dm._i1 + smallsize))
666 margin_i = np.minimum(tmp, np.min(p_dm._i1))
667 tmp = p_geom.ssize - (np.max(p_dm._j1 + smallsize))
668 margin_j = np.minimum(tmp, np.min(p_dm._j1))
670 p_dm._xpos = f_xpos * scaleToCompass + offsetXToCompass
671 p_dm._ypos = f_ypos * scaleToCompass + offsetYToCompass
673 p_dm._n1 = int(np.minimum(margin_i, margin_j))
674 p_dm._n2 = p_geom.ssize - 1 - p_dm._n1
682 def make_tiptilt_dm(p_dm: conf.Param_dm, patchDiam: int, p_geom: conf.Param_geom,
684 """Compute the influence functions for a tip-tilt DM
687 p_dm: (Param_dm) : dm settings
689 patchDiam: (int) : patchDiam for dm size
691 p_geom: (Param_geom) : geom settings
693 diam: (float) : telescope diameter
695 influ: (np.ndarray(dims=3,dtype=np.float64)) : cube of the IF
698 dim = max(p_dm._n2 - p_dm._n1 + 1, p_geom._mpupil.shape[0])
702 influ = dm_util.make_zernike(nzer + 1, dim, patchDiam, p_geom.cent - p_dm._n1 - 0.5,
703 p_geom.cent - p_dm._n1 - 0.5, 1)[:, :, 1:]
706 current = influ[dim // 2 - 1, dim // 2 - 1, 0] - \
707 influ[dim // 2 - 2, dim // 2 - 2, 0]
708 fact = p_dm.unitpervolt * diam / p_geom.pupdiam * 4.848 / current
711 p_dm._ntotact = influ.shape[2]
712 p_dm._influsize = influ.shape[0]
718 def make_kl_dm(p_dm: conf.Param_dm, patchDiam: int, p_geom: conf.Param_geom,
719 cobs: float) ->
None:
720 """Compute the influence function for a Karhunen-Loeve DM
723 p_dm: (Param_dm) : dm settings
725 patchDiam: (int) : patchDiam for dm size
727 p_geom: (Param_geom) : geom settings
729 cobs: (float) : telescope cobs
732 dim = p_geom._mpupil.shape[0]
734 print(
"KL type: ", p_dm.type_kl)
737 nr = np.int64(5.0 * np.sqrt(52))
738 npp = np.int64(10.0 * nr)
740 nr = np.int64(5.0 * np.sqrt(p_dm.nkl))
741 npp = np.int64(10.0 * nr)
743 radp = kl_util.make_radii(cobs, nr)
745 kers = kl_util.make_kernels(cobs, nr, radp, p_dm.type_kl, p_dm.outscl)
747 evals, nord, npo, ordd, rabas = kl_util.gkl_fcom(kers, cobs, p_dm.nkl)
749 azbas = kl_util.make_azimuth(nord, npp)
751 ncp, ncmar, px, py, cr, cp, pincx, pincy, pincw, ap = kl_util.set_pctr(
752 patchDiam, nr, npp, p_dm.nkl, cobs, nord)
754 p_dm._ntotact = p_dm.nkl
763 p_dm._i1 = np.zeros((p_dm.nkl), dtype=np.int32) + \
764 (dim - patchDiam) // 2
765 p_dm._j1 = np.zeros((p_dm.nkl), dtype=np.int32) + \
766 (dim - patchDiam) // 2
767 p_dm._ntotact = p_dm.nkl
771 def comp_dmgeom(p_dm: conf.Param_dm, p_geom: conf.Param_geom):
772 """Compute the geometry of a DM : positions of actuators and influence functions
775 dm: (Param_dm) : dm settings
777 geom: (Param_geom) : geom settings
779 smallsize = p_dm._influsize
781 dm_dim = int(p_dm._n2 - p_dm._n1 + 1)
782 mpup_dim = p_geom._mpupil.shape[0]
784 if (dm_dim < mpup_dim):
785 print(
'DM support is smaller than mpupil')
786 offs = (mpup_dim - dm_dim) // 2
788 print(
'DM support is larger than mpupil')
792 indgen = np.tile(np.arange(smallsize, dtype=np.int32), (smallsize, 1))
794 tmpx = np.tile(indgen, (nact, 1, 1))
795 tmpy = np.tile(indgen.T, (nact, 1, 1))
797 tmpx += offs + p_dm._i1[:,
None,
None]
798 tmpy += offs + p_dm._j1[:,
None,
None]
800 tmp = tmpx + mpup_dim * tmpy
803 tmp[tmpx < 0] = mpup_dim * mpup_dim + 10
804 tmp[tmpy < 0] = mpup_dim * mpup_dim + 10
805 tmp[tmpx > dm_dim - 1] = mpup_dim * mpup_dim + 10
806 tmp[tmpy > dm_dim - 1] = mpup_dim * mpup_dim + 10
807 itmps = np.argsort(tmp.flatten()).astype(np.int32)
808 tmps = tmp.flatten()[itmps].astype(np.int32)
809 itmps = itmps[np.where(itmps > -1)]
811 istart = np.zeros((mpup_dim * mpup_dim), dtype=np.int32)
812 npts = np.zeros((mpup_dim * mpup_dim), dtype=np.int32)
814 tmps_unique, cpt = np.unique(tmps, return_counts=
True)
815 if (tmps_unique > npts.size - 1).any():
816 tmps_unique = tmps_unique[:-1]
819 for i
in range(tmps_unique.size):
820 npts[tmps_unique[i]] = cpt[i]
821 istart[1:] = np.cumsum(npts[:-1])
823 p_dm._influpos = itmps[:np.sum(npts)].astype(np.int32)
827 p_dm._ninflu = npts.astype(np.int32)
828 p_dm._influstart = istart.astype(np.int32)
840 def correct_dm(context, dms: Dms, p_dms: list, p_controller: conf.Param_controller,
841 p_geom: conf.Param_geom, imat: np.ndarray =
None, dataBase: dict = {},
842 use_DB: bool =
False):
843 """Correct the geometry of the DMs using the imat (filter unseen actuators)
846 context: (carmaWrap_context): context
847 dms: (Dms) : Dms object
848 p_dms: (list of Param_dm) : dms settings
849 p_controller: (Param_controller) : controller settings
850 p_geom: (Param_geom) : geom settings
851 imat: (np.ndarray) : interaction matrix
852 dataBase: (dict): dictionary containing paths to files to load
853 use_DB: (bool): dataBase use flag
855 print(
"Filtering unseen actuators... ")
857 resp = np.sqrt(np.sum(imat**2, axis=0))
859 ndm = p_controller.ndm.size
862 for nmc
in range(ndm):
863 nm = p_controller.ndm[nmc]
864 nactu_nm = p_dms[nm]._ntotact
866 if (p_dms[nm].type == scons.DmType.PZT):
868 influpos, ninflu, influstart, i1, j1, ok = h5u.load_dm_geom_from_dataBase(
870 p_dms[nm].set_ntotact(ok.shape[0])
871 p_dms[nm].set_influ(p_dms[nm]._influ[:, :, ok.tolist()])
872 p_dms[nm].set_xpos(p_dms[nm]._xpos[ok])
873 p_dms[nm].set_ypos(p_dms[nm]._ypos[ok])
874 p_dms[nm]._influpos = influpos
875 p_dms[nm]._ninflu = ninflu
876 p_dms[nm]._influstart = influstart
880 tmp = resp[inds:inds + p_dms[nm]._ntotact]
881 ok = np.where(tmp > p_dms[nm].thresh * np.max(tmp))[0]
882 nok = np.where(tmp <= p_dms[nm].thresh * np.max(tmp))[0]
884 p_dms[nm].set_ntotact(ok.shape[0])
885 p_dms[nm].set_influ(p_dms[nm]._influ[:, :, ok.tolist()])
886 p_dms[nm].set_xpos(p_dms[nm]._xpos[ok])
887 p_dms[nm].set_ypos(p_dms[nm]._ypos[ok])
888 p_dms[nm].set_i1(p_dms[nm]._i1[ok])
889 p_dms[nm].set_j1(p_dms[nm]._j1[ok])
893 h5u.save_dm_geom_in_dataBase(nmc, p_dms[nm]._influpos,
895 p_dms[nm]._influstart, p_dms[nm]._i1,
898 dim = max(p_dms[nm]._n2 - p_dms[nm]._n1 + 1, p_geom._mpupil.shape[0])
899 ninflupos = p_dms[nm]._influpos.size
900 n_npts = p_dms[nm]._ninflu.size
902 dms.insert_dm(context, p_dms[nm].type, p_dms[nm].alt, dim,
903 p_dms[nm]._ntotact, p_dms[nm]._influsize, ninflupos, n_npts,
904 p_dms[nm].push4imat, 0, p_dms[nm].dx / p_geom._pixsize,
905 p_dms[nm].dy / p_geom._pixsize, p_dms[nm].theta, p_dms[nm].G,
906 context.active_device, nm)
907 dms.d_dms[nm].pzt_loadarrays(p_dms[nm]._influ, p_dms[nm]._influpos.astype(
908 np.int32), p_dms[nm]._ninflu, p_dms[nm]._influstart, p_dms[nm]._i1,
917 makePetalDm(p_dm, p_geom, pupAngleDegree)
919 The function builds a DM, segmented in petals according to the pupil
920 shape. The petals will be adapted to the EELT case only.
922 <p_geom> : compass object p_geom. The function requires the object p_geom
923 in order to know what is the pupil mask, and what is the mpupil.
924 <p_dm> : compass petal dm object p_dm to be created. The function will
925 transform/modify in place the attributes of the object p_dm.
926 <pupAngleDegree> : rotation/clocking angle of the pupil in degrees
930 p_dm._n1 = p_geom._n1
931 p_dm._n2 = p_geom._n2
933 p_dm._influsize = smallsize
934 p_dm.set_ntotact(nbSeg)
937 p_dm._xpos = i1 + smallsize / 2 + p_dm._n1
938 p_dm._ypos = j1 + smallsize / 2 + p_dm._n1
947 <pupImage> : image of the pupil
948 <pupAngleDegree> : rotation angle of the pupil in degrees
950 La fonction renvoie des fn d'influence en forme de petale d'apres
951 une image de la pupille, qui est supposee etre segmentee.
954 influ, i1, j1, smallsize, nbSeg = make_petal_dm_core(pupImage, 0.0)
960 from scipy.ndimage
import label
961 from scipy.ndimage.morphology
import binary_opening
962 s = np.ones((2, 2), dtype=bool)
963 segments, nbSeg = label(binary_opening(pupImage, s))
972 for i
in range(nbSeg):
973 petal = segments == (i + 1)
974 profil = np.sum(petal, axis=1) != 0
975 extent = np.sum(profil).astype(np.int32)
976 i1t.append(np.min(np.where(profil)[0]))
977 i2t.append(np.max(np.where(profil)[0]))
978 if extent > smallsize:
981 profil = np.sum(petal, axis=0) != 0
982 extent = np.sum(profil).astype(np.int32)
983 j1t.append(np.min(np.where(profil)[0]))
984 j2t.append(np.max(np.where(profil)[0]))
985 if extent > smallsize:
992 influ = np.zeros((smallsize, smallsize, nbSeg), dtype=np.float32)
994 npt = pupImage.shape[0]
995 i0 = j0 = npt / 2 - 0.5
996 petalMap =
build_petals(nbSeg, pupAngleDegree, i0, j0, npt)
997 ii1 = np.zeros(nbSeg)
998 jj1 = np.zeros(nbSeg)
999 for i
in range(nbSeg):
1000 ip = (smallsize - i2t[i] + i1t[i] - 1) // 2
1001 jp = (smallsize - j2t[i] + j1t[i] - 1) // 2
1002 i1 = np.maximum(i1t[i] - ip, 0)
1003 j1 = np.maximum(j1t[i] - jp, 0)
1004 if (j1 + smallsize) > npt:
1005 j1 = npt - smallsize
1006 if (i1 + smallsize) > npt:
1007 i1 = npt - smallsize
1009 k = petalMap[i1 + smallsize // 2, j1 + smallsize // 2]
1010 petal = (petalMap == k)
1011 influ[:, :, k] = petal[i1:i1 + smallsize, j1:j1 + smallsize]
1015 return influ, ii1, jj1, int(smallsize), nbSeg
1020 Makes an image npt x npt of <nbSeg> regularly spaced angular segments
1021 centred on (i0, j0).
1022 Origin of angles is set by <pupAngleDegree>.
1024 The segments are oriented as defined in document "Standard Coordinates
1025 and Basic Conventions", ESO-193058.
1026 This document states that the X axis lies in the middle of a petal, i.e.
1027 that the axis Y is along the spider.
1028 The separation angle between segments are [-30, 30, 90, 150, -150, -90].
1029 For this reason, an <esoOffsetAngle> = -pi/6 is introduced in the code.
1032 pupAngleDegree = 5.0 # pupil rotation/clocking angle
1035 p = build_petals(nbSeg, pupAngleDegree, i0, j0, npt)
1038 rot = pupAngleDegree * np.pi / 180.0
1041 esoOffsetAngle = -np.pi / 6
1042 x = np.arange(npt) - i0
1043 y = np.arange(npt) - j0
1044 X, Y = np.meshgrid(x, y, indexing=
'ij')
1045 theta = (np.arctan2(Y, X) - rot + 2 * np.pi - esoOffsetAngle) % (2 * np.pi)
1048 angleStep = 2 * np.pi / nbSeg
1049 startAngle = np.arange(nbSeg) * angleStep
1050 endAngle = np.roll(startAngle, -1)
1051 endAngle[-1] = 2 * np.pi
1052 petalMap = np.zeros((npt, npt), dtype=int)
1053 for i
in range(nbSeg):
1054 nn = np.where(np.logical_and(theta >= startAngle[i], theta < endAngle[i]))
Parameter classes for COMPASS.
Numerical constants for shesha and config enumerations for safe-typing.
def make_pzt_dm(conf.Param_dm p_dm, conf.Param_geom p_geom, float cobs, float pupAngle)
Compute the actuators positions and the influence functions for a pzt DM.
def makePetalDm(p_dm, p_geom, pupAngleDegree)
def make_petal_dm_core(pupImage, pupAngleDegree)
<pupImage> : image of the pupil <pupAngleDegree> : rotation angle of the pupil in degrees
def comp_dmgeom(conf.Param_dm p_dm, conf.Param_geom p_geom)
Compute the geometry of a DM : positions of actuators and influence functions.
Dms dm_init(carmaWrap_context context, List[conf.Param_dm] p_dms, conf.Param_tel p_tel, conf.Param_geom p_geom, List[conf.Param_wfs] p_wfss=None)
Create and initialize a Dms object on the gpu.
def build_petals(nbSeg, pupAngleDegree, i0, j0, npt)
Makes an image npt x npt of <nbSeg> regularly spaced angular segments centred on (i0,...
None make_kl_dm(conf.Param_dm p_dm, int patchDiam, conf.Param_geom p_geom, float cobs)
Compute the influence function for a Karhunen-Loeve DM.
def init_custom_dm(conf.Param_dm p_dm, conf.Param_geom p_geom, float diam)
Read Fits for influence pzt fonction and form.
def make_tiptilt_dm(conf.Param_dm p_dm, int patchDiam, conf.Param_geom p_geom, float diam)
Compute the influence functions for a tip-tilt DM.
def dm_init_standalone(carmaWrap_context context, list p_dms, conf.Param_geom p_geom, diam=1., cobs=0., pupAngle=0., wfs_xpos=[0], wfs_ypos=[0])
Create and initialize a Dms object on the gpu.
def correct_dm(context, Dms dms, list p_dms, conf.Param_controller p_controller, conf.Param_geom p_geom, np.ndarray imat=None, dict dataBase={}, bool use_DB=False)
Correct the geometry of the DMs using the imat (filter unseen actuators)