COMPASS  5.0.0
End-to-end AO simulation tool using GPU acceleration
groot.py
1 """
2 GROOT (Gpu-based Residual errOr cOvariance maTrix)
3 Python module for modelization of error covariance matrix
4 """
5 import numpy as np
6 import h5py
7 from shesha.sutra_wrap import carmaWrap_context, Groot
8 import time
9 import sys
10 import os
11 from tqdm import tqdm
12 
13 #from guardians import gamora
14 from guardians import drax, starlord
15 import matplotlib.pyplot as plt
16 plt.ion()
17 
18 gpudevices = np.array([0, 1, 2, 3], dtype=np.int32)
19 cxt = carmaWrap_context.get_instance_ngpu(gpudevices.size, gpudevices)
20 
21 
22 def compute_Cerr(filename, modal=True, ctype="float", speed=None, H=None, theta=None,
23  r0=None, L0=None):
24  """ Returns the residual error covariance matrix using GROOT from a ROKET file
25  :parameter:
26  filename : (string) : full path to the ROKET file
27  modal : (bool) : if True (default), Cerr is returned in the Btt modal basis,
28  in the actuator basis if False
29  ctype : (string) : "float" or "double"
30  speed: (np.ndarray(ndim=1, dtype=np.float32)): (optionnal) wind speed for each layer [m/s]
31  H: (np.ndarray(ndim=1, dtype=np.float32)): (optionnal) altitude of each layer [m]
32  theta: (np.ndarray(ndim=1, dtype=np.float32)): (optionnal) wind direction for each layer [rad]
33  r0: (float): (optionnal) Fried parameter @ 0.5 µm [m]
34  L0: (np.ndarray(ndim=1, dtype=np.float32)): (optionnal) Outer scale [m]
35 
36  :return:
37  Cerr : (np.ndarray(dim=2, dtype=np.float32)) : residual error covariance matrix
38  """
39  f = h5py.File(filename, 'r')
40  Lambda_tar = f.attrs["_Param_target__Lambda"][0]
41  Lambda_wfs = f.attrs["_Param_wfs__Lambda"]
42  dt = f.attrs["_Param_loop__ittime"]
43  gain = f.attrs["_Param_controller__gain"]
44  wxpos = f.attrs["_Param_wfs__xpos"][0]
45  wypos = f.attrs["_Param_wfs__ypos"][0]
46  if r0 is None:
47  r0 = f.attrs["_Param_atmos__r0"]
48  r0 = r0 * (Lambda_tar / Lambda_wfs)**(6. / 5.)
49  RASC = 180. / np.pi * 3600.
50  xpos = f["dm.xpos"][:]
51  ypos = f["dm.ypos"][:]
52  p2m = f.attrs["_Param_tel__diam"] / f.attrs["_Param_geom__pupdiam"]
53  pupshape = int(2**np.ceil(np.log2(f.attrs["_Param_geom__pupdiam"]) + 1))
54  xactu = (xpos - pupshape / 2) * p2m
55  yactu = (ypos - pupshape / 2) * p2m
56  if H is None:
57  H = f.attrs["_Param_atmos__alt"]
58  if L0 is None:
59  L0 = f.attrs["_Param_atmos__L0"]
60  if speed is None:
61  speed = f.attrs["_Param_atmos__windspeed"]
62  if theta is None:
63  theta = (f.attrs["_Param_atmos__winddir"] * np.pi / 180.)
64  frac = f.attrs["_Param_atmos__frac"]
65 
66  Htheta = np.linalg.norm([wxpos, wypos]) / RASC * H
67  vdt = speed * dt / gain
68  angleht = np.arctan2(wypos, wxpos)
69  fc = 1 / (2 * (xactu[1] - xactu[0]))
70  scale = (1 / r0)**(5 / 3.) * frac * (Lambda_tar / (2 * np.pi))**2
71  Nact = f["Nact"][:]
72  Nact = np.linalg.inv(Nact)
73  P = f["P"][:]
74  Btt = f["Btt"][:]
75  Tf = Btt[:-2, :-2].dot(P[:-2, :-2])
76  IF, T = drax.get_IF(filename)
77  IF = IF.T
78  T = T.T
79  N = IF.shape[0]
80  deltaTT = T.T.dot(T) / N
81  deltaF = IF.T.dot(T) / N
82  pzt2tt = np.linalg.inv(deltaTT).dot(deltaF.T)
83 
84  if (ctype == "float"):
85  groot = Groot(cxt, cxt.active_device, Nact.shape[0],
86  int(f.attrs["_Param_atmos__nscreens"]), angleht,
87  vdt.astype(np.float32), Htheta.astype(np.float32), L0, theta,
88  scale.astype(np.float32), pzt2tt.astype(np.float32),
89  Tf.astype(np.float32), Nact.astype(np.float32),
90  xactu.astype(np.float32), yactu.astype(np.float32), fc)
91  else:
92  raise TypeError("Unknown ctype : must be float")
93  tic = time.time()
94  groot.compute_Cerr()
95  Cerr = np.array(groot.d_Cerr)
96  cov_err_groot = np.zeros((Nact.shape[0] + 2, Nact.shape[0] + 2))
97  cov_err_groot[:-2, :-2] = Cerr
98  cov_err_groot[-2:, -2:] = np.array(groot.d_TT)
99  tac = time.time()
100  print("Cee computed in : %.2f seconds" % (tac - tic))
101  if (modal):
102  cov_err_groot = P.dot(cov_err_groot).dot(P.T)
103 
104  f.close()
105  return cov_err_groot
106 
107 
108 def compute_Cerr_cpu(filename, modal=True):
109  """ Returns the residual error covariance matrix using CPU version of GROOT
110  from a ROKET file
111  :parameter:
112  filename : (string) : full path to the ROKET file
113  modal : (bool) : if True (default), Cerr is returned in the Btt modal basis,
114  in the actuator basis if False
115  :return:
116  Cerr : (np.ndarray(dim=2, dtype=np.float32)) : residual error covariance matrix
117  """
118  f = h5py.File(filename, 'r')
119 
120  tabx, taby = starlord.tabulateIj0()
121  Lambda_tar = f.attrs["_Param_target__Lambda"][0]
122  Lambda_wfs = f.attrs["_Param_wfs__Lambda"]
123  dt = f.attrs["_Param_loop__ittime"]
124  gain = f.attrs["_Param_controller__gain"]
125  wxpos = f.attrs["_Param_wfs__xpos"][0]
126  wypos = f.attrs["_Param_wfs__ypos"][0]
127  r0 = f.attrs["_Param_atmos__r0"] * (Lambda_tar / Lambda_wfs)**(6. / 5.)
128  RASC = 180. / np.pi * 3600.
129  xpos = f["dm.xpos"][:]
130  ypos = f["dm.ypos"][:]
131  p2m = f.attrs["_Param_tel__diam"] / f.attrs["_Param_geom__pupdiam"]
132  pupshape = int(2**np.ceil(np.log2(f.attrs["_Param_geom__pupdiam"]) + 1))
133  xactu = (xpos - pupshape / 2) * p2m
134  yactu = (ypos - pupshape / 2) * p2m
135  Ccov = np.zeros((xpos.size, xpos.size))
136  Caniso = np.zeros((xpos.size, xpos.size))
137  Cbp = np.zeros((xpos.size, xpos.size))
138  xx = np.tile(xactu, (xactu.shape[0], 1))
139  yy = np.tile(yactu, (yactu.shape[0], 1))
140  xij = xx - xx.T
141  yij = yy - yy.T
142 
143  for l in range(f.attrs["_Param_atmos__nscreens"]):
144  H = f.attrs["_Param_atmos__alt"][l]
145  L0 = f.attrs["_Param_atmos__L0"][l]
146  speed = f.attrs["_Param_atmos__windspeed"][l]
147  theta = f.attrs["_Param_atmos__winddir"][l] * np.pi / 180.
148  frac = f.attrs["_Param_atmos__frac"][l]
149 
150  Htheta = np.linalg.norm([wxpos, wypos]) / RASC * H
151  vdt = speed * dt / gain
152  # Covariance matrices models on actuators space
153  M = np.zeros((xpos.size, xpos.size))
154  Mvdt = M.copy()
155  Mht = M.copy()
156  Mhvdt = M.copy()
157  angleht = np.arctan2(wypos, wxpos)
158  fc = xactu[1] - xactu[0]
159 
160  M = np.linalg.norm([xij, yij], axis=0)
161  Mvdt = np.linalg.norm([xij - vdt * np.cos(theta), yij - vdt * np.sin(theta)],
162  axis=0)
163  Mht = np.linalg.norm(
164  [xij - Htheta * np.cos(angleht), yij - Htheta * np.sin(angleht)], axis=0)
165  Mhvdt = np.linalg.norm([
166  xij - vdt * np.cos(theta) - Htheta * np.cos(angleht),
167  yij - vdt * np.sin(theta) - Htheta * np.sin(angleht)
168  ], axis=0)
169 
170  Ccov += 0.5 * (starlord.dphi_lowpass(Mhvdt, fc, L0, tabx, taby) -
171  starlord.dphi_lowpass(Mht, fc, L0, tabx, taby) - starlord.
172  dphi_lowpass(Mvdt, fc, L0, tabx, taby) + starlord.dphi_lowpass(
173  M, fc, L0, tabx, taby)) * (1. / r0)**(5. / 3.) * frac
174 
175  Caniso += 0.5 * (
176  starlord.dphi_lowpass(Mht, fc, L0, tabx, taby) - starlord.dphi_lowpass(
177  M, fc, L0, tabx, taby)) * (1. / r0)**(5. / 3.) * frac
178  Cbp += 0.5 * (starlord.dphi_lowpass(Mvdt, fc, L0, tabx, taby) - starlord.
179  dphi_lowpass(M, fc, L0, tabx, taby)) * (1. / r0)**(5. / 3.) * frac
180 
181  Sp = (Lambda_tar / (2 * np.pi))**2
182  Ctt = (Caniso + Caniso.T) * Sp
183  Ctt += ((Cbp + Cbp.T) * Sp)
184  Ctt += ((Ccov + Ccov.T) * Sp)
185 
186  P = f["P"][:]
187  Btt = f["Btt"][:]
188  Tf = Btt[:-2, :-2].dot(P[:-2, :-2])
189 
190  IF, T = drax.get_IF(filename)
191  IF = IF.T
192  T = T.T
193  N = IF.shape[0]
194  deltaTT = T.T.dot(T) / N
195  deltaF = IF.T.dot(T) / N
196  pzt2tt = np.linalg.inv(deltaTT).dot(deltaF.T)
197 
198  Nact = f["Nact"][:]
199  N1 = np.linalg.inv(Nact)
200  Ctt = N1.dot(Ctt).dot(N1)
201  ttcomp = pzt2tt.dot(Ctt).dot(pzt2tt.T)
202  Ctt = Tf.dot(Ctt).dot(Tf.T)
203  cov_err = np.zeros((Ctt.shape[0] + 2, Ctt.shape[0] + 2))
204  cov_err[:-2, :-2] = Ctt
205  cov_err[-2:, -2:] = ttcomp
206  if (modal):
207  cov_err = P.dot(cov_err).dot(P.T)
208  f.close()
209 
210  return cov_err
211 
212 
213 def test_Cerr(filename):
214  """ Compute PSF of aniso and bandwidth from GROOT model and ROKET to compare
215 
216  :parameters:
217  filename:(str):path to the ROKET file
218  """
219  C = drax.get_covmat_contrib(filename, ["bandwidth", "tomography"])
220  Cerr = compute_Cerr(filename)
221  _, _, psfr, _ = gamora.psf_rec_Vii(filename, covmodes=C.astype(np.float32),
222  fitting=False)
223  _, _, psf, _ = gamora.psf_rec_Vii(filename, cov=Cerr.astype(np.float32),
224  fitting=False)
225  drax.cutsPSF(filename, psfr, psf)
226  print("PSFR SR: ", psfr.max())
227  print("PSF SR: ", psf.max())
228  psf = drax.ensquare_PSF(filename, psf, 20)
229  psfr = drax.ensquare_PSF(filename, psfr, 20)
230  plt.matshow(np.log10(np.abs(psfr)))
231  plt.colorbar()
232  plt.title("PSF_R")
233  plt.matshow(
234  np.log10(np.abs(psf)), vmax=np.log10(np.abs(psfr)).max(),
235  vmin=np.log10(np.abs(psfr)).min())
236  plt.colorbar()
237  plt.title("PSF")
238  plt.matshow(
239  np.log10(np.abs(psfr - psf)), vmax=np.log10(np.abs(psfr)).max(),
240  vmin=np.log10(np.abs(psfr)).min())
241  plt.colorbar()
242  plt.title("PSF_R - PSF")
243 
244  return psf, psfr
245 
246 
247 def compare_GPU_vs_CPU(filename):
248  """ Compare results of GROOT vs its CPU version in terms of execution time
249  and precision on the PSF renconstruction
250  :parameter:
251  filename : (string) : full path to the ROKET file
252 
253  """
254  timer = ch.carmaWrap_timer()
255 
256  ch.threadSync()
257  timer.start()
258  ch.threadSync()
259  synctime = timer.stop()
260  timer.reset()
261 
262  timer.start()
263  cov_err_gpu_s = compute_Cerr(filename)
264  ch.threadSync()
265  gpu_time_s = timer.stop() - synctime
266  timer.reset()
267 
268  timer.start()
269  cov_err_gpu_d = compute_Cerr(filename, ctype="double")
270  ch.threadSync()
271  gpu_time_d = timer.stop() - synctime
272  timer.reset()
273 
274  tic = time.time()
275  cov_err_cpu = compute_Cerr_cpu(filename)
276  tac = time.time()
277  cpu_time = tac - tic
278 
279  otftel, otf2, psf_cpu, gpu = gamora.psf_rec_Vii(filename, fitting=False,
280  cov=cov_err_cpu.astype(np.float32))
281  otftel, otf2, psf_gpu_s, gpu = gamora.psf_rec_Vii(
282  filename, fitting=False, cov=cov_err_gpu_s.astype(np.float32))
283  otftel, otf2, psf_gpu_d, gpu = gamora.psf_rec_Vii(
284  filename, fitting=False, cov=cov_err_gpu_d.astype(np.float32))
285 
286  print("-----------------------------------------")
287  print("CPU time : ", cpu_time, " s ")
288  print("GPU time simple precision : ", gpu_time_s, " s ")
289  print("GPU time double precision : ", gpu_time_d, " s ")
290  print("Max absolute difference in PSFs simple precision : ",
291  np.abs(psf_cpu - psf_gpu_s).max())
292  print("Max absolute difference in PSFs double precision : ",
293  np.abs(psf_cpu - psf_gpu_d).max())
294  gamora.cutsPSF(filename, psf_cpu, psf_gpu_s)
295  gamora.cutsPSF(filename, psf_cpu, psf_gpu_d)
296 
297 
298 def compute_Ca_cpu(filename, modal=True):
299  """ Returns the aliasing error covariance matrix using CPU version of GROOT
300  from a ROKET file
301  :parameter:
302  filename : (string) : full path to the ROKET file
303  modal : (bool) : if True (default), Ca is returned in the Btt modal basis,
304  in the actuator basis if False
305  :return:
306  Ca : (np.ndarray(dim=2, dtype=np.float32)) : aliasing error covariance matrix
307  """
308  f = h5py.File(filename, 'r')
309  nsub = f["R"][:].shape[1] // 2
310  nssp = f.attrs["_Param_wfs__nxsub"][0]
311  validint = f.attrs["_Param_tel__cobs"]
312  x = np.linspace(-1, 1, nssp)
313  x, y = np.meshgrid(x, x)
314  r = np.sqrt(x * x + y * y)
315  rorder = np.sort(r.reshape(nssp * nssp))
316  ncentral = nssp * nssp - np.sum(r >= validint)
317  validext = rorder[ncentral + nsub]
318  valid = (r < validext) & (r >= validint)
319  ivalid = np.where(valid)
320  xvalid = ivalid[0] + 1
321  yvalid = ivalid[1] + 1
322  ivalid = (xvalid, yvalid)
323  d = f.attrs["_Param_tel__diam"] / (f.attrs["_Param_dm__nact"][0] - 1)
324  r0 = f.attrs["_Param_atmos__r0"] * (f.attrs["_Param_target__Lambda"] / 0.5)**(
325  6. / 5.)
326  RASC = 180 / np.pi * 3600.
327 
328  scale = 0.23 * (d / r0)**(5 / 3.) * \
329  (f.attrs["_Param_target__Lambda"] * 1e-6 / (2 * np.pi * d))**2 * RASC**2
330 
331  mask = np.zeros((nssp + 2, nssp + 2))
332  Ca = np.identity(nsub * 2)
333 
334  for k in range(nsub):
335  mask *= 0
336  mask[xvalid[k], yvalid[k]] = 1
337  mask[xvalid[k], yvalid[k] - 1] = -0.5
338  mask[xvalid[k], yvalid[k] + 1] = -0.5
339  Ca[k, :nsub] = mask[ivalid].flatten()
340  mask *= 0
341  mask[xvalid[k], yvalid[k]] = 1
342  mask[xvalid[k] - 1, yvalid[k]] = -0.5
343  mask[xvalid[k] + 1, yvalid[k]] = -0.5
344  Ca[k + nsub, nsub:] = mask[ivalid].flatten()
345 
346  R = f["R"][:]
347  Ca = R.dot(Ca * scale).dot(R.T)
348  if (modal):
349  P = f["P"][:]
350  Ca = P.dot(Ca).dot(P.T)
351  f.close()
352  return Ca
353 
354 
355 def compute_Cn_cpu(filename, model="data", modal=True):
356  """ Returns the noise error covariance matrix using CPU version of GROOT
357  from a ROKET file
358  :parameter:
359  filename : (string) : full path to the ROKET file
360  modal : (bool) : if True (default), Cn is returned in the Btt modal basis,
361  in the actuator basis if False
362  :return:
363  Cn : (np.ndarray(dim=2, dtype=np.float32)) : noise error covariance matrix
364  """
365  f = h5py.File(filename, 'r')
366  if (model == "data"):
367  N = f["noise"][:]
368  Cn = N.dot(N.T) / N.shape[1]
369  if modal:
370  P = f["P"][:]
371  Cn = P.dot(Cn).dot(P.T)
372  else:
373  nslopes = f["R"][:].shape[1]
374  Cn = np.zeros(nslopes)
375  noise = f.attrs["_Param_wfs__noise"][0]
376  RASC = 180 / np.pi * 3600.
377  if (noise >= 0):
378  Nph = f.attrs["_Param_wfs__zerop"] * 10 ** (-0.4 * f.attrs["_Param_wfs__gsmag"]) * \
379  f.attrs["_Param_wfs__optthroughput"] * \
380  (f.attrs["_Param_tel__diam"] / f.attrs["_Param_wfs__nxsub"]
381  ) ** 2. * f.attrs["_Param_loop__ittime"]
382 
383  r0 = (f.attrs["_Param_wfs__Lambda"] / 0.5)**(
384  6.0 / 5.0) * f.attrs["_Param_atmos__r0"]
385 
386  sig = (np.pi ** 2 / 2) * (1 / Nph) * \
387  (1. / r0) ** 2 # Photon noise in m^-2
388  # Noise variance in arcsec^2
389  sig = sig * (
390  (f.attrs["_Param_wfs__Lambda"] * 1e-6) / (2 * np.pi))**2 * RASC**2
391 
392  Ns = f.attrs["_Param_wfs__npix"] # Number of pixel
393  Nd = (f.attrs["_Param_wfs__Lambda"] *
394  1e-6) * RASC / f.attrs["_Param_wfs__pixsize"]
395  sigphi = (np.pi ** 2 / 3.0) * (1 / Nph ** 2) * (f.attrs["_Param_wfs__noise"]) ** 2 * \
396  Ns ** 2 * (Ns / Nd) ** 2 # Phase variance in m^-2
397  # Noise variance in arcsec^2
398  sigsh = sigphi * \
399  ((f.attrs["_Param_wfs__Lambda"] * 1e-6) / (2 * np.pi)) ** 2 * RASC ** 2
400 
401  Cn[:len(sig)] = sig + sigsh
402  Cn[len(sig):] = sig + sigsh
403 
404  Cn = np.diag(Cn)
405  R = f["R"][:]
406  Cn = R.dot(Cn).dot(R.T)
407  if (modal):
408  P = f["P"][:]
409  Cn = P.dot(Cn).dot(P.T)
410  f.close()
411  return Cn
412 
413 
414 def compute_OTF_fitting(filename, otftel):
415  """
416  Modelize the OTF due to the fitting using dphi_highpass
417 
418  :parameters:
419  filename: (str) : ROKET hdf5 file path
420  otftel: (np.ndarray) : Telescope OTF
421  :return:
422  otf_fit: (np.ndarray) : Fitting OTF
423  psf_fit (np.ndarray) : Fitting PSF
424  """
425  f = h5py.File(filename, 'r')
426  r0 = f.attrs["_Param_atmos__r0"] * (f.attrs["_Param_target__Lambda"][0] / 0.5)**(
427  6. / 5.)
428  ratio_lambda = 2 * np.pi / f.attrs["_Param_target__Lambda"][0]
429  # Telescope OTF
430  spup = drax.get_pup(filename)
431  mradix = 2
432  fft_size = mradix**int((np.log(2 * spup.shape[0]) / np.log(mradix)) + 1)
433  mask = np.ones((fft_size, fft_size))
434  mask[np.where(otftel < 1e-5)] = 0
435 
436  x = np.arange(fft_size) - fft_size / 2
437  pixsize = f.attrs["_Param_tel__diam"] / f.attrs["_Param_geom__pupdiam"]
438  x = x * pixsize
439  r = np.sqrt(x[:, None] * x[:, None] + x[None, :] * x[None, :])
440  tabx, taby = starlord.tabulateIj0()
441  dphi = np.fft.fftshift(
442  starlord.dphi_highpass(
443  r, f.attrs["_Param_tel__diam"] / (f.attrs["_Param_dm__nact"][0] - 1),
444  tabx, taby) * (1 / r0)**(5 / 3.)) # * den * ratio_lambda**2 * mask
445  otf_fit = np.exp(-0.5 * dphi) * mask
446  otf_fit = otf_fit / otf_fit.max()
447 
448  psf_fit = np.fft.fftshift(np.real(np.fft.ifft2(otftel * otf_fit)))
449  psf_fit *= (fft_size * fft_size / float(np.where(spup)[0].shape[0]))
450 
451  f.close()
452  return otf_fit, psf_fit
453 
454 
455 def compute_PSF(filename):
456  """
457  Modelize the PSF using GROOT model for aniso and bandwidth, Gendron model for aliasing,
458  dphi_highpass for fitting, noise extracted from datas. Non linearity not taken into account
459 
460  :parameters:
461  filename: (str) : ROKET hdf5 file path
462  :return:
463  psf: (np.ndarray) : PSF
464  """
465  tic = time.time()
466  spup = drax.get_pup(filename)
467  Cab = compute_Cerr(filename)
468  Cn = compute_Cn_cpu(filename)
469  Ca = compute_Calias(filename)
470  Cee = Cab + Cn + Ca
471  otftel, otf2, psf, gpu = gamora.psf_rec_Vii(filename, fitting=False,
472  cov=(Cee).astype(np.float32))
473  otf_fit, psf_fit = compute_OTF_fitting(filename, otftel)
474  psf = np.fft.fftshift(np.real(np.fft.ifft2(otf_fit * otf2 * otftel)))
475  psf *= (psf.shape[0] * psf.shape[0] / float(np.where(spup)[0].shape[0]))
476  tac = time.time()
477  print("PSF computed in ", tac - tic, " seconds")
478 
479  return psf
480 
481 
482 def compute_Calias_gpu(filename, slopes_space=False, modal=True, npts=3):
483  f = h5py.File(filename, 'r')
484  nsub = f["R"][:].shape[1] // 2
485  nssp = f.attrs["_Param_wfs__nxsub"][0]
486  npix = f.attrs["_Param_wfs__npix"][0]
487  validint = f.attrs["_Param_tel__cobs"]
488  x = np.linspace(-1, 1, nssp)
489  x, y = np.meshgrid(x, x)
490  r = np.sqrt(x * x + y * y)
491  rorder = np.sort(r.reshape(nssp * nssp))
492  ncentral = nssp * nssp - np.sum(r >= validint, dtype=np.int32)
493  validext = rorder[ncentral + nsub]
494  valid = (r < validext) & (r >= validint)
495  ivalid = np.where(valid)
496  r0 = f.attrs["_Param_atmos__r0"]
497  Lambda_wfs = f.attrs["_Param_wfs__Lambda"][0]
498  d = f.attrs["_Param_tel__diam"] / nssp
499  RASC = 180 / np.pi * 3600
500  scale = 0.5 * (1 / r0)**(5 / 3)
501  c = (RASC * Lambda_wfs * 1e-6 / 2 / np.pi) / d**2
502  h = d / (npts - 1)
503  x = (np.arange(nssp) - nssp / 2) * d
504  x, y = np.meshgrid(x, x)
505  x = x[ivalid].astype(np.float32)
506  y = y[ivalid].astype(np.float32)
507  fc = 1 / (2 * d) #/ npix
508  scale = scale * c**2 * (h / 3)**2
509  coeff = simpson_coeff(npts)
510  weights = np.zeros(npts)
511  for k in range(npts):
512  weights[k] = (coeff[k:] * coeff[:npts - k]).sum()
513  groot = Groot(cxt, cxt.active_device, nsub, weights.astype(np.float32), scale, x, y,
514  fc, d, npts)
515  groot.compute_Calias()
516  CaXX = np.array(groot.d_CaXX)
517  Ca = np.zeros((2 * CaXX.shape[0], 2 * CaXX.shape[0]))
518  Ca[:CaXX.shape[0], :CaXX.shape[0]] = CaXX
519  Ca[CaXX.shape[0]:, CaXX.shape[0]:] = np.array(groot.d_CaYY)
520  if not slopes_space:
521  R = f["R"][:]
522  Ca = R.dot(Ca).dot(R.T)
523  if modal:
524  P = f["P"][:]
525  Ca = P.dot(Ca).dot(P.T)
526  f.close()
527 
528  return Ca
529 
530 
531 def compute_Calias(filename, slopes_space=False, modal=True, npts=3):
532  """ Returns the aliasing slopes covariance matrix using CPU version of GROOT
533  from a ROKET file and a model based on structure function
534  :parameter:
535  filename : (string) : full path to the ROKET file
536  slopes_space: (bool): (optionnal) if True, return the covariance matrix in the slopes space
537  modal: (bool): (optionnal) if True, return the covariance matrix in the modal space
538  :return:
539  Ca : (np.ndarray(dim=2, dtype=np.float32)) : aliasing error covariance matrix
540  """
541 
542  f = h5py.File(filename, 'r')
543  tabx, taby = starlord.tabulateIj0()
544  nsub = f["R"][:].shape[1] // 2
545  nssp = f.attrs["_Param_wfs__nxsub"][0]
546  npix = f.attrs["_Param_wfs__npix"][0]
547  validint = f.attrs["_Param_tel__cobs"]
548  x = np.linspace(-1, 1, nssp)
549  x, y = np.meshgrid(x, x)
550  r = np.sqrt(x * x + y * y)
551  rorder = np.sort(r.reshape(nssp * nssp))
552  ncentral = nssp * nssp - np.sum(r >= validint, dtype=np.int32)
553  validext = rorder[ncentral + nsub]
554  valid = (r < validext) & (r >= validint)
555  ivalid = np.where(valid)
556  r0 = f.attrs["_Param_atmos__r0"]
557  Lambda_wfs = f.attrs["_Param_wfs__Lambda"][0]
558  d = f.attrs["_Param_tel__diam"] / nssp
559  RASC = 180 / np.pi * 3600
560  scale = 0.5 * (1 / r0)**(5 / 3)
561  c = (RASC * Lambda_wfs * 1e-6 / 2 / np.pi) / d**2
562  x = (np.arange(nssp) - nssp / 2) * d
563  x, y = np.meshgrid(x, x)
564  x = x[ivalid]
565  y = y[ivalid]
566  fc = d #/ npix
567  xx = np.tile(x, (nsub, 1))
568  yy = np.tile(y, (nsub, 1))
569  # Ca = compute_Calias_element(xx, yy, fc, d, nsub, tabx, taby)
570  # Ca += compute_Calias_element(xx, yy, fc, d, nsub, tabx, taby, xoff=0.5)
571  # Ca += compute_Calias_element(xx, yy, fc, d, nsub, tabx, taby, xoff=-0.5)
572  # Ca += compute_Calias_element(xx, yy, fc, d, nsub, tabx, taby, yoff=0.5)
573  # Ca += compute_Calias_element(xx, yy, fc, d, nsub, tabx, taby, yoff=-0.5)
574  # Ca = Ca * scale / 5
575  Ca = np.zeros((2 * nsub, 2 * nsub))
576  coeff = simpson_coeff(npts)
577  # for k in tqdm(range(npts)):
578  # weight = (coeff[k:] * coeff[:npts - k]).sum()
579  # Ca += compute_Calias_element_XX(xx, yy, fc, d, nsub, tabx, taby, yoff=k /
580  # (npts - 1)) * weight
581  # Ca += compute_Calias_element_YY(xx, yy, fc, d, nsub, tabx, taby, xoff=k /
582  # (npts - 1)) * weight
583  # if k > 0:
584  # Ca += compute_Calias_element_XX(xx, yy, fc, d, nsub, tabx, taby, yoff=-k /
585  # (npts - 1)) * weight
586  # Ca += compute_Calias_element_YY(xx, yy, fc, d, nsub, tabx, taby, xoff=-k /
587  # (npts - 1)) * weight
588  if (npts > 1):
589  h = d / (npts - 1)
590  else:
591  h = 1
592  for k in tqdm(range(npts)):
593  for p in tqdm(range(npts)):
594  Ca += (compute_Calias_element_XX(xx, yy, fc, d, nsub, tabx, taby,
595  yoff=(k - p) * h) * coeff[k] * coeff[p])
596  Ca += (compute_Calias_element_YY(xx, yy, fc, d, nsub, tabx, taby,
597  xoff=(k - p) * h) * coeff[k] * coeff[p])
598 
599  if not slopes_space:
600  R = f["R"][:]
601  Ca = R.dot(Ca).dot(R.T)
602  if modal:
603  P = f["P"][:]
604  Ca = P.dot(Ca).dot(P.T)
605  f.close()
606 
607  return Ca * scale * c**2 * (h / 3)**2
608 
609 
610 def simpson_coeff(n):
611  """
612  Returns the n weights to apply for a Simpson integration on n elements
613  :parameters:
614  n: (int): number of elements, must be odd
615  :return:
616  coeff: (np.array[ndims=1,dtype=np.int64]): simpson coefficients
617  """
618  if (n == 1):
619  coeff = np.ones(n)
620  else:
621  if (n % 2):
622  coeff = np.ones(n)
623  coeff[1::2] = 4
624  coeff[2:-1:2] = 2
625  else:
626  raise ValueError("n must be odd")
627 
628  return coeff
629 
630 
631 def compute_Calias_element_XX(xx, yy, fc, d, nsub, tabx, taby, xoff=0, yoff=0):
632  """
633  Compute the element of the aliasing covariance matrix
634 
635  :parameters:
636  Ca: (np.ndarray(ndim=2, dtype=np.float32)): aliasing covariance matrix to fill
637  xx: (np.ndarray(ndim=2, dtype=np.float32)): X positions of the WFS subap
638  yy: (np.ndarray(ndim=2, dtype=np.float32)): Y positions of the WFS subap
639  fc: (float): cut-off frequency for structure function
640  d: (float): subap diameter
641  nsub: (int): number of subap
642  tabx: (np.ndarray(ndim=1, dtype=np.float32)): X tabulation for dphi
643  taby: (np.ndarray(ndim=1, dtype=np.float32)): Y tabulation for dphi
644  xoff: (float) : (optionnal) offset to apply on the WFS xpos (units of d)
645  yoff: (float) : (optionnal) offset to apply on the WFS ypos (units of d)
646  """
647  xx = xx - xx.T #+ xoff * d
648  yy = yy - yy.T #+ yoff * d
649  #xx = np.triu(xx) - np.triu(xx, -1).T
650  #yy = np.triu(yy) - np.triu(yy, -1).T
651  Ca = np.zeros((2 * nsub, 2 * nsub))
652 
653  # XX covariance
654  AB = np.linalg.norm([xx, yy + yoff], axis=0)
655  Ab = np.linalg.norm([xx - d, yy + yoff], axis=0)
656  aB = np.linalg.norm([xx + d, yy + yoff], axis=0)
657  ab = AB
658 
659  Ca[:nsub, :nsub] += starlord.dphi_highpass(
660  Ab, fc, tabx, taby) + starlord.dphi_highpass(
661  aB, fc, tabx, taby) - 2 * starlord.dphi_highpass(AB, fc, tabx, taby)
662 
663  return Ca
664 
665 
666 def compute_Calias_element_YY(xx, yy, fc, d, nsub, tabx, taby, xoff=0, yoff=0):
667  """
668  Compute the element of the aliasing covariance matrix
669 
670  :parameters:
671  Ca: (np.ndarray(ndim=2, dtype=np.float32)): aliasing covariance matrix to fill
672  xx: (np.ndarray(ndim=2, dtype=np.float32)): X positions of the WFS subap
673  yy: (np.ndarray(ndim=2, dtype=np.float32)): Y positions of the WFS subap
674  fc: (float): cut-off frequency for structure function
675  d: (float): subap diameter
676  nsub: (int): number of subap
677  tabx: (np.ndarray(ndim=1, dtype=np.float32)): X tabulation for dphi
678  taby: (np.ndarray(ndim=1, dtype=np.float32)): Y tabulation for dphi
679  xoff: (float) : (optionnal) offset to apply on the WFS xpos (units of d)
680  yoff: (float) : (optionnal) offset to apply on the WFS ypos (units of d)
681  """
682  xx = xx - xx.T #+ xoff * d
683  yy = yy - yy.T #+ yoff * d
684  #xx = np.triu(xx) - np.triu(xx, -1).T
685  #yy = np.triu(yy) - np.triu(yy, -1).T
686  Ca = np.zeros((2 * nsub, 2 * nsub))
687 
688  # YY covariance
689  CD = np.linalg.norm([xx + xoff, yy], axis=0)
690  Cd = np.linalg.norm([xx + xoff, yy - d], axis=0)
691  cD = np.linalg.norm([xx + xoff, yy + d], axis=0)
692  cd = CD
693 
694  Ca[nsub:, nsub:] += starlord.dphi_highpass(
695  Cd, fc, tabx, taby) + starlord.dphi_highpass(
696  cD, fc, tabx, taby) - 2 * starlord.dphi_highpass(CD, fc, tabx, taby)
697 
698  return Ca
699 
700 
701 def compute_Calias_element_XY(xx, yy, fc, d, nsub, tabx, taby, xoff=0, yoff=0):
702  """
703  Compute the element of the aliasing covariance matrix
704 
705  :parameters:
706  Ca: (np.ndarray(ndim=2, dtype=np.float32)): aliasing covariance matrix to fill
707  xx: (np.ndarray(ndim=2, dtype=np.float32)): X positions of the WFS subap
708  yy: (np.ndarray(ndim=2, dtype=np.float32)): Y positions of the WFS subap
709  fc: (float): cut-off frequency for struture function
710  d: (float): subap diameter
711  nsub: (int): number of subap
712  tabx: (np.ndarray(ndim=1, dtype=np.float32)): X tabulation for dphi
713  taby: (np.ndarray(ndim=1, dtype=np.float32)): Y tabulation for dphi
714  xoff: (float) : (optionnal) offset to apply on the WFS xpos (units of d)
715  yoff: (float) : (optionnal) offset to apply on the WFS ypos (units of d)
716  """
717  xx = xx - xx.T + xoff * d
718  yy = yy - yy.T + yoff * d
719  Ca = np.zeros((2 * nsub, 2 * nsub))
720 
721  # YY covariance
722  aD = np.linalg.norm([xx + d / 2, yy + d / 2], axis=0)
723  ad = np.linalg.norm([xx + d / 2, yy - d / 2], axis=0)
724  Ad = np.linalg.norm([xx - d / 2, yy - d / 2], axis=0)
725  AD = np.linalg.norm([xx - d / 2, yy + d / 2], axis=0)
726 
727  Ca[nsub:, :nsub] = 0.25 * (
728  starlord.dphi_highpass(Ad, d, tabx, taby) + starlord.dphi_highpass(
729  aD, d, tabx, taby) - starlord.dphi_highpass(AD, d, tabx, taby) -
730  starlord.dphi_highpass(ad, d, tabx, taby))
731  Ca[:nsub, nsub:] = Ca[nsub:, :nsub].copy()
732  return Ca
733 
734 
735 def compute_Calias_element(xx, yy, fc, d, nsub, tabx, taby, xoff=0, yoff=0):
736  """
737  Compute the element of the aliasing covariance matrix
738 
739  :parameters:
740  Ca: (np.ndarray(ndim=2, dtype=np.float32)): aliasing covariance matrix to fill
741  xx: (np.ndarray(ndim=2, dtype=np.float32)): X positions of the WFS subap
742  yy: (np.ndarray(ndim=2, dtype=np.float32)): Y positions of the WFS subap
743  fc: (float): cut-off frequency for structure function
744  d: (float): subap diameter
745  nsub: (int): number of subap
746  tabx: (np.ndarray(ndim=1, dtype=np.float32)): X tabulation for dphi
747  taby: (np.ndarray(ndim=1, dtype=np.float32)): Y tabulation for dphi
748  xoff: (float) : (optionnal) offset to apply on the WFS xpos (units of d)
749  yoff: (float) : (optionnal) offset to apply on the WFS ypos (units of d)
750  """
751  xx = xx - xx.T + xoff * d
752  yy = yy - yy.T + yoff * d
753  Ca = np.zeros((2 * nsub, 2 * nsub))
754 
755  # XX covariance
756  AB = np.linalg.norm([xx, yy], axis=0)
757  Ab = np.linalg.norm([xx - d, yy], axis=0)
758  aB = np.linalg.norm([xx + d, yy], axis=0)
759  ab = AB
760 
761  Ca[:nsub, :nsub] += starlord.dphi_highpass(
762  Ab, fc, tabx, taby) + starlord.dphi_highpass(
763  aB, fc, tabx, taby) - 2 * starlord.dphi_highpass(AB, fc, tabx, taby)
764 
765  # YY covariance
766  CD = AB
767  Cd = np.linalg.norm([xx, yy - d], axis=0)
768  cD = np.linalg.norm([xx, yy + d], axis=0)
769  cd = CD
770 
771  Ca[nsub:, nsub:] += starlord.dphi_highpass(
772  Cd, fc, tabx, taby) + starlord.dphi_highpass(
773  cD, fc, tabx, taby) - 2 * starlord.dphi_highpass(CD, fc, tabx, taby)
774 
775  # XY covariance
776 
777  # aD = np.linalg.norm([xx + d/2, yy + d/2], axis=0)
778  # ad = np.linalg.norm([xx + d/2, yy - d/2], axis=0)
779  # Ad = np.linalg.norm([xx - d/2, yy - d/2], axis=0)
780  # AD = np.linalg.norm([xx - d/2, yy + d/2], axis=0)
781  #
782  # Ca[nsub:,:nsub] = 0.25 * (starlord.dphi_highpass(Ad, d, tabx, taby)
783  # + starlord.dphi_highpass(aD, d, tabx, taby)
784  # - starlord.dphi_highpass(AD, d, tabx, taby)
785  # - starlord.dphi_highpass(ad, d, tabx, taby)) * (1 / r0)**(5. / 3.)
786  # Ca[:nsub,nsub:] = Ca[nsub:,:nsub].copy()
787  return Ca
788 
789 
790 def compute_dCmm(filename, ws=None, wd=None, dk=1):
791  """ Returns the derivative slopes covariance matrix using CPU version of GROOT
792  from a ROKET file and a model based on structure function
793  :parameter:
794  filename : (string) : full path to the ROKET file
795  ws: (np.array[ndim=1, dtype=np.float32]): wind speed per layer [m/s]
796  wd: (np.array[ndim=1, dtype=np.float32]): wind direction per layer [deg]
797  dk: (int): slopes shift [iterations]
798  :return:
799  dCmm : (np.ndarray(dim=2, dtype=np.float32)) : d/dt(slopes)*slopes
800  """
801 
802  f = h5py.File(filename, 'r')
803  if ws is None:
804  ws = f.attrs["_Param_atmos__windspeed"]
805  if wd is None:
806  wd = f.attrs["_Param_atmos__winddir"]
807  dt = f.attrs["_Param_loop__ittime"] * dk
808  L0 = f.attrs["_Param_atmos__L0"]
809  frac = f.attrs["_Param_atmos__frac"]
810  nsub = f["R"][:].shape[1] // 2
811  nssp = f.attrs["_Param_wfs__nxsub"][0]
812  validint = f.attrs["_Param_tel__cobs"]
813  x = np.linspace(-1, 1, nssp)
814  x, y = np.meshgrid(x, x)
815  r = np.sqrt(x * x + y * y)
816  rorder = np.sort(r.reshape(nssp * nssp))
817  ncentral = nssp * nssp - np.sum(r >= validint, dtype=np.int32)
818  validext = rorder[ncentral + nsub]
819  valid = (r < validext) & (r >= validint)
820  ivalid = np.where(valid)
821  r0 = f.attrs["_Param_atmos__r0"]
822  Lambda_wfs = f.attrs["_Param_wfs__Lambda"][0]
823  d = f.attrs["_Param_tel__diam"] / nssp
824  RASC = 180 / np.pi * 3600
825  scale = 0.5 * (1 / r0)**(5 / 3) * (RASC * Lambda_wfs * 1e-6 / 2 / np.pi)**2 / d**2
826  x = (np.arange(nssp) - nssp / 2) * d
827  x, y = np.meshgrid(x, x)
828  x = x[ivalid]
829  y = y[ivalid]
830  xx = np.tile(x, (nsub, 1))
831  yy = np.tile(y, (nsub, 1))
832  f.close()
833  dCmm = np.zeros((2 * nsub, 2 * nsub))
834  for k in range(ws.size):
835  dCmm += frac[k] * compute_dCmm_element(xx, yy, d, nsub, ws[k], wd[k], dt, L0[k])
836 
837  return dCmm * scale
838 
839 
840 def compute_dCmm_element(xx, yy, d, nsub, ws, wd, dt, L0):
841  """
842  Compute the element of the derivative slopes covariance matrix
843 
844  :parameters:
845  xx: (np.ndarray(ndim=2, dtype=np.float32)): X positions of the WFS subap
846  yy: (np.ndarray(ndim=2, dtype=np.float32)): Y positions of the WFS subap
847  d: (float): subap diameter
848  nsub: (int): number of subap
849  ws: (float): wind speed per layer [m/s]
850  wd: (float): wind direction per layer [deg]
851  dt: (float): iteration time [s]
852  L0: (float): outer scale [m]
853  """
854  xij = xx - xx.T
855  yij = yy - yy.T
856  dCmm = np.zeros((2 * nsub, 2 * nsub))
857  vdt = ws * dt
858  wd = wd / 180 * np.pi
859 
860  # XX covariance
861  AB = np.linalg.norm([-xij + vdt * np.cos(wd), -yij + vdt * np.sin(wd)], axis=0)
862  Ab = np.linalg.norm([-xij - d + vdt * np.cos(wd), -yij + vdt * np.sin(wd)], axis=0)
863  aB = np.linalg.norm([-xij + d + vdt * np.cos(wd), -yij + vdt * np.sin(wd)], axis=0)
864 
865  dCmm[:nsub, :nsub] += starlord.rodconan(Ab, L0) + starlord.rodconan(
866  aB, L0) - 2 * starlord.rodconan(AB, L0)
867 
868  AB = np.linalg.norm([xij + vdt * np.cos(wd), yij + vdt * np.sin(wd)], axis=0)
869  Ab = np.linalg.norm([xij - d + vdt * np.cos(wd), yij + vdt * np.sin(wd)], axis=0)
870  aB = np.linalg.norm([xij + d + vdt * np.cos(wd), yij + vdt * np.sin(wd)], axis=0)
871 
872  dCmm[:nsub, :nsub] -= (starlord.rodconan(Ab, L0) + starlord.rodconan(aB, L0) -
873  2 * starlord.rodconan(AB, L0))
874 
875  # YY covariance
876  CD = np.linalg.norm([-xij + vdt * np.cos(wd), -yij + vdt * np.sin(wd)], axis=0)
877  Cd = np.linalg.norm([-xij + vdt * np.cos(wd), -yij - d + vdt * np.sin(wd)], axis=0)
878  cD = np.linalg.norm([-xij + vdt * np.cos(wd), -yij + d + vdt * np.sin(wd)], axis=0)
879 
880  dCmm[nsub:, nsub:] += starlord.rodconan(Cd, L0) + starlord.rodconan(
881  cD, L0) - 2 * starlord.rodconan(CD, L0)
882 
883  CD = np.linalg.norm([xij + vdt * np.cos(wd), yij + vdt * np.sin(wd)], axis=0)
884  Cd = np.linalg.norm([xij + vdt * np.cos(wd), yij - d + vdt * np.sin(wd)], axis=0)
885  cD = np.linalg.norm([xij + vdt * np.cos(wd), yij + d + vdt * np.sin(wd)], axis=0)
886 
887  dCmm[nsub:, nsub:] -= (starlord.rodconan(Cd, L0) + starlord.rodconan(cD, L0) -
888  2 * starlord.rodconan(CD, L0))
889  # XY covariance
890 
891  # aD = np.linalg.norm([xx + d/2, yy + d/2], axis=0)
892  # ad = np.linalg.norm([xx + d/2, yy - d/2], axis=0)
893  # Ad = np.linalg.norm([xx - d/2, yy - d/2], axis=0)
894  # AD = np.linalg.norm([xx - d/2, yy + d/2], axis=0)
895  #
896  # dCmm[nsub:,:nsub] = 0.25 * (starlord.rodconan(Ad, d, tabx, taby)
897  # + starlord.dphi_highpass(aD, d, tabx, taby)
898  # - starlord.dphi_highpass(AD, d, tabx, taby)
899  # - starlord.dphi_highpass(ad, d, tabx, taby)) * (1 / r0)**(5. / 3.)
900  # dCmm[:nsub,nsub:] = dCmm[nsub:,:nsub].copy()
901  return 0.25 * dCmm
guardians.groot.simpson_coeff
def simpson_coeff(n)
Definition: groot.py:617
guardians.groot.compute_Calias
def compute_Calias(filename, slopes_space=False, modal=True, npts=3)
Returns the aliasing slopes covariance matrix using CPU version of GROOT from a ROKET file and a mode...
Definition: groot.py:540
guardians.groot.compute_PSF
def compute_PSF(filename)
Definition: groot.py:464
shesha.sutra_wrap
Definition: sutra_wrap.py:1
guardians.groot.test_Cerr
def test_Cerr(filename)
Compute PSF of aniso and bandwidth from GROOT model and ROKET to compare.
Definition: groot.py:218
guardians.groot.compute_Ca_cpu
def compute_Ca_cpu(filename, modal=True)
Returns the aliasing error covariance matrix using CPU version of GROOT from a ROKET file :parameter:...
Definition: groot.py:307
guardians.groot.compute_dCmm_element
def compute_dCmm_element(xx, yy, d, nsub, ws, wd, dt, L0)
Compute the element of the derivative slopes covariance matrix.
Definition: groot.py:853
guardians.groot.compute_Cn_cpu
def compute_Cn_cpu(filename, model="data", modal=True)
Returns the noise error covariance matrix using CPU version of GROOT from a ROKET file :parameter: fi...
Definition: groot.py:364
guardians.groot.compute_Calias_element_XY
def compute_Calias_element_XY(xx, yy, fc, d, nsub, tabx, taby, xoff=0, yoff=0)
Compute the element of the aliasing covariance matrix.
Definition: groot.py:716
guardians.groot.compute_Calias_element_XX
def compute_Calias_element_XX(xx, yy, fc, d, nsub, tabx, taby, xoff=0, yoff=0)
Compute the element of the aliasing covariance matrix.
Definition: groot.py:646
guardians.groot.compute_Calias_element_YY
def compute_Calias_element_YY(xx, yy, fc, d, nsub, tabx, taby, xoff=0, yoff=0)
Compute the element of the aliasing covariance matrix.
Definition: groot.py:681
guardians.groot.compare_GPU_vs_CPU
def compare_GPU_vs_CPU(filename)
Compare results of GROOT vs its CPU version in terms of execution time and precision on the PSF renco...
Definition: groot.py:253
guardians.groot.compute_Calias_element
def compute_Calias_element(xx, yy, fc, d, nsub, tabx, taby, xoff=0, yoff=0)
Compute the element of the aliasing covariance matrix.
Definition: groot.py:750
guardians.groot.compute_Cerr_cpu
def compute_Cerr_cpu(filename, modal=True)
Returns the residual error covariance matrix using CPU version of GROOT from a ROKET file :parameter:...
Definition: groot.py:117
guardians.groot.compute_Calias_gpu
def compute_Calias_gpu(filename, slopes_space=False, modal=True, npts=3)
Definition: groot.py:482
guardians.groot.compute_Cerr
def compute_Cerr(filename, modal=True, ctype="float", speed=None, H=None, theta=None, r0=None, L0=None)
Returns the residual error covariance matrix using GROOT from a ROKET file :parameter: filename (stri...
Definition: groot.py:37
guardians.groot.compute_dCmm
def compute_dCmm(filename, ws=None, wd=None, dk=1)
Returns the derivative slopes covariance matrix using CPU version of GROOT from a ROKET file and a mo...
Definition: groot.py:800
guardians.groot.compute_OTF_fitting
def compute_OTF_fitting(filename, otftel)
Modelize the OTF due to the fitting using dphi_highpass.
Definition: groot.py:424