COMPASS  5.4.4
End-to-end AO simulation tool using GPU acceleration
dm_util.py
1 
37 
38 import numpy as np
39 from astropy.io import fits
40 
41 import shesha.constants as scons
42 from shesha.constants import CONST
43 
44 from . import utilities as util
45 
46 from typing import List, Union
47 
48 
49 def dim_dm_support(cent: float, extent: int, ssize: int):
50  """ Compute the DM support dimensions
51 
52  Args:
53 
54  cent : (float): center of the pupil
55 
56  extent: (float): size of the DM support
57 
58  ssize: (int): size of ipupil support
59  """
60  n1 = np.floor(cent - extent / 2)
61  n2 = np.ceil(cent + extent / 2)
62  if (n1 < 1):
63  n1 = 1
64  if (n2 > ssize):
65  n2 = ssize
66 
67  return int(n1), int(n2)
68 
69 
70 def dim_dm_patch(pupdiam: int, diam: float, type: bytes, alt: float,
71  xpos_wfs: List[float], ypos_wfs: List[float]):
72  """ compute patchDiam for DM
73 
74  Args:
75 
76  pupdiam: (int) : pupil diameter
77 
78  diam: (float) : telescope diameter
79 
80  type: (bytes) : type of dm
81 
82  alt: (float) : altitude of dm
83 
84  xpos_wfs: (list) : list of wfs xpos
85 
86  ypos_wfs: (list) : list of wfs ypos
87  """
88 
89  if len(xpos_wfs) == 0:
90  norms = 0. # type: Union[float, List[float]]
91  else:
92  norms = [
93  np.linalg.norm([xpos_wfs[w], ypos_wfs[w]]) for w in range(len(xpos_wfs))
94  ]
95  if ((type == scons.DmType.PZT) or (type == scons.DmType.TT)):
96  pp = (diam / pupdiam)
97  elif (type == scons.DmType.KL):
98  pp = (pupdiam)
99  else:
100  raise TypeError("This type of DM doesn't exist ")
101 
102  patchDiam = int(pupdiam + 2 * np.max(norms) * CONST.ARCSEC2RAD * np.abs(alt) / (pp))
103  return patchDiam
104 
105 
106 def createSquarePattern(pitch: float, nxact: int):
107  """
108  Creates a list of M=nxact^2 actuator positions spread over an square grid.
109  Coordinates are centred around (0,0).
110 
111  Args:
112 
113  pitch: (float) : distance in pixels between 2 adjacent actus
114 
115  nxact: (int) : number of actu across the pupil diameter
116 
117  :return:
118 
119  xy: (np.ndarray(dims=2,dtype=np.float32)) : xy[M,2] list of coodinates
120  """
121 
122  xy = np.tile(np.arange(nxact) - (nxact - 1.) / 2., (nxact, 1)).astype(np.float32)
123  xy = np.array([xy.flatten(), xy.T.flatten()]) * pitch
124  xy = np.float32(xy)
125  return xy
126 
127 
128 def createHexaPattern(pitch: float, supportSize: int):
129  """
130  Creates a list of M actuator positions spread over an hexagonal grid.
131  The number M is the number of points of this grid, it cannot be
132  known before the procedure is called.
133  Coordinates are centred around (0,0).
134  The support that limits the grid is a square [-supportSize/2, supportSize/2].
135 
136  Args:
137 
138  pitch: (float) : distance in pixels between 2 adjacent actus
139 
140  supportSize: (int) : size in pixels of the support over which the coordinate list
141  should be returned.
142 
143  :return:
144 
145  xy: (np.ndarray(dims=2,dtype=np.float32)) : xy[2,M] list of coordinates
146  """
147  V3 = np.sqrt(3)
148  nx = int(np.ceil((supportSize / 2.0) / pitch) + 1)
149  x = pitch * (np.arange(2 * nx + 1, dtype=np.float32) - nx)
150  Nx = x.shape[0]
151  ny = int(np.ceil((supportSize / 2.0) / pitch / V3) + 1)
152  y = (V3 * pitch) * (np.arange(2 * ny + 1, dtype=np.float32) - ny)
153  Ny = y.shape[0]
154  x = np.tile(x, (Ny, 1)).flatten()
155  y = np.tile(y, (Nx, 1)).T.flatten()
156  x = np.append(x, x + pitch / 2.)
157  y = np.append(y, y + pitch * V3 / 2.)
158  xy = np.float32(np.array([y, x]))
159  return xy
160 
161 
162 def createDoubleHexaPattern(pitch: float, supportSize: int, pupAngleDegree: float):
163  """
164  Creates a list of M actuator positions spread over an hexagonal grid.
165  The number M is the number of points of this grid, it cannot be
166  known before the procedure is called.
167  Coordinates are centred around (0,0).
168  The support of the grid is a square [-supportSize/2,vsupportSize/2].
169 
170  Args:
171 
172  pitch: (float) : distance in pixels between 2 adjacent actus
173  supportSize: (int) : size in pixels of the support over which the coordinate list
174  should be returned.
175  pupAngleDegree: (float) : Rotation angle of the DM
176 
177  :return:
178 
179  xy: (np.ndarray(dims=2,dtype=np.float32)) : xy[2,M] list of coodinates
180  """
181  # ici on cree le pattern hexa. Par simplicite on replique le meme code
182  # que dans createHexaPattern(). On a donc un reseau "pointe selon X"
183  V3 = np.sqrt(3)
184  pi = np.pi
185  ny = int(np.ceil((supportSize / 2.0) / pitch) + 1)
186  y = pitch * (np.arange(2 * ny + 1, dtype=np.float32) - ny)
187  nx = int(np.ceil((supportSize / 2.0) / pitch / V3) + 1)
188  x = (V3 * pitch) * (np.arange(2 * nx + 1, dtype=np.float32) - nx)
189  # LA ligne de code qui change tout: on shifte le reseau de 1 pitch en X.
190  x = x + pitch
191  x, y = np.meshgrid(x, y, indexing='ij')
192  x = x.flatten()
193  y = y.flatten()
194  x = np.append(x, x + pitch * V3 / 2.)
195  y = np.append(y, y + pitch / 2.)
196 
197  # classement dans l'ordre kivabien
198  u = x + 1e-3 * y
199  idx = np.argsort(u)
200  x = x[idx]
201  y = y[idx]
202 
203  # on selection 1/6ieme du reseau, entre -30 et 30 degres
204  th = np.arctan2(y, x)
205  nn = np.where(((th < pi / 6) & (th > -pi / 6)))
206  x = x[nn]
207  y = y[nn]
208 
209  # on va maintenant repliquer ce reseau 6 fois, et le rotationnant a chaque
210  # fois de 60°. Note:
211  X = np.array([])
212  Y = np.array([])
213  for k in range(6):
214  xx = np.cos(k * pi / 3) * x + np.sin(k * pi / 3) * y
215  yy = -np.sin(k * pi / 3) * x + np.cos(k * pi / 3) * y
216  X = np.r_[X, xx]
217  Y = np.r_[Y, yy]
218 
219  # Rotation matrices pour suivre l'angle pupille
220  rot = pupAngleDegree * np.pi / 180.0
221  mrot = np.array([[np.cos(rot), -np.sin(rot)], [np.sin(rot), np.cos(rot)]])
222  XY = np.dot(mrot, [X, Y])
223  return np.float32(XY)
224 
225 
226 def filterActuWithPupil(actuPos: np.ndarray, pupil: np.ndarray,
227  threshold: float) -> np.ndarray:
228  '''
229  Select actuators based on their distance to the nearest pupil pixel
230  The implementation proposed here is easy but limits it precision to
231  an integer roundoff of the threshold
232 
233  actuPos: 2 x nActu np.array[float]: actuator position list - pupil pixel units
234  pupil: nPup x nPup np.ndarray[bool]: pupil mask
235  threshold: float: max allowed distance - pupil pixel units
236  '''
237 
238  # Gen a dilation mask of expected size
239  from scipy.ndimage.morphology import binary_dilation
240  k = np.ceil(threshold)
241  i, j = np.meshgrid(np.arange(2 * k + 1), np.arange(2 * k + 1), indexing='ij')
242  disk = ((i - k)**2 + (j - k)**2)**.5 <= k
243  dilatedPupil = binary_dilation(pupil, disk)
244  actuIsIn = dilatedPupil[(np.round(actuPos[0]).astype(np.int32),
245  np.round(actuPos[1]).astype(np.int32))]
246  return actuPos[:, actuIsIn]
247 
248 
249 def select_actuators(xc: np.ndarray, yc: np.ndarray, nxact: int, pitch: int, cobs: float,
250  margin_in: float, margin_out: float, N=None):
251  """
252  Select the "valid" actuators according to the system geometry
253 
254  Args:
255 
256  xc: actuators x positions (origine in center of mirror)
257 
258  yc: actuators y positions (origine in center of mirror)
259 
260  nxact:
261 
262  pitch:
263 
264  cobs:
265 
266  margin_in:
267 
268  margin_out:
269 
270  N:
271 
272  :return:
273 
274  liste_fin: actuator indice selection for xpos/ypos
275 
276 
277  """
278  # the following determine if an actuator is to be considered or not
279  # relative to the pitchmargin parameter.
280  dis = np.sqrt(xc**2 + yc**2)
281 
282  # test Margin_in
283  rad_in = (((nxact - 1) / 2) * cobs - margin_in) * pitch
284 
285  if N is None:
286  if (margin_out is None):
287  margin_out = 1.44
288  rad_out = ((nxact - 1.) / 2. + margin_out) * pitch
289 
290  valid_actus = np.where((dis <= rad_out) * (dis >= rad_in))[0]
291 
292  else:
293  valid_actus = np.where(dis >= rad_in)[0]
294  indsort = np.argsort(dis[valid_actus])
295 
296  if (N > valid_actus.size):
297  print('Too many actuators wanted, restricted to ', valid_actus.size)
298  else:
299  valid_actus = np.sort(indsort[:N])
300 
301  return valid_actus
302 
303 
304 def make_zernike(nzer: int, size: int, diameter: int, xc=-1., yc=-1., ext=0):
305  """Compute the zernike modes
306 
307  Args:
308 
309  nzer: (int) : number of modes
310 
311  size: (int) : size of the screen
312 
313  diameter: (int) : pupil diameter
314 
315  xc: (float) : (optional) x-position of the center
316 
317  yc: (float) : (optional) y-position of the center
318 
319  ext: (int) : (optional) extension
320 
321  :return:
322 
323  z : (np.ndarray(ndims=3,dtype=np.float64)) : zernikes modes
324  """
325  m = 0
326  n = 0
327 
328  if (xc == -1):
329  xc = size / 2
330  if (yc == -1):
331  yc = size / 2
332 
333  radius = (diameter) / 2.
334  zr = util.dist(size, xc, yc).astype(np.float32).T / radius
335  zmask = np.zeros((zr.shape[0], zr.shape[1], nzer), dtype=np.float32)
336  zmaskmod = np.zeros((zr.shape[0], zr.shape[1], nzer), dtype=np.float32)
337 
338  zmask[:, :, 0] = (zr <= 1).astype(np.float32)
339  zmaskmod[:, :, 0] = (zr <= 1.2).astype(np.float32)
340 
341  for i in range(1, nzer):
342  zmask[:, :, i] = zmask[:, :, 0]
343  zmaskmod[:, :, i] = zmaskmod[:, :, 0]
344 
345  zrmod = zr * zmaskmod[:, :, 0]
346 
347  zr = zr * zmask[:, :, 0]
348 
349  x = np.tile(np.arange(size).astype(np.float32), (size, 1))
350  zteta = np.arctan2(x - yc, x.T - xc).astype(np.float32)
351 
352  z = np.zeros((size, size, nzer), dtype=np.float32)
353 
354  for zn in range(nzer):
355  n, m = zernumero(zn + 1)
356 
357  if ext:
358  for i in range((n - m) // 2 + 1):
359  z[:, :, zn] = z[:, :, zn] + (-1.) ** i * zrmod ** (n - 2. * i) * float(np.math.factorial(n - i)) / \
360  float(np.math.factorial(i) * np.math.factorial((n + m) // 2 - i) *
361  np.math.factorial((n - m) // 2 - i))
362  else:
363  for i in range((n - m) // 2 + 1):
364  z[:, :, zn] = z[:, :, zn] + (-1.) ** i * zr ** (n - 2. * i) * float(np.math.factorial(n - i)) / \
365  float(np.math.factorial(i) * np.math.factorial((n + m) // 2 - i) *
366  np.math.factorial((n - m) // 2 - i))
367 
368  if ((zn + 1) % 2 == 1):
369  if (m == 0):
370  z[:, :, zn] = z[:, :, zn] * np.sqrt(n + 1.)
371  else:
372  z[:, :, zn] = z[:, :, zn] * \
373  np.sqrt(2. * (n + 1)) * np.sin(m * zteta)
374  else:
375  if (m == 0):
376  z[:, :, zn] = z[:, :, zn] * np.sqrt(n + 1.)
377  else:
378  z[:, :, zn] = z[:, :, zn] * \
379  np.sqrt(2. * (n + 1)) * np.cos(m * zteta)
380 
381  if (ext):
382  return z * zmaskmod
383  else:
384  return z * zmask
385 
386 
387 def zernumero(zn: int):
388  """
389  Returns the radial degree and the azimuthal number of zernike
390  number zn, according to Noll numbering (Noll, JOSA, 1976)
391 
392  Args:
393 
394  zn: (int) : zernike number
395 
396  :returns:
397 
398  rd: (int) : radial degrees
399 
400  an: (int) : azimuthal numbers
401 
402  """
403  j = 0
404  for n in range(101):
405  for m in range(n + 1):
406  if ((n - m) % 2 == 0):
407  j = j + 1
408  if (j == zn):
409  return n, m
410  if (m != 0):
411  j = j + 1
412  if (j == zn):
413  return n, m
414 
415 
416 dm_fits_content="""The DM FITS file is compatible with COMPASS DM database.
417  The primary header contains the keywords:
418  * PIXSIZE : the size of the pixels on the maps in meters.
419 
420  * XCENTER, YCENTER are the coordinates of the centre of the pupil, expressed in pixels, in a reference frame conformable to (i,j) coords. The translation from pixels to meters can be done using:
421  meters = (pixels - XCENTER) * PIXSIZE
422 
423  * PUPM is the diameter of pupil stop (meters).
424 
425  * Additionally the header provides the user with:
426  PITCHM is the size of the DM pitch in meters (may be handy is some cases and useful when using this file with COMPASS software)
427 
428  This FITS file contains 3 extensions:
429  * Extension 'I1_J1' are the coordinate (i,j) of the first pixel for each of the 2D maps (see Extension 2), so that they can be inserted in a larger map.
430 
431  * Extension 'INFLU' are the 2D maps of the influence functions.
432 
433  * Extension 'XPOS_YPOS' are the coordinates (xpos, ypos) of the physical location of the actuator, in pixels. This data is provided for information only and does not directly participate to build the DM. The present coordinates are positions in M1 space, i.e. include the distorsion due to telescope optics.
434 """
435 
436 def add_doc_content(*content):
437  """adds content to a docstring (to be used as decorator)"""
438  def dec(obj):
439  obj.__doc__ = obj.__doc__.format(content)
440  return obj
441  return dec
442 
443 def write_dm_custom_fits(file_name, i1, j1, influ_cube, xpos, ypos, xcenter, ycenter,pixsize,pupm, *,pitchm=None):
444  """Write a custom_dm fits file based on user provided data (see args)
445 
446  Args:
447  file_name : (string) : name of the custom dm fits file
448 
449  i1 : (np.ndarray) : x coordinate of the first pixel for each of the 2D maps
450 
451  j1 : (np.ndarray) : y coordinate of the first pixel for each of the 2D maps
452 
453  influ_cube : (np.ndarray) : 2D maps of the influence functions
454 
455  xpos : (np.ndarray) : x coordinate of the physical location of the actuator
456 
457  ypos : (np.ndarray) : y coordinate of the physical location of the actuator
458 
459  xcenter : (float) : x coordinate of the centre of the pupil, expressed in pixels
460 
461  ycenter : (float) : y coordinate of the centre of the pupil, expressed in pixels
462 
463  pixsize : (float) : size of the pixels on the maps in meters
464 
465  pupm : (float) : diameter of pupil stop (meters)
466 
467  Kwargs:
468  pitchm : (float) : size of the DM pitch in meters. Defaults to None.
469 
470  Returns:
471  (HDUList) : custom_dm data
472  """
473  fits_version=1.2
474  primary_hdu = fits.PrimaryHDU()
475  primary_hdu.header['VERSION'] = (fits_version,'file format version')
476  primary_hdu.header['XCENTER'] = (xcenter ,'DM centre along X in pixels')
477  primary_hdu.header['YCENTER'] = (ycenter ,'DM centre along Y in pixels')
478  primary_hdu.header['PIXSIZE'] = (pixsize ,'pixel size (meters)')
479  primary_hdu.header['PUPM'] = (pupm ,'nominal pupil diameter (meters)')
480  if(pitchm is not None):
481  primary_hdu.header['PITCHM'] = (pitchm,'DM pitch (meters)')
482 
483  for line in dm_fits_content.splitlines():
484  primary_hdu.header.add_comment(line)
485 
486  image_hdu = fits.ImageHDU(np.c_[i1 , j1 ].T, name="I1_J1")
487  image_hdu2 = fits.ImageHDU(influ_cube, name="INFLU")
488  image_hdu3 = fits.ImageHDU(np.c_[xpos, ypos].T, name="XPOS_YPOS")
489 
490  dm_custom = fits.HDUList([primary_hdu, image_hdu, image_hdu2, image_hdu3])
491 
492  dm_custom.writeto(file_name,overwrite=1)
493  return dm_custom
494 
495 @add_doc_content(dm_fits_content)
496 def export_custom_dm(file_name, p_dm, p_geom, *, p_tel=None):
497  """Return an HDUList (FITS) with the data required to create a COMPASS custom_dm
498 
499  {}
500 
501  Args:
502  p_dm : (Param_dm) : dm settings
503 
504  p_geom : (Param_geom) : geometry settings
505 
506  Kwargs:
507  file_name : (string) : if set, the HDU is written to the file specified by this variable
508 
509  p_tel : (Param_tel) : telescope settings, used to provide the diameter (if not provided, the default diameter id obtained from the p_geom as pupdiam*pixsize)
510 
511  Returns:
512  (HDUList) : custom_dm data
513  """
514 
515  pixsize = p_geom.get_pixsize()
516  diam = p_geom.pupdiam * p_geom._pixsize
517  if(p_tel is not None):
518  diam = p_tel.diam
519  xpos = p_dm._xpos
520  ypos = p_dm._ypos
521  i1 = p_dm._i1 + p_dm._n1
522  j1 = p_dm._j1 + p_dm._n1
523  influ = p_dm._influ / p_dm.unitpervolt
524 
525  xcenter = p_geom.cent
526  ycenter = p_geom.cent
527 
528  pitchm=None
529  if p_dm._pitch :
530  pitchm = p_dm._pitch*pixsize
531  dm_custom = write_dm_custom_fits(file_name,i1,j1,influ,xpos,ypos,xcenter,ycenter,pixsize,diam,pitchm=pitchm)
532 
533  return dm_custom
Numerical constants for shesha and config enumerations for safe-typing.
Definition: constants.py:1
def createDoubleHexaPattern(float pitch, int supportSize, float pupAngleDegree)
Creates a list of M actuator positions spread over an hexagonal grid.
Definition: dm_util.py:180
def dim_dm_support(float cent, int extent, int ssize)
Compute the DM support dimensions.
Definition: dm_util.py:59
def select_actuators(np.ndarray xc, np.ndarray yc, int nxact, int pitch, float cobs, float margin_in, float margin_out, N=None)
Select the "valid" actuators according to the system geometry.
Definition: dm_util.py:277
def export_custom_dm(file_name, p_dm, p_geom, *p_tel=None)
Return an HDUList (FITS) with the data required to create a COMPASS custom_dm.
Definition: dm_util.py:515
def write_dm_custom_fits(file_name, i1, j1, influ_cube, xpos, ypos, xcenter, ycenter, pixsize, pupm, *pitchm=None)
Write a custom_dm fits file based on user provided data (see args)
Definition: dm_util.py:472
def zernumero(int zn)
Returns the radial degree and the azimuthal number of zernike number zn, according to Noll numbering ...
Definition: dm_util.py:402
def make_zernike(int nzer, int size, int diameter, xc=-1., yc=-1., ext=0)
Compute the zernike modes.
Definition: dm_util.py:324
def add_doc_content(*content)
adds content to a docstring (to be used as decorator)
Definition: dm_util.py:437
np.ndarray filterActuWithPupil(np.ndarray actuPos, np.ndarray pupil, float threshold)
Select actuators based on their distance to the nearest pupil pixel The implementation proposed here ...
Definition: dm_util.py:236
def dim_dm_patch(int pupdiam, float diam, bytes type, float alt, List[float] xpos_wfs, List[float] ypos_wfs)
compute patchDiam for DM
Definition: dm_util.py:87
def createHexaPattern(float pitch, int supportSize)
Creates a list of M actuator positions spread over an hexagonal grid.
Definition: dm_util.py:146
def createSquarePattern(float pitch, int nxact)
Creates a list of M=nxact^2 actuator positions spread over an square grid.
Definition: dm_util.py:120