COMPASS  5.0.0
End-to-end AO simulation tool using GPU acceleration
dm_init.py
1 
37 
38 import shesha.config as conf
39 import shesha.constants as scons
40 
41 from shesha.constants import CONST
42 
43 from shesha.util import dm_util, influ_util, kl_util
44 from shesha.util import hdf5_util as h5u
45 
46 import numpy as np
47 
48 import pandas as pd
49 from scipy import interpolate
50 from shesha.sutra_wrap import carmaWrap_context, Dms
51 
52 from typing import List
53 
54 from tqdm import tqdm
55 
56 
57 def dm_init(context: carmaWrap_context, p_dms: List[conf.Param_dm],
58  p_tel: conf.Param_tel, p_geom: conf.Param_geom,
59  p_wfss: List[conf.Param_wfs] = None, keepAllActu: bool = False) -> Dms:
60  """Create and initialize a Dms object on the gpu
61 
62  :parameters:
63  context: (carmaWrap_context): context
64  p_dms: (list of Param_dms) : dms settings
65  p_tel: (Param_tel) : telescope settings
66  p_geom: (Param_geom) : geom settings
67  p_wfss: (list of Param_wfs) : wfs settings
68  :return:
69  Dms: (Dms): Dms object
70  """
71  max_extent = 0
72  if (p_wfss is not None):
73  xpos_wfs = []
74  ypos_wfs = []
75  for i in range(len(p_wfss)):
76  xpos_wfs.append(p_wfss[i].xpos)
77  ypos_wfs.append(p_wfss[i].ypos)
78  else:
79  xpos_wfs = [0]
80  ypos_wfs = [0]
81  if (len(p_dms) != 0):
82  dms = Dms()
83  types_dm = [p_dm.type for p_dm in p_dms]
84  if scons.DmType.TT in types_dm:
85  first_TT = types_dm.index(scons.DmType.TT)
86  if np.any(np.array(types_dm[first_TT:]) != scons.DmType.TT):
87  raise RuntimeError("TT must be defined at the end of the dms parameters")
88 
89  for i in range(len(p_dms)):
90  max_extent = _dm_init(context, dms, p_dms[i], xpos_wfs, ypos_wfs, p_geom,
91  p_tel.diam, p_tel.cobs, p_tel.pupangle, max_extent,
92  keepAllActu=keepAllActu)
93 
94  return dms
95 
96 
97 def _dm_init(context: carmaWrap_context, dms: Dms, p_dm: conf.Param_dm, xpos_wfs: list,
98  ypos_wfs: list, p_geom: conf.Param_geom, diam: float, cobs: float,
99  pupAngle: float, max_extent: int, keepAllActu: bool = False):
100  """ inits a Dms object on the gpu
101 
102  :parameters:
103  context: (carmaWrap_context): context
104  dms: (Dms) : dm object
105 
106  p_dm: (Param_dms) : dm settings
107 
108  xpos_wfs: (list) : list of wfs xpos
109 
110  ypos_wfs: (list) : list of wfs ypos
111 
112  p_geom: (Param_geom) : geom settings
113 
114  diam: (float) : diameter of telescope
115 
116  cobs: (float) : cobs of telescope
117 
118  max_extent: (int) : maximum dimension of all dms
119 
120  :return:
121  max_extent: (int) : new maximum dimension of all dms
122 
123  """
124 
125  if (p_dm.pupoffset is not None):
126  p_dm._puppixoffset = p_dm.pupoffset / diam * p_geom.pupdiam
127  # For patchDiam
128  patchDiam = dm_util.dim_dm_patch(p_geom.pupdiam, diam, p_dm.type, p_dm.alt, xpos_wfs,
129  ypos_wfs)
130 
131  if (p_dm.type == scons.DmType.PZT):
132  if p_dm.file_influ_hdf5 == None:
133  p_dm._pitch = patchDiam / float(p_dm.nact - 1)
134  # + 2.5 pitch each side
135  extent = p_dm._pitch * (p_dm.nact + p_dm.pzt_extent)
136  p_dm._n1, p_dm._n2 = dm_util.dim_dm_support(p_geom.cent, extent,
137  p_geom.ssize)
138 
139  # calcul defaut influsize
140  make_pzt_dm(p_dm, p_geom, cobs, pupAngle, keepAllActu=keepAllActu)
141  else:
142  init_custom_dm(p_dm, p_geom, diam)
143 
144  # max_extent
145  max_extent = max(max_extent, p_dm._n2 - p_dm._n1 + 1)
146 
147  dim = max(p_dm._n2 - p_dm._n1 + 1, p_geom._mpupil.shape[0])
148  p_dm._dim_screen = dim
149  ninflupos = p_dm._influpos.size
150  n_npts = p_dm._ninflu.size #// 2
151  dms.add_dm(context, p_dm.type, p_dm.alt, dim, p_dm._ntotact, p_dm._influsize,
152  ninflupos, n_npts, p_dm.push4imat, 0, context.active_device)
153  #infludata = p_dm._influ.flatten()[p_dm._influpos]
154  dms.d_dms[-1].pzt_loadarrays(p_dm._influ, p_dm._influpos.astype(np.int32),
155  p_dm._ninflu, p_dm._influstart, p_dm._i1, p_dm._j1)
156 
157  elif (p_dm.type == scons.DmType.TT):
158 
159  if (p_dm.alt == 0) and (max_extent != 0):
160  extent = int(max_extent * 1.05)
161  if (extent % 2 != 0):
162  extent += 1
163  else:
164  extent = p_geom.pupdiam + 16
165  p_dm._n1, p_dm._n2 = dm_util.dim_dm_support(p_geom.cent, extent, p_geom.ssize)
166  # max_extent
167  max_extent = max(max_extent, p_dm._n2 - p_dm._n1 + 1)
168 
169  dim = p_dm._n2 - p_dm._n1 + 1
170  make_tiptilt_dm(p_dm, patchDiam, p_geom, diam)
171  p_dm._dim_screen = dim
172  dms.add_dm(context, p_dm.type, p_dm.alt, dim, 2, dim, 1, 1, p_dm.push4imat, 0,
173  context.active_device)
174  dms.d_dms[-1].tt_loadarrays(p_dm._influ)
175 
176  elif (p_dm.type == scons.DmType.KL):
177 
178  extent = p_geom.pupdiam + 16
179  p_dm._n1, p_dm._n2 = dm_util.dim_dm_support(p_geom.cent, extent, p_geom.ssize)
180  # max_extent
181  max_extent = max(max_extent, p_dm._n2 - p_dm._n1 + 1)
182 
183  dim = p_dm._n2 - p_dm._n1 + 1
184 
185  make_kl_dm(p_dm, patchDiam, p_geom, cobs)
186 
187  ninflu = p_dm.nkl
188  p_dm._dim_screen = dim
189 
190  dms.add_dm(context, p_dm.type, p_dm.alt, dim, p_dm.nkl, p_dm._ncp, p_dm._nr,
191  p_dm._npp, p_dm.push4imat, p_dm._ord.max(), context.active_device)
192 
193  dms.d_dms[-1].kl_loadarrays(p_dm._rabas, p_dm._azbas, p_dm._ord, p_dm._cr,
194  p_dm._cp)
195 
196  else:
197 
198  raise TypeError("This type of DM doesn't exist ")
199  # Verif
200  # res1 = pol2car(*y_dm(n)._klbas,gkl_sfi(*y_dm(n)._klbas, 1));
201  # res2 = yoga_getkl(g_dm,0.,1);
202 
203  return max_extent
204 
205 
206 def _dm_init_factorized(context: carmaWrap_context, dms: Dms, p_dm: conf.Param_dm,
207  xpos_wfs: list, ypos_wfs: list, p_geom: conf.Param_geom,
208  diam: float, cobs: float, pupAngle: float, max_extent: int,
209  keepAllActu: bool = False):
210  """ inits a Dms object on the gpu
211  NOTE: This is the
212 
213  :parameters:
214  context: (carmaWrap_context): context
215  dms: (Dms) : dm object
216 
217  p_dm: (Param_dms) : dm settings
218 
219  xpos_wfs: (list) : list of wfs xpos
220 
221  ypos_wfs: (list) : list of wfs ypos
222 
223  p_geom: (Param_geom) : geom settings
224 
225  diam: (float) : diameter of telescope
226 
227  cobs: (float) : cobs of telescope
228 
229  max_extent: (int) : maximum dimension of all dms
230 
231  :return:
232  max_extent: (int) : new maximum dimension of all dms
233 
234  """
235 
236  if (p_dm.pupoffset is not None):
237  p_dm._puppixoffset = p_dm.pupoffset / diam * p_geom.pupdiam
238  # For patchDiam
239  patchDiam = dm_util.dim_dm_patch(p_geom.pupdiam, diam, p_dm.type, p_dm.alt, xpos_wfs,
240  ypos_wfs)
241 
242  if (p_dm.type == scons.DmType.PZT) and p_dm.file_influ_hdf5 is not None:
243  init_custom_dm(p_dm, p_geom, diam)
244  else:
245  if (p_dm.type == scons.DmType.PZT):
246  p_dm._pitch = patchDiam / float(p_dm.nact - 1)
247  # + 2.5 pitch each side
248  extent = p_dm._pitch * (p_dm.nact + p_dm.pzt_extent)
249 
250  # calcul defaut influsize
251  make_pzt_dm(p_dm, p_geom, cobs, pupAngle, keepAllActu=keepAllActu)
252 
253  elif (p_dm.type == scons.DmType.TT):
254  if (p_dm.alt == 0) and (max_extent != 0):
255  extent = int(max_extent * 1.05)
256  if (extent % 2 != 0):
257  extent += 1
258  else:
259  extent = p_geom.pupdiam + 16
260 
261  elif (p_dm.type == scons.DmType.KL):
262  extent = p_geom.pupdiam + 16
263  else:
264  raise TypeError("This type of DM doesn't exist ")
265 
266  # Verif
267  # res1 = pol2car(*y_dm(n)._klbas,gkl_sfi(*y_dm(n)._klbas, 1));
268  # res2 = yoga_getkl(g_dm,0.,1);
269 
270  p_dm._n1, p_dm._n2 = dm_util.dim_dm_support(p_geom.cent, extent, p_geom.ssize)
271 
272  # max_extent
273  max_extent = max(max_extent, p_dm._n2 - p_dm._n1 + 1)
274 
275  dim = max(p_dm._n2 - p_dm._n1 + 1, p_geom._mpupil.shape[0])
276 
277  if (p_dm.type == scons.DmType.PZT):
278  ninflupos = p_dm._influpos.size
279  n_npts = p_dm._ninflu.size #// 2
280  dms.add_dm(context, p_dm.type, p_dm.alt, dim, p_dm._ntotact, p_dm._influsize,
281  ninflupos, n_npts, p_dm.push4imat, 0, context.active_device)
282  #infludata = p_dm._influ.flatten()[p_dm._influpos]
283  dms.d_dms[-1].pzt_loadarrays(p_dm._influ, p_dm._influpos.astype(np.int32),
284  p_dm._ninflu, p_dm._influstart, p_dm._i1, p_dm._j1)
285  elif (p_dm.type == scons.DmType.TT):
286  make_tiptilt_dm(p_dm, patchDiam, p_geom, diam)
287  dms.add_dm(context, p_dm.type, p_dm.alt, dim, 2, dim, 1, 1, p_dm.push4imat, 0,
288  context.active_device)
289  dms.d_dms[-1].tt_loadarrays(p_dm._influ)
290  elif (p_dm.type == scons.DmType.KL):
291  make_kl_dm(p_dm, patchDiam, p_geom, cobs)
292  ninflu = p_dm.nkl
293 
294  dms.add_dm(context, p_dm.type, p_dm.alt, dim, p_dm.nkl, p_dm._ncp, p_dm._nr,
295  p_dm._npp, p_dm.push4imat, p_dm._ord.max(), context.active_device)
296  dms.d_dms[-1].kl_loadarrays(p_dm._rabas, p_dm._azbas, p_dm._ord, p_dm._cr,
297  p_dm._cp)
298 
299  return max_extent
300 
301 
302 def dm_init_standalone(context: carmaWrap_context, p_dms: list, p_geom: conf.Param_geom,
303  diam=1., cobs=0., pupAngle=0., wfs_xpos=[0], wfs_ypos=[0]):
304  """Create and initialize a Dms object on the gpu
305 
306  :parameters:
307  p_dms: (list of Param_dms) : dms settings
308 
309  p_geom: (Param_geom) : geom settings
310 
311  diam: (float) : diameter of telescope (default 1.)
312 
313  cobs: (float) : cobs of telescope (default 0.)
314 
315  pupAngle: (float) : pupil rotation angle (degrees, default 0.)
316 
317  wfs_xpos: (array) : guide star x position on sky (in arcsec).
318 
319  wfs_ypos: (array) : guide star y position on sky (in arcsec).
320 
321  """
322  max_extent = [0]
323  if (len(p_dms) != 0):
324  dms = Dms()
325  for i in range(len(p_dms)):
326  _dm_init(context, dms, p_dms[i], wfs_xpos, wfs_ypos, p_geom, diam, cobs,
327  pupAngle, max_extent)
328  return dms
329 
330 
331 def make_pzt_dm(p_dm: conf.Param_dm, p_geom: conf.Param_geom, cobs: float,
332  pupAngle: float, keepAllActu: bool = False):
333  """Compute the actuators positions and the influence functions for a pzt DM.
334  NOTE: if the DM is in altitude, central obstruction is forced to 0
335 
336  :parameters:
337  p_dm: (Param_dm) : dm parameters
338 
339  p_geom: (Param_geom) : geometry parameters
340 
341  cobs: (float) : telescope central obstruction
342 
343  :return:
344  influ: (np.ndarray(dims=3, dtype=np.float64)) : cube of the IF for each actuator
345 
346  """
347  # best parameters, as determined by a multi-dimensional fit
348  #(see coupling3.i)
349  coupling = p_dm.coupling
350 
351  # prepare to compute IF on partial (local) support of size <smallsize>
352  pitch = p_dm._pitch
353  smallsize = 0
354 
355  # Petal DM (segmentation of M4)
356  if (p_dm.influ_type == scons.InfluType.PETAL):
357  makePetalDm(p_dm, p_geom, pupAngle)
358  return
359 
360  if (p_dm.influ_type == scons.InfluType.RADIALSCHWARTZ):
361  smallsize = influ_util.makeRadialSchwartz(pitch, coupling)
362  elif (p_dm.influ_type == scons.InfluType.SQUARESCHWARTZ):
363  smallsize = influ_util.makeSquareSchwartz(pitch, coupling)
364  elif (p_dm.influ_type == scons.InfluType.BLACKNUTT):
365  smallsize = influ_util.makeBlacknutt(pitch, coupling)
366  elif (p_dm.influ_type == scons.InfluType.GAUSSIAN):
367  smallsize = influ_util.makeGaussian(pitch, coupling)
368  elif (p_dm.influ_type == scons.InfluType.BESSEL):
369  smallsize = influ_util.makeBessel(pitch, coupling, p_dm.type_pattern)
370  elif (p_dm.influ_type == scons.InfluType.DEFAULT):
371  smallsize = influ_util.makeRigaut(pitch, coupling)
372  else:
373  print("ERROR influtype not recognized ")
374  p_dm._influsize = smallsize
375 
376  # compute location (x,y and i,j) of each actuator:
377  nxact = p_dm.nact
378 
379  if p_dm.type_pattern is None:
380  p_dm.type_pattern = scons.PatternType.SQUARE
381 
382  if p_dm.type_pattern == scons.PatternType.HEXA:
383  print("Pattern type : hexa")
384  cub = dm_util.createHexaPattern(pitch, p_geom.pupdiam * 1.1)
385  keepAllActu = True
386  elif p_dm.type_pattern == scons.PatternType.HEXAM4:
387  print("Pattern type : hexaM4")
388  keepAllActu = True
389  cub = dm_util.createDoubleHexaPattern(pitch, p_geom.pupdiam * 1.1, pupAngle)
390  if p_dm.margin_out is not None:
391  print(f'p_dm.margin_out={p_dm.margin_out} is being '
392  'used for pupil-based actuator filtering')
393  pup_side = p_geom._ipupil.shape[0]
394  cub_off = dm_util.filterActuWithPupil(cub + pup_side // 2 - 0.5,
395  p_geom._ipupil,
396  p_dm.margin_out * p_dm.get_pitch())
397  cub = cub_off - pup_side // 2 + 0.5
398  p_dm.set_ntotact(cub.shape[1])
399  elif p_dm.type_pattern == scons.PatternType.SQUARE:
400  print("Pattern type : square")
401  cub = dm_util.createSquarePattern(pitch, nxact + 4)
402  else:
403  raise ValueError("This pattern does not exist for pzt dm")
404 
405  if keepAllActu:
406  inbigcirc = np.arange(cub.shape[1])
407  else:
408  if (p_dm.alt > 0):
409  cobs = 0
410  inbigcirc = dm_util.select_actuators(cub[0, :], cub[1, :], p_dm.nact,
411  p_dm._pitch, cobs, p_dm.margin_in,
412  p_dm.margin_out, p_dm._ntotact)
413  p_dm._ntotact = inbigcirc.size
414 
415  # print(('inbigcirc',inbigcirc.shape))
416 
417  # converting to array coordinates:
418  cub += p_geom.cent
419 
420  # filtering actuators outside of a disk radius = rad (see above)
421  cubval = cub[:, inbigcirc]
422  ntotact = cubval.shape[1]
423  #pfits.writeto("cubeval.fits", cubval)
424  xpos = cubval[0, :]
425  ypos = cubval[1, :]
426  i1t = (cubval[0, :] - smallsize / 2 - 0.5 - p_dm._n1).astype(np.int32)
427  j1t = (cubval[1, :] - smallsize / 2 - 0.5 - p_dm._n1).astype(np.int32)
428 
429  p_dm._xpos = xpos
430  p_dm._ypos = ypos
431  p_dm._i1 = i1t
432  p_dm._j1 = j1t
433 
434  # Allocate array of influence functions
435 
436  influ = np.zeros((smallsize, smallsize, ntotact), dtype=np.float32)
437  # Computation of influence function for each actuator
438 
439  print("Computing Influence Function type : ", p_dm.influ_type)
440 
441  for i in tqdm(range(ntotact)):
442 
443  i1 = i1t[i]
444  x = np.tile(np.arange(i1, i1 + smallsize, dtype=np.float32),
445  (smallsize, 1)).T # pixel coords in ref frame "dm support"
446  # pixel coords in ref frame "pupil support"
447  x += p_dm._n1
448  # pixel coords in local ref frame
449  x -= xpos[i]
450 
451  j1 = j1t[i]
452  y = np.tile(np.arange(j1, j1 + smallsize, dtype=np.float32),
453  (smallsize, 1)) # idem as X, in Y
454  y += p_dm._n1
455  y -= ypos[i]
456  # print("Computing Influence Function #%d/%d \r" % (i, ntotact), end=' ')
457 
458  if (p_dm.influ_type == scons.InfluType.RADIALSCHWARTZ):
459  influ[:, :, i] = influ_util.makeRadialSchwartz(pitch, coupling, x=x, y=y)
460  elif (p_dm.influ_type == scons.InfluType.SQUARESCHWARTZ):
461  influ[:, :, i] = influ_util.makeSquareSchwartz(pitch, coupling, x=x, y=y)
462  elif (p_dm.influ_type == scons.InfluType.BLACKNUTT):
463  influ[:, :, i] = influ_util.makeBlacknutt(pitch, coupling, x=x, y=y)
464  elif (p_dm.influ_type == scons.InfluType.GAUSSIAN):
465  influ[:, :, i] = influ_util.makeGaussian(pitch, coupling, x=x, y=y)
466  elif (p_dm.influ_type == scons.InfluType.BESSEL):
467  influ[:, :, i] = influ_util.makeBessel(pitch, coupling, x=x, y=y,
468  patternType=p_dm.type_pattern)
469  elif (p_dm.influ_type == scons.InfluType.DEFAULT):
470  influ[:, :, i] = influ_util.makeRigaut(pitch, coupling, x=x, y=y)
471  else:
472  print("ERROR influtype not recognized (defaut or gaussian or bessel)")
473 
474  if (p_dm._puppixoffset is not None):
475  xpos += p_dm._puppixoffset[0]
476  ypos += p_dm._puppixoffset[1]
477 
478  if p_dm.segmented_mirror:
479  # mpupil to ipupil shift
480  # Good centering assumptions
481  s = (p_geom._ipupil.shape[0] - p_geom._mpupil.shape[0]) // 2
482 
483  from skimage.morphology import label
484  k = 0
485  for i in tqdm(range(ntotact)):
486  # Pupil area corresponding to influ data
487  i1, j1 = i1t[i] + s - smallsize // 2, j1t[i] + s - smallsize // 2
488  pupilSnapshot = p_geom._ipupil[i1:i1 + smallsize, j1:j1 + smallsize]
489  if np.all(pupilSnapshot): # We have at least one non-pupil pixel
490  continue
491  labels, num = label(pupilSnapshot, background=0, return_num=True)
492  if num <= 1:
493  continue
494  k += 1
495  maxPerArea = np.array([
496  (influ[:, :, i] * (labels == k).astype(np.float32)).max()
497  for k in range(1, num + 1)
498  ])
499  influ[:, :, i] *= (labels == np.argmax(maxPerArea) + 1).astype(np.float32)
500  print(f'{k} cross-spider influence functions trimmed.')
501 
502  influ = influ * float(p_dm.unitpervolt / np.max(influ))
503 
504  p_dm._influ = influ
505 
506  comp_dmgeom(p_dm, p_geom)
507 
508  dim = max(p_geom._mpupil.shape[0], p_dm._n2 - p_dm._n1 + 1)
509  off = (dim - p_dm._influsize) // 2
510 
511 
512 def init_custom_dm(p_dm: conf.Param_dm, p_geom: conf.Param_geom, diam: float):
513  """Read Fits for influence pzt fonction and form
514 
515  :parameters:
516  p_dm: (Param_dm) : dm settings
517 
518  p_geom: (Param_geom) : geom settings
519 
520  diam: (float) : tel diameter
521 
522  Conversion. There are sereval coordinate systems.
523  Some are coming from the input fits file, others from compass.
524  Those systems differ by scales and offsets.
525  Coord systems from the input fits file:
526  - they all have the same scale: the coordinates are expressed in
527  pixels
528  - one system is the single common frame where fits data are described [i]
529  - one is local to the minimap [l]
530  Coord systems from compass:
531  - they all have the same scale (different from fits one) expressed in
532  pixels
533  - one system is attached to ipupil (i=image, largest support) [i]
534  - one system is attached to mpupil (m=medium, medium support) [m]
535  - one system is local to minimap [l)]
536 
537  Variables will be named using
538  f, c: define either fits or compass
539 
540  """
541  from astropy.io import fits as pfits
542 
543  # read fits file
544  hdul = pfits.open(p_dm.file_influ_hdf5)
545  print("Read influence function from fits file : ", p_dm.file_influ_hdf5)
546 
547  xC = hdul[0].header['XCENTER']
548  yC = hdul[0].header['YCENTER']
549  pitchPix = hdul[0].header['PITCHPX']
550  pitchMeters = hdul[0].header['PITCHM']
551  hi_i1, hi_j1 = hdul[1].data
552  influ = hdul[2].data
553  xpos, ypos = hdul[3].data
554 
555  # facteur de conversion des coordonnees du fichier Fits vers les pixels
556  # de compass
557  # pitchPixCompass = p_dm._pitch # valeur du pitch entre actus en pixels de compass
558  pitchPixCompass = pitchMeters / p_geom._pixsize
559  scaleToCompass = pitchPixCompass / pitchPix
560 
561 
565  offsetXToCompass = p_geom.cent - xC * scaleToCompass
566  offsetYToCompass = p_geom.cent - yC * scaleToCompass
567 
568 
569  iSize, jSize, ntotact = influ.shape
570  if iSize != jSize:
571  raise ('Error')
572 
573  ess1 = np.ceil((hi_i1+iSize) * scaleToCompass + offsetXToCompass) \
574  - np.floor(hi_i1 * scaleToCompass + offsetXToCompass)
575  ess1 = np.max(ess1)
576 
577  ess2 = np.ceil((hi_j1+jSize) * scaleToCompass + offsetYToCompass) \
578  - np.floor(hi_j1 * scaleToCompass + offsetYToCompass)
579  ess2 = np.max(ess2)
580  smallsize = np.maximum(ess1, ess2).astype(int)
581 
582  # Allocate influence function maps and other arrays
583  p_dm._ntotact = ntotact
584  p_dm._influsize = np.int(smallsize)
585  p_dm._i1 = np.zeros(ntotact, dtype=np.int32)
586  p_dm._j1 = np.zeros(ntotact, dtype=np.int32)
587  p_dm._xpos = np.zeros(ntotact, dtype=np.float32)
588  p_dm._ypos = np.zeros(ntotact, dtype=np.float32)
589  p_dm._influ = np.zeros((smallsize, smallsize, ntotact), dtype=np.float32)
590 
591  # loop on actuators
592  for i in range(ntotact):
593  # coordonnes en pixels de la minicarte dans le systeme de depart
594  hi_x = np.arange(iSize) + hi_i1[i]
595  hi_y = np.arange(jSize) + hi_j1[i]
596 
597  # transfert des coord d'origine vers systeme compass
598  ci_x = hi_x * scaleToCompass + offsetXToCompass
599  ci_y = hi_y * scaleToCompass + offsetYToCompass
600 
601  # Creation des coord de destination dans systeme compass
602  ci_i1 = hi_i1[
603  i] * scaleToCompass + offsetXToCompass # itruc = truc dans repere ipup
604  ci_j1 = hi_j1[i] * scaleToCompass + offsetYToCompass
605  ci_i1 = np.floor(ci_i1).astype(np.int32)
606  ci_j1 = np.floor(ci_j1).astype(np.int32)
607  ci_xpix = ci_i1 + np.arange(smallsize)
608  ci_ypix = ci_j1 + np.arange(smallsize)
609  # WARNING: les xpos et ypos sont approximatifs !! Bon pour debug only ...
610  # p_dm._xpos[i] = ci_i1 + smallsize/2.0
611  # p_dm._ypos[i] = ci_j1 + smallsize/2.0
612 
613  f = interpolate.interp2d(ci_y, ci_x, influ[:, :, i], kind='cubic')
614  temp = f(ci_ypix, ci_xpix) * p_dm.unitpervolt
615  # temp[np.where(temp<1e-6)] = 0.
616  p_dm._influ[:, :, i] = temp
617 
618  # ....
619  p_dm._i1[i] = ci_i1
620  p_dm._j1[i] = ci_j1
621 
622  tmp = p_geom.ssize - (np.max(p_dm._i1 + smallsize))
623  margin_i = np.minimum(tmp, np.min(p_dm._i1))
624  tmp = p_geom.ssize - (np.max(p_dm._j1 + smallsize))
625  margin_j = np.minimum(tmp, np.min(p_dm._j1))
626 
627  p_dm._xpos = xpos * scaleToCompass + offsetXToCompass
628  p_dm._ypos = ypos * scaleToCompass + offsetYToCompass
629 
630  p_dm._n1 = int(np.minimum(margin_i, margin_j))
631  p_dm._n2 = p_geom.ssize - 1 - p_dm._n1
632 
633  p_dm._i1 -= p_dm._n1
634  p_dm._j1 -= p_dm._n1
635 
636  comp_dmgeom(p_dm, p_geom)
637 
638 
639 def make_tiptilt_dm(p_dm: conf.Param_dm, patchDiam: int, p_geom: conf.Param_geom,
640  diam: float):
641  """Compute the influence functions for a tip-tilt DM
642 
643  :parameters:
644  p_dm: (Param_dm) : dm settings
645 
646  patchDiam: (int) : patchDiam for dm size
647 
648  p_geom: (Param_geom) : geom settings
649 
650  diam: (float) : telescope diameter
651  :return:
652  influ: (np.ndarray(dims=3,dtype=np.float64)) : cube of the IF
653 
654  """
655  dim = max(p_dm._n2 - p_dm._n1 + 1, p_geom._mpupil.shape[0])
656  #norms = [np.linalg.norm([w.xpos, w.ypos]) for w in p_wfs]
657 
658  nzer = 2
659  influ = dm_util.make_zernike(nzer + 1, dim, patchDiam, p_geom.cent - p_dm._n1 + 1,
660  p_geom.cent - p_dm._n1 + 1, 1)[:, :, 1:]
661 
662  # normalization factor: one unit of tilt gives 1 arcsec:
663  current = influ[dim // 2 - 1, dim // 2 - 1, 0] - \
664  influ[dim // 2 - 2, dim // 2 - 2, 0]
665  fact = p_dm.unitpervolt * diam / p_geom.pupdiam * 4.848 / current
666 
667  influ = influ * fact
668  p_dm._ntotact = influ.shape[2]
669  p_dm._influsize = influ.shape[0]
670  p_dm._influ = influ
671 
672  return influ
673 
674 
675 def make_kl_dm(p_dm: conf.Param_dm, patchDiam: int, p_geom: conf.Param_geom,
676  cobs: float) -> None:
677  """Compute the influence function for a Karhunen-Loeve DM
678 
679  :parameters:
680  p_dm: (Param_dm) : dm settings
681 
682  patchDiam: (int) : patchDiam for dm size
683 
684  p_geom: (Param_geom) : geom settings
685 
686  cobs: (float) : telescope cobs
687 
688  """
689  dim = p_geom._mpupil.shape[0]
690 
691  print("KL type: ", p_dm.type_kl)
692 
693  if (p_dm.nkl < 13):
694  nr = np.long(5.0 * np.sqrt(52)) # one point per degree
695  npp = np.long(10.0 * nr)
696  else:
697  nr = np.long(5.0 * np.sqrt(p_dm.nkl))
698  npp = np.long(10.0 * nr)
699 
700  radp = kl_util.make_radii(cobs, nr)
701 
702  kers = kl_util.make_kernels(cobs, nr, radp, p_dm.type_kl, p_dm.outscl)
703 
704  evals, nord, npo, ordd, rabas = kl_util.gkl_fcom(kers, cobs, p_dm.nkl)
705 
706  azbas = kl_util.make_azimuth(nord, npp)
707 
708  ncp, ncmar, px, py, cr, cp, pincx, pincy, pincw, ap = kl_util.set_pctr(
709  patchDiam, nr, npp, p_dm.nkl, cobs, nord)
710 
711  p_dm._ntotact = p_dm.nkl
712  p_dm._nr = nr # number of radial points
713  p_dm._npp = npp # number of elements
714  p_dm._ord = ordd # the radial orders of the basis
715  p_dm._rabas = rabas # the radial array of the basis
716  p_dm._azbas = azbas # the azimuthal array of the basis
717  p_dm._ncp = ncp # dim of grid
718  p_dm._cr = cr # radial coord in cartesien grid
719  p_dm._cp = cp # phi coord in cartesien grid
720  p_dm._i1 = np.zeros((p_dm.nkl), dtype=np.int32) + \
721  (dim - patchDiam) // 2
722  p_dm._j1 = np.zeros((p_dm.nkl), dtype=np.int32) + \
723  (dim - patchDiam) // 2
724  p_dm._ntotact = p_dm.nkl
725  p_dm.ap = ap
726 
727 
728 def comp_dmgeom(p_dm: conf.Param_dm, p_geom: conf.Param_geom):
729  """Compute the geometry of a DM : positions of actuators and influence functions
730 
731  :parameters:
732  dm: (Param_dm) : dm settings
733 
734  geom: (Param_geom) : geom settings
735  """
736  smallsize = p_dm._influsize
737  nact = p_dm._ntotact
738  dm_dim = int(p_dm._n2 - p_dm._n1 + 1)
739  mpup_dim = p_geom._mpupil.shape[0]
740 
741  if (dm_dim < mpup_dim):
742  offs = (mpup_dim - dm_dim) // 2
743  else:
744  offs = 0
745  mpup_dim = dm_dim
746 
747  indgen = np.tile(np.arange(smallsize, dtype=np.int32), (smallsize, 1))
748 
749  tmpx = np.tile(indgen, (nact, 1, 1))
750  tmpy = np.tile(indgen.T, (nact, 1, 1))
751 
752  tmpx += offs + p_dm._i1[:, None, None]
753  tmpy += offs + p_dm._j1[:, None, None]
754 
755  tmp = tmpx + mpup_dim * tmpy
756 
757  # bug in limit of def zone -10 destoe influpos for all actuator
758  tmp[tmpx < 0] = mpup_dim * mpup_dim + 10 # -10
759  tmp[tmpy < 0] = mpup_dim * mpup_dim + 10 # -10
760  tmp[tmpx > dm_dim - 1] = mpup_dim * mpup_dim + 10 # -10
761  tmp[tmpy > dm_dim - 1] = mpup_dim * mpup_dim + 10 # -10
762  itmps = np.argsort(tmp.flatten()).astype(np.int32)
763  tmps = tmp.flatten()[itmps].astype(np.int32)
764  itmps = itmps[np.where(itmps > -1)]
765 
766  istart = np.zeros((mpup_dim * mpup_dim), dtype=np.int32)
767  npts = np.zeros((mpup_dim * mpup_dim), dtype=np.int32)
768 
769  tmps_unique, cpt = np.unique(tmps, return_counts=True)
770  if (tmps_unique > npts.size - 1).any():
771  tmps_unique = tmps_unique[:-1]
772  cpt = cpt[:-1]
773 
774  for i in range(tmps_unique.size):
775  npts[tmps_unique[i]] = cpt[i]
776  istart[1:] = np.cumsum(npts[:-1])
777 
778  p_dm._influpos = itmps[:np.sum(npts)].astype(np.int32)
779  # infludata = p_dm._influ.flatten()[p_dm._influpos]
780  # p_dm._influ = infludata[:,None,None]
781  # p_dm._influpos = p_dm._influpos / (smallsize * smallsize)
782  p_dm._ninflu = npts.astype(np.int32)
783  p_dm._influstart = istart.astype(np.int32)
784 
785  p_dm._i1 += offs
786  p_dm._j1 += offs
787 
788  # ninflu = np.zeros((istart.size * 2))
789  # ninflu[::2] = istart.astype(np.int32)
790  # ninflu[1::2] = npts.astype(np.int32)
791 
792  # p_dm._ninflu = ninflu
793 
794 
795 def correct_dm(context, dms: Dms, p_dms: list, p_controller: conf.Param_controller,
796  p_geom: conf.Param_geom, imat: np.ndarray = None, dataBase: dict = {},
797  use_DB: bool = False):
798  """Correct the geometry of the DMs using the imat (filter unseen actuators)
799 
800  :parameters:
801  context: (carmaWrap_context): context
802  dms: (Dms) : Dms object
803  p_dms: (list of Param_dm) : dms settings
804  p_controller: (Param_controller) : controller settings
805  p_geom: (Param_geom) : geom settings
806  imat: (np.ndarray) : interaction matrix
807  dataBase: (dict): dictionary containing paths to files to load
808  use_DB: (bool): dataBase use flag
809  """
810  print("Filtering unseen actuators... ")
811  if imat is not None:
812  resp = np.sqrt(np.sum(imat**2, axis=0))
813 
814  ndm = p_controller.ndm.size
815  inds = 0
816 
817  for nmc in range(ndm):
818  nm = p_controller.ndm[nmc]
819  nactu_nm = p_dms[nm]._ntotact
820  # filter actuators only in stackarray mirrors:
821  if (p_dms[nm].type == scons.DmType.PZT):
822  if "dm" in dataBase:
823  influpos, ninflu, influstart, i1, j1, ok = h5u.load_dm_geom_from_dataBase(
824  dataBase, nmc)
825  p_dms[nm].set_ntotact(ok.shape[0])
826  p_dms[nm].set_influ(p_dms[nm]._influ[:, :, ok.tolist()])
827  p_dms[nm].set_xpos(p_dms[nm]._xpos[ok])
828  p_dms[nm].set_ypos(p_dms[nm]._ypos[ok])
829  p_dms[nm]._influpos = influpos
830  p_dms[nm]._ninflu = ninflu
831  p_dms[nm]._influstart = influstart
832  p_dms[nm]._i1 = i1
833  p_dms[nm]._j1 = j1
834  else:
835  tmp = resp[inds:inds + p_dms[nm]._ntotact]
836  ok = np.where(tmp > p_dms[nm].thresh * np.max(tmp))[0]
837  nok = np.where(tmp <= p_dms[nm].thresh * np.max(tmp))[0]
838 
839  p_dms[nm].set_ntotact(ok.shape[0])
840  p_dms[nm].set_influ(p_dms[nm]._influ[:, :, ok.tolist()])
841  p_dms[nm].set_xpos(p_dms[nm]._xpos[ok])
842  p_dms[nm].set_ypos(p_dms[nm]._ypos[ok])
843  p_dms[nm].set_i1(p_dms[nm]._i1[ok])
844  p_dms[nm].set_j1(p_dms[nm]._j1[ok])
845 
846  comp_dmgeom(p_dms[nm], p_geom)
847  if use_DB:
848  h5u.save_dm_geom_in_dataBase(nmc, p_dms[nm]._influpos,
849  p_dms[nm]._ninflu,
850  p_dms[nm]._influstart, p_dms[nm]._i1,
851  p_dms[nm]._j1, ok)
852 
853  dim = max(p_dms[nm]._n2 - p_dms[nm]._n1 + 1, p_geom._mpupil.shape[0])
854  ninflupos = p_dms[nm]._influpos.size
855  n_npts = p_dms[nm]._ninflu.size
856  dms.remove_dm(nm)
857  dms.insert_dm(context, p_dms[nm].type, p_dms[nm].alt, dim,
858  p_dms[nm]._ntotact, p_dms[nm]._influsize, ninflupos, n_npts,
859  p_dms[nm].push4imat, 0, p_dms[nm].dx / p_geom._pixsize,
860  p_dms[nm].dy / p_geom._pixsize, p_dms[nm].theta, p_dms[nm].G,
861  context.active_device, nm)
862  dms.d_dms[nm].pzt_loadarrays(p_dms[nm]._influ, p_dms[nm]._influpos.astype(
863  np.int32), p_dms[nm]._ninflu, p_dms[nm]._influstart, p_dms[nm]._i1,
864  p_dms[nm]._j1)
865 
866  inds += nactu_nm
867  print("Done")
868 
869 
870 def makePetalDm(p_dm, p_geom, pupAngleDegree):
871  '''
872  makePetalDm(p_dm, p_geom, pupAngleDegree)
873 
874  The function builds a DM, segmented in petals according to the pupil
875  shape. The petals will be adapted to the EELT case only.
876 
877  <p_geom> : compass object p_geom. The function requires the object p_geom
878  in order to know what is the pupil mask, and what is the mpupil.
879  <p_dm> : compass petal dm object p_dm to be created. The function will
880  transform/modify in place the attributes of the object p_dm.
881 
882 
883  '''
884  p_dm._n1 = p_geom._n1
885  p_dm._n2 = p_geom._n2
886  influ, i1, j1, smallsize, nbSeg = make_petal_dm_core(p_geom._mpupil, pupAngleDegree)
887  p_dm._influsize = smallsize
888  p_dm.set_ntotact(nbSeg)
889  p_dm._i1 = i1
890  p_dm._j1 = j1
891  p_dm._xpos = i1 + smallsize / 2 + p_dm._n1
892  p_dm._ypos = j1 + smallsize / 2 + p_dm._n1
893  p_dm._influ = influ
894 
895  # generates the arrays of indexes for the GPUs
896  comp_dmgeom(p_dm, p_geom)
897 
898 
899 def make_petal_dm_core(pupImage, pupAngleDegree):
900  """
901  <pupImage> : image of the pupil
902 
903  La fonction renvoie des fn d'influence en forme de petale d'apres
904  une image de la pupille, qui est supposee etre segmentee.
905 
906 
907  influ, i1, j1, smallsize, nbSeg = make_petal_dm_core(pupImage, 0.0)
908  """
909  # Splits the pupil into connex areas.
910  # <segments> is the map of the segments, <nbSeg> in their number.
911  # binary_opening() allows us to suppress individual pixels that could
912  # be identified as relevant connex areas
913  from scipy.ndimage.measurements import label
914  from scipy.ndimage.morphology import binary_opening
915  s = np.ones((2, 2), dtype=np.bool)
916  segments, nbSeg = label(binary_opening(pupImage, s))
917 
918  # Faut trouver le plus petit support commun a tous les
919  # petales : on determine <smallsize>
920  smallsize = 0
921  i1t = [] # list of starting indexes of influ functions
922  j1t = []
923  i2t = [] # list of ending indexes of influ functions
924  j2t = []
925  for i in range(nbSeg):
926  petal = segments == (i + 1) # identification (boolean) of a given segment
927  profil = np.sum(petal, axis=1) != 0
928  extent = np.sum(profil).astype(np.int32)
929  i1t.append(np.min(np.where(profil)[0]))
930  i2t.append(np.max(np.where(profil)[0]))
931  if extent > smallsize:
932  smallsize = extent
933 
934  profil = np.sum(petal, axis=0) != 0
935  extent = np.sum(profil).astype(np.int32)
936  j1t.append(np.min(np.where(profil)[0]))
937  j2t.append(np.max(np.where(profil)[0]))
938  if extent > smallsize:
939  smallsize = extent
940 
941  # extension de la zone minimale pour avoir un peu de marge
942  smallsize += 2
943 
944  # Allocate array of influence functions
945  influ = np.zeros((smallsize, smallsize, nbSeg), dtype=np.float32)
946 
947  npt = pupImage.shape[0]
948  i0 = j0 = npt / 2 - 0.5
949  petalMap = build_petals(nbSeg, pupAngleDegree, i0, j0, npt)
950  ii1 = np.zeros(nbSeg)
951  jj1 = np.zeros(nbSeg)
952  for i in range(nbSeg):
953  ip = (smallsize - i2t[i] + i1t[i] - 1) // 2
954  jp = (smallsize - j2t[i] + j1t[i] - 1) // 2
955  i1 = np.maximum(i1t[i] - ip, 0)
956  j1 = np.maximum(j1t[i] - jp, 0)
957  if (j1 + smallsize) > npt:
958  j1 = npt - smallsize
959  if (i1 + smallsize) > npt:
960  i1 = npt - smallsize
961  #petal = segments==(i+1) # determine le segment pupille veritable
962  k = petalMap[i1 + smallsize // 2, j1 + smallsize // 2]
963  petal = (petalMap == k)
964  influ[:, :, k] = petal[i1:i1 + smallsize, j1:j1 + smallsize]
965  ii1[k] = i1
966  jj1[k] = j1
967 
968  return influ, ii1, jj1, int(smallsize), nbSeg
969 
970 
971 def build_petals(nbSeg, pupAngleDegree, i0, j0, npt):
972  """
973  Makes an image npt x npt of <nbSeg> regularly spaced angular segments
974  centred on (i0, j0).
975  Origin of angles is set by <pupAngleDegree>.
976 
977  The segments are oriented as defined in document "Standard Coordinates
978  and Basic Conventions", ESO-193058.
979  This document states that the X axis lies in the middle of a petal, i.e.
980  that the axis Y is along the spider.
981  The separation angle between segments are [-30, 30, 90, 150, -150, -90].
982  For this reason, an <esoOffsetAngle> = -pi/6 is introduced in the code.
983 
984  nbSeg = 6
985  pupAngleDegree = 5.0
986  i0 = j0 = 112.3
987  npt = 222
988  p = build_petals(nbSeg, pupAngleDegree, i0, j0, npt)
989  """
990  # conversion to radians
991  rot = pupAngleDegree * np.pi / 180.0
992 
993  # building coordinate maps
994  esoOffsetAngle = -np.pi / 6 # -30°, ESO definition.
995  x = np.arange(npt) - i0
996  y = np.arange(npt) - j0
997  X, Y = np.meshgrid(x, y, indexing='ij')
998  theta = (np.arctan2(Y, X) - rot + 2 * np.pi - esoOffsetAngle) % (2 * np.pi)
999 
1000  # Compute separation angle between segments: start and end.
1001  angleStep = 2 * np.pi / nbSeg
1002  startAngle = np.arange(nbSeg) * angleStep
1003  endAngle = np.roll(startAngle, -1)
1004  endAngle[-1] = 2 * np.pi # last angle is 0.00 and must be replaced by 2.pi
1005  petalMap = np.zeros((npt, npt), dtype=int)
1006  for i in range(nbSeg):
1007  nn = np.where(np.logical_and(theta >= startAngle[i], theta < endAngle[i]))
1008  petalMap[nn] = i
1009  return petalMap
shesha.init.dm_init.make_petal_dm_core
def make_petal_dm_core(pupImage, pupAngleDegree)
<pupImage> : image of the pupil
Definition: dm_init.py:912
shesha.util
Utilities functions.
Definition: shesha/shesha/util/__init__.py:1
shesha.sutra_wrap
Definition: sutra_wrap.py:1
shesha.init.dm_init.dm_init_standalone
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.
Definition: dm_init.py:324
shesha.constants
Numerical constants for shesha and config enumerations for safe-typing.
Definition: constants.py:1
shesha.init.dm_init.comp_dmgeom
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.
Definition: dm_init.py:739
shesha.config
Parameter classes for COMPASS.
Definition: shesha/shesha/config/__init__.py:1
shesha.init.dm_init.make_pzt_dm
def make_pzt_dm(conf.Param_dm p_dm, conf.Param_geom p_geom, float cobs, float pupAngle, bool keepAllActu=False)
Compute the actuators positions and the influence functions for a pzt DM.
Definition: dm_init.py:349
shesha.init.dm_init.dm_init
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, bool keepAllActu=False)
Create and initialize a Dms object on the gpu.
Definition: dm_init.py:68
shesha.init.dm_init.make_tiptilt_dm
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.
Definition: dm_init.py:657
shesha.init.dm_init.build_petals
def build_petals(nbSeg, pupAngleDegree, i0, j0, npt)
Definition: dm_init.py:993
shesha.init.dm_init.correct_dm
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)
Definition: dm_init.py:811
shesha.init.dm_init.makePetalDm
def makePetalDm(p_dm, p_geom, pupAngleDegree)
makePetalDm(p_dm, p_geom, pupAngleDegree)
Definition: dm_init.py:887
shesha.init.dm_init.make_kl_dm
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.
Definition: dm_init.py:691
shesha.init.dm_init.init_custom_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.
Definition: dm_init.py:544