COMPASS  5.0.0
End-to-end AO simulation tool using GPU acceleration
tomo.py
1 
37 
38 import numpy as np
39 
40 import shesha.config as conf
41 import shesha.constants as scons
42 from shesha.constants import CONST
43 from shesha.sutra_wrap import Sensors, Dms, Rtc_FFF as Rtc, Atmos
44 
45 import typing
46 from typing import List
47 
48 
49 def do_tomo_matrices(ncontrol: int, rtc: Rtc, p_wfss: List[conf.Param_wfs], dms: Dms,
50  atmos: Atmos, wfs: Sensors, p_controller: conf.Param_controller,
51  p_geom: conf.Param_geom, p_dms: list, p_tel: conf.Param_tel,
52  p_atmos: conf.Param_atmos):
53  """ Compute Cmm and Cphim matrices for the MV controller on GPU
54 
55  :parameters:
56 
57  ncontrol: (int): controller index
58 
59  rtc: (Rtc) : rtc object
60 
61  p_wfss: (list of Param_wfs) : wfs settings
62 
63  dms: (Dms) : Dms object
64 
65  atmos: (Atmos) : Atmos object
66 
67  wfs: (Sensors) : Sensors object
68 
69  p_controller: (Param_controller): controller settings
70 
71  p_geom: (Param_geom) : geom settings
72 
73  p_dms: (list of Param_dms) : dms settings
74 
75  p_tel: (Param_tel) : telescope settings
76 
77  p_atmos: (Param_atmos) : atmos settings
78  """
79  nvalidperwfs = np.array([o._nvalid for o in p_wfss], dtype=np.int64)
80  # Bring bottom left corner of valid subapertures in ipupil
81  ipup = p_geom._ipupil
82  spup = p_geom._spupil
83  s2ipup = (ipup.shape[0] - spup.shape[0]) / 2.
84  # Total number of subapertures
85  nvalid = sum([nvalidperwfs[j] for j in p_controller.nwfs])
86  ind = 0
87  # X-position of the bottom left corner of each valid subaperture
88  X = np.zeros(nvalid, dtype=np.float64)
89  # Y-position of the bottom left corner of each subaperture
90  Y = np.zeros(nvalid, dtype=np.float64)
91 
92  for k in p_controller.nwfs:
93  posx = p_wfss[k]._validpuppixx + s2ipup
94  # X-position of the bottom left corner of each valid subaperture
95  # bring the origin in the center of ipupil and 0-index it
96  posx = posx - ipup.shape[0] / 2 - 1
97  posy = p_wfss[k]._validpuppixy + s2ipup
98  posy = posy.T - ipup.shape[0] / 2 - 1
99  p2m = (p_tel.diam / p_wfss[k].nxsub) / \
100  p_wfss[k]._pdiam # Size of one pixel in meters
101  posx *= p2m # Positions in meters
102  posy *= p2m
103  X[ind:ind + p_wfss[k]._nvalid] = posx
104  Y[ind:ind + p_wfss[k]._nvalid] = posy
105  ind += p_wfss[k]._nvalid
106 
107  # Get the total number of pzt DM and actuators to control
108  nactu = 0
109  npzt = 0
110  for k in p_controller.ndm:
111  if (p_dms[k].type == scons.DmType.PZT):
112  nactu += p_dms[k]._ntotact
113  npzt += 1
114 
115  Xactu = np.zeros(nactu, dtype=np.float64) # X-position actuators in ipupil
116  Yactu = np.zeros(nactu, dtype=np.float64) # Y-position actuators in ipupil
117  k2 = np.zeros(npzt, dtype=np.float64) # k2 factors for computation
118  pitch = np.zeros(npzt, dtype=np.float64)
119  alt_DM = np.zeros(npzt, dtype=np.float64)
120  ind = 0
121  indk = 0
122  for k in p_controller.ndm:
123  if (p_dms[k].type == scons.DmType.PZT):
124  p2m = p_tel.diam / p_geom.pupdiam
125  # Conversion in meters in the center of ipupil
126  actu_x = (p_dms[k]._xpos - ipup.shape[0] / 2) * p2m
127  actu_y = (p_dms[k]._ypos - ipup.shape[0] / 2) * p2m
128  pitch[indk] = actu_x[1] - actu_x[0]
129  k2[indk] = p_wfss[0].Lambda / 2. / np.pi / p_dms[k].unitpervolt
130  alt_DM[indk] = p_dms[k].alt
131  Xactu[ind:ind + p_dms[k]._ntotact] = actu_x
132  Yactu[ind:ind + p_dms[k]._ntotact] = actu_y
133 
134  ind += p_dms[k]._ntotact
135  indk += 1
136 
137  # Select a DM for each layer of atmos
138  NlayersDM = np.zeros(npzt, dtype=np.int64) # Useless for now
139  indlayersDM = selectDMforLayers(p_atmos, p_controller, p_dms)
140  # print("indlayer = ",indlayersDM)
141 
142  # Get FoV
143  # conf.RAD2ARCSEC = 180.0/np.pi * 3600.
144  wfs_distance = np.zeros(len(p_controller.nwfs), dtype=np.float64)
145  ind = 0
146  for k in p_controller.nwfs:
147  wfs_distance[ind] = np.sqrt(p_wfss[k].xpos**2 + p_wfss[k].ypos**2)
148  ind += 1
149  FoV = np.max(wfs_distance) / CONST.RAD2ARCSEC
150 
151  # WFS postions in rad
152  alphaX = np.zeros(len(p_controller.nwfs))
153  alphaY = np.zeros(len(p_controller.nwfs))
154 
155  ind = 0
156  for k in p_controller.nwfs:
157  alphaX[ind] = p_wfss[k].xpos / CONST.RAD2ARCSEC
158  alphaY[ind] = p_wfss[k].ypos / CONST.RAD2ARCSEC
159  ind += 1
160 
161  L0_d = np.copy(p_atmos.L0).astype(np.float64)
162  frac_d = np.copy(p_atmos.frac * (p_atmos.r0**(-5.0 / 3.0))).astype(np.float64)
163 
164  print("Computing Cphim...")
165  rtc.d_control[ncontrol].compute_Cphim(atmos, wfs, dms, L0_d, frac_d, alphaX, alphaY,
166  X, Y, Xactu, Yactu, p_tel.diam, k2, NlayersDM,
167  indlayersDM, FoV, pitch,
168  alt_DM.astype(np.float64))
169  print("Done")
170 
171  print("Computing Cmm...")
172  rtc.d_control[ncontrol].compute_Cmm(atmos, wfs, L0_d, frac_d, alphaX, alphaY,
173  p_tel.diam, p_tel.cobs)
174  print("Done")
175 
176  Nact = np.zeros([nactu, nactu], dtype=np.float32)
177  F = np.zeros([nactu, nactu], dtype=np.float32)
178  ind = 0
179  for k in range(len(p_controller.ndm)):
180  if (p_dms[k].type == "pzt"):
181  Nact[ind:ind + p_dms[k]._ntotact, ind:ind +
182  p_dms[k]._ntotact] = create_nact_geom(p_dms[k])
183  F[ind:ind + p_dms[k]._ntotact, ind:ind +
184  p_dms[k]._ntotact] = create_piston_filter(p_dms[k])
185  ind += p_dms[k]._ntotact
186  rtc.d_control[ncontrol].filter_cphim(F, Nact)
187 
188 
189 def selectDMforLayers(p_atmos: conf.Param_atmos, p_controller: conf.Param_controller,
190  p_dms: list):
191  """ For each atmos layer, select the DM which have to handle it in the Cphim computation for MV controller
192 
193  :parameters:
194 
195  p_atmos : (Param_atmos) : atmos parameters
196 
197  p_controller : (Param_controller) : controller parameters
198 
199  p_dms :(list of Param_dm) : dms parameters
200 
201  :return:
202 
203  indlayersDM : (np.array(dtype=np.int32)) : for each atmos layer, the Dm number corresponding to it
204  """
205  indlayersDM = np.zeros(p_atmos.nscreens, dtype=np.int64)
206  for i in range(p_atmos.nscreens):
207  mindif = 1e6
208  for j in p_controller.ndm:
209  alt_diff = np.abs(p_dms[j].alt - p_atmos.alt[i])
210  if (alt_diff < mindif):
211  indlayersDM[i] = j
212  mindif = alt_diff
213 
214  return indlayersDM
215 
216 
217 def create_nact_geom(p_dm: conf.Param_dm):
218  """ Compute the DM coupling matrix
219 
220  :param:
221 
222  p_dm : (Param_dm) : dm parameters
223 
224  :return:
225 
226  Nact : (np.array(dtype=np.float64)) : the DM coupling matrix
227  """
228  nactu = p_dm._ntotact
229  Nact = np.zeros([nactu, nactu], dtype=np.float32)
230  coupling = p_dm.coupling
231  dim = p_dm._n2 - p_dm._n1 + 1
232  mask = np.zeros([dim, dim], dtype=np.float32)
233  shape = np.zeros([dim, dim], dtype=np.float32)
234 
235  for i in range(len(p_dm._i1)):
236  mask[p_dm._j1[i]][p_dm._i1[i]] = 1
237 
238  mask_act = np.where(mask)
239 
240  pitch = int(p_dm._pitch)
241 
242  for i in range(nactu):
243  shape *= 0
244  # Diagonal
245  shape[p_dm._j1[i]][p_dm._i1[i]] = 1
246  # Left, right, above and under the current actuator
247  shape[p_dm._j1[i]][p_dm._i1[i] - pitch] = coupling
248  shape[p_dm._j1[i] - pitch][p_dm._i1[i]] = coupling
249  shape[p_dm._j1[i]][p_dm._i1[i] + pitch] = coupling
250  shape[p_dm._j1[i] + pitch][p_dm._i1[i]] = coupling
251  # Diagonals of the current actuators
252  shape[p_dm._j1[i] - pitch][p_dm._i1[i] - pitch] = coupling**2
253  shape[p_dm._j1[i] - pitch][p_dm._i1[i] + pitch] = coupling**2
254  shape[p_dm._j1[i] + pitch][p_dm._i1[i] + pitch] = coupling**2
255  shape[p_dm._j1[i] + pitch][p_dm._i1[i] - pitch] = coupling**2
256 
257  Nact[:, i] = shape[mask_act]
258 
259  return Nact
260 
261 
262 def create_piston_filter(p_dm: conf.Param_dm):
263  """ Create the piston filter matrix
264 
265  :parameters:
266 
267  p_dm: (Param_dm): dm settings
268  """
269  nactu = p_dm._ntotact
270  F = np.ones([nactu, nactu], dtype=np.float32)
271  F = F * (-1.0 / nactu)
272  for i in range(nactu):
273  F[i][i] = 1 - 1.0 / nactu
274  return F
shesha.ao.tomo.create_piston_filter
def create_piston_filter(conf.Param_dm p_dm)
Create the piston filter matrix.
Definition: tomo.py:268
shesha.ao.tomo.create_nact_geom
def create_nact_geom(conf.Param_dm p_dm)
Compute the DM coupling matrix.
Definition: tomo.py:227
shesha.sutra_wrap
Definition: sutra_wrap.py:1
shesha.constants
Numerical constants for shesha and config enumerations for safe-typing.
Definition: constants.py:1
shesha.ao.tomo.do_tomo_matrices
def do_tomo_matrices(int ncontrol, Rtc rtc, List[conf.Param_wfs] p_wfss, Dms dms, Atmos atmos, Sensors wfs, conf.Param_controller p_controller, conf.Param_geom p_geom, list p_dms, conf.Param_tel p_tel, conf.Param_atmos p_atmos)
Compute Cmm and Cphim matrices for the MV controller on GPU.
Definition: tomo.py:75
shesha.ao.tomo.selectDMforLayers
def selectDMforLayers(conf.Param_atmos p_atmos, conf.Param_controller p_controller, list p_dms)
For each atmos layer, select the DM which have to handle it in the Cphim computation for MV controlle...
Definition: tomo.py:203
shesha.config
Parameter classes for COMPASS.
Definition: shesha/shesha/config/__init__.py:1