2 GROOT (Gpu-based Residual errOr cOvariance maTrix)
3 Python module for modelization of error covariance matrix
14 from guardians
import drax, starlord
15 import matplotlib.pyplot
as plt
18 gpudevices = np.array([0, 1, 2, 3], dtype=np.int32)
19 cxt = carmaWrap_context.get_instance_ngpu(gpudevices.size, gpudevices)
22 def compute_Cerr(filename, modal=True, ctype="float", speed=None, H=None, theta=None,
24 """ Returns the residual error covariance matrix using GROOT from a ROKET file
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]
37 Cerr : (np.ndarray(dim=2, dtype=np.float32)) : residual error covariance matrix
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]
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
57 H = f.attrs[
"_Param_atmos__alt"]
59 L0 = f.attrs[
"_Param_atmos__L0"]
61 speed = f.attrs[
"_Param_atmos__windspeed"]
63 theta = (f.attrs[
"_Param_atmos__winddir"] * np.pi / 180.)
64 frac = f.attrs[
"_Param_atmos__frac"]
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
72 Nact = np.linalg.inv(Nact)
75 Tf = Btt[:-2, :-2].dot(P[:-2, :-2])
76 IF, T = drax.get_IF(filename)
80 deltaTT = T.T.dot(T) / N
81 deltaF = IF.T.dot(T) / N
82 pzt2tt = np.linalg.inv(deltaTT).dot(deltaF.T)
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)
92 raise TypeError(
"Unknown ctype : must be float")
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)
100 print(
"Cee computed in : %.2f seconds" % (tac - tic))
102 cov_err_groot = P.dot(cov_err_groot).dot(P.T)
109 """ Returns the residual error covariance matrix using CPU version of GROOT
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
116 Cerr : (np.ndarray(dim=2, dtype=np.float32)) : residual error covariance matrix
118 f = h5py.File(filename,
'r')
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))
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]
150 Htheta = np.linalg.norm([wxpos, wypos]) / RASC * H
151 vdt = speed * dt / gain
153 M = np.zeros((xpos.size, xpos.size))
157 angleht = np.arctan2(wypos, wxpos)
158 fc = xactu[1] - xactu[0]
160 M = np.linalg.norm([xij, yij], axis=0)
161 Mvdt = np.linalg.norm([xij - vdt * np.cos(theta), yij - vdt * np.sin(theta)],
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)
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
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
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)
188 Tf = Btt[:-2, :-2].dot(P[:-2, :-2])
190 IF, T = drax.get_IF(filename)
194 deltaTT = T.T.dot(T) / N
195 deltaF = IF.T.dot(T) / N
196 pzt2tt = np.linalg.inv(deltaTT).dot(deltaF.T)
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
207 cov_err = P.dot(cov_err).dot(P.T)
214 """ Compute PSF of aniso and bandwidth from GROOT model and ROKET to compare
217 filename:(str):path to the ROKET file
219 C = drax.get_covmat_contrib(filename, [
"bandwidth",
"tomography"])
221 _, _, psfr, _ = gamora.psf_rec_Vii(filename, covmodes=C.astype(np.float32),
223 _, _, psf, _ = gamora.psf_rec_Vii(filename, cov=Cerr.astype(np.float32),
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)))
234 np.log10(np.abs(psf)), vmax=np.log10(np.abs(psfr)).max(),
235 vmin=np.log10(np.abs(psfr)).min())
239 np.log10(np.abs(psfr - psf)), vmax=np.log10(np.abs(psfr)).max(),
240 vmin=np.log10(np.abs(psfr)).min())
242 plt.title(
"PSF_R - PSF")
248 """ Compare results of GROOT vs its CPU version in terms of execution time
249 and precision on the PSF renconstruction
251 filename : (string) : full path to the ROKET file
254 timer = ch.carmaWrap_timer()
259 synctime = timer.stop()
265 gpu_time_s = timer.stop() - synctime
271 gpu_time_d = timer.stop() - synctime
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))
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)
299 """ Returns the aliasing error covariance matrix using CPU version of GROOT
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
306 Ca : (np.ndarray(dim=2, dtype=np.float32)) : aliasing error covariance matrix
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)**(
326 RASC = 180 / np.pi * 3600.
328 scale = 0.23 * (d / r0)**(5 / 3.) * \
329 (f.attrs[
"_Param_target__Lambda"] * 1e-6 / (2 * np.pi * d))**2 * RASC**2
331 mask = np.zeros((nssp + 2, nssp + 2))
332 Ca = np.identity(nsub * 2)
334 for k
in range(nsub):
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()
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()
347 Ca = R.dot(Ca * scale).dot(R.T)
350 Ca = P.dot(Ca).dot(P.T)
356 """ Returns the noise error covariance matrix using CPU version of GROOT
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
363 Cn : (np.ndarray(dim=2, dtype=np.float32)) : noise error covariance matrix
365 f = h5py.File(filename,
'r')
366 if (model ==
"data"):
368 Cn = N.dot(N.T) / N.shape[1]
371 Cn = P.dot(Cn).dot(P.T)
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.
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"]
383 r0 = (f.attrs[
"_Param_wfs__Lambda"] / 0.5)**(
384 6.0 / 5.0) * f.attrs[
"_Param_atmos__r0"]
386 sig = (np.pi ** 2 / 2) * (1 / Nph) * \
390 (f.attrs[
"_Param_wfs__Lambda"] * 1e-6) / (2 * np.pi))**2 * RASC**2
392 Ns = f.attrs[
"_Param_wfs__npix"]
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
399 ((f.attrs[
"_Param_wfs__Lambda"] * 1e-6) / (2 * np.pi)) ** 2 * RASC ** 2
401 Cn[:len(sig)] = sig + sigsh
402 Cn[len(sig):] = sig + sigsh
406 Cn = R.dot(Cn).dot(R.T)
409 Cn = P.dot(Cn).dot(P.T)
416 Modelize the OTF due to the fitting using dphi_highpass
419 filename: (str) : ROKET hdf5 file path
420 otftel: (np.ndarray) : Telescope OTF
422 otf_fit: (np.ndarray) : Fitting OTF
423 psf_fit (np.ndarray) : Fitting PSF
425 f = h5py.File(filename,
'r')
426 r0 = f.attrs[
"_Param_atmos__r0"] * (f.attrs[
"_Param_target__Lambda"][0] / 0.5)**(
428 ratio_lambda = 2 * np.pi / f.attrs[
"_Param_target__Lambda"][0]
430 spup = drax.get_pup(filename)
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
436 x = np.arange(fft_size) - fft_size / 2
437 pixsize = f.attrs[
"_Param_tel__diam"] / f.attrs[
"_Param_geom__pupdiam"]
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.))
445 otf_fit = np.exp(-0.5 * dphi) * mask
446 otf_fit = otf_fit / otf_fit.max()
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]))
452 return otf_fit, psf_fit
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
461 filename: (str) : ROKET hdf5 file path
463 psf: (np.ndarray) : PSF
466 spup = drax.get_pup(filename)
471 otftel, otf2, psf, gpu = gamora.psf_rec_Vii(filename, fitting=
False,
472 cov=(Cee).astype(np.float32))
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]))
477 print(
"PSF computed in ", tac - tic,
" seconds")
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
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)
508 scale = scale * c**2 * (h / 3)**2
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,
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)
522 Ca = R.dot(Ca).dot(R.T)
525 Ca = P.dot(Ca).dot(P.T)
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
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
539 Ca : (np.ndarray(dim=2, dtype=np.float32)) : aliasing error covariance matrix
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)
567 xx = np.tile(x, (nsub, 1))
568 yy = np.tile(y, (nsub, 1))
575 Ca = np.zeros((2 * nsub, 2 * nsub))
592 for k
in tqdm(range(npts)):
593 for p
in tqdm(range(npts)):
595 yoff=(k - p) * h) * coeff[k] * coeff[p])
597 xoff=(k - p) * h) * coeff[k] * coeff[p])
601 Ca = R.dot(Ca).dot(R.T)
604 Ca = P.dot(Ca).dot(P.T)
607 return Ca * scale * c**2 * (h / 3)**2
612 Returns the n weights to apply for a Simpson integration on n elements
614 n: (int): number of elements, must be odd
616 coeff: (np.array[ndims=1,dtype=np.int64]): simpson coefficients
626 raise ValueError(
"n must be odd")
633 Compute the element of the aliasing covariance matrix
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)
651 Ca = np.zeros((2 * nsub, 2 * nsub))
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)
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)
668 Compute the element of the aliasing covariance matrix
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)
686 Ca = np.zeros((2 * nsub, 2 * nsub))
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)
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)
703 Compute the element of the aliasing covariance matrix
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)
717 xx = xx - xx.T + xoff * d
718 yy = yy - yy.T + yoff * d
719 Ca = np.zeros((2 * nsub, 2 * nsub))
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)
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()
737 Compute the element of the aliasing covariance matrix
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)
751 xx = xx - xx.T + xoff * d
752 yy = yy - yy.T + yoff * d
753 Ca = np.zeros((2 * nsub, 2 * nsub))
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)
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)
767 Cd = np.linalg.norm([xx, yy - d], axis=0)
768 cD = np.linalg.norm([xx, yy + d], axis=0)
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)
791 """ Returns the derivative slopes covariance matrix using CPU version of GROOT
792 from a ROKET file and a model based on structure function
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]
799 dCmm : (np.ndarray(dim=2, dtype=np.float32)) : d/dt(slopes)*slopes
802 f = h5py.File(filename,
'r')
804 ws = f.attrs[
"_Param_atmos__windspeed"]
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)
830 xx = np.tile(x, (nsub, 1))
831 yy = np.tile(y, (nsub, 1))
833 dCmm = np.zeros((2 * nsub, 2 * nsub))
834 for k
in range(ws.size):
842 Compute the element of the derivative slopes covariance matrix
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]
856 dCmm = np.zeros((2 * nsub, 2 * nsub))
858 wd = wd / 180 * np.pi
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)
865 dCmm[:nsub, :nsub] += starlord.rodconan(Ab, L0) + starlord.rodconan(
866 aB, L0) - 2 * starlord.rodconan(AB, L0)
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)
872 dCmm[:nsub, :nsub] -= (starlord.rodconan(Ab, L0) + starlord.rodconan(aB, L0) -
873 2 * starlord.rodconan(AB, L0))
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)
880 dCmm[nsub:, nsub:] += starlord.rodconan(Cd, L0) + starlord.rodconan(
881 cD, L0) - 2 * starlord.rodconan(CD, L0)
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)
887 dCmm[nsub:, nsub:] -= (starlord.rodconan(Cd, L0) + starlord.rodconan(cD, L0) -
888 2 * starlord.rodconan(CD, L0))