2 DRAX (Dedicated functions for Roket file Analysis and eXploitation)
4 Useful functions for ROKET file exploitation
10 import matplotlib.pyplot
as plt
12 from scipy.sparse
import csr_matrix
15 def variance(f, contributors, method="Default"):
16 """ Return the error variance of specified contributors
18 f : (h5py.File) : roket hdf5 file opened with h5py
19 contributors : (list of string) : list of the contributors
20 method : (optional, default="Default") : if "Independence", the
21 function returns ths sum of the contributors variances.
22 If "Default", it returns the variance of the contributors sum
26 swap = np.arange(nmodes) - 2
27 swap[0:2] = [nmodes - 2, nmodes - 1]
28 if (method ==
"Default"):
29 err = f[contributors[0]][:] * 0.
30 for c
in contributors:
32 return np.var(P.dot(err), axis=1)
34 elif (method ==
"Independence"):
37 for c
in contributors:
38 v += np.var(P.dot(f[c][:]), axis=1)
42 raise TypeError(
"Wrong method input")
46 """ Return the variance computed from the sum of contributors of roket
47 files fs, ponderated by frac
49 fs : (list) : list of hdf5 files opened with h5py
50 frac_per_layer : (dict) : frac for each layer
51 contributors : (list of string) : list of the contributors
53 v : (np.array(dim=1)) : variance vector
58 swap = np.arange(nmodes) - 2
59 swap[0:2] = [nmodes - 2, nmodes - 1]
60 err = f[contributors[0]][:] * 0.
62 frac = frac_per_layer[f.attrs[
"_Param_atmos__.alt"][0]]
63 for c
in contributors:
64 err += np.sqrt(frac) * f[c][:]
66 return np.var(P.dot(err), axis=1)
70 """ Returns the cumulative Strehl ratio over the modes from the variance
73 v : (np.array(dim=1)) : variance vector
75 s : (np.array(dim=1)) : cumulative SR
78 s = np.exp(-s * (2 * np.pi / Lambda_tar)**2)
85 Compute the SR over the modes from the variance
89 filename: (str): path to the ROKET file
91 f = h5py.File(filename,
'r')
93 "noise",
"aliasing",
"tomography",
"filtered modes",
"non linearity",
96 if (list(f.attrs.keys()).count(
"_Param_target__Lambda")):
97 Lambda = f.attrs[
"_Param_target__Lambda"][0]
100 nactus = f[
"noise"][:].shape[0]
101 niter = f[
"noise"][:].shape[1]
104 swap = np.arange(nmodes) - 2
105 swap[0:2] = [nmodes - 2, nmodes - 1]
106 data = np.zeros((nmodes, niter))
107 data2 = np.zeros(nmodes)
110 data += np.dot(P, f[i][:])
111 data2 += np.var(np.dot(P, f[i][:]), axis=1)
113 data = np.var(data, axis=1)
114 data = np.cumsum(data[swap])
115 data = np.exp(-data * (2 * np.pi / Lambda)**2)
116 data2 = np.cumsum(data2[swap])
117 data2 = np.exp(-data2 * (2 * np.pi / Lambda)**2)
118 data *= np.exp(-f[
"fitting"].value)
119 data2 *= np.exp(-f[
"fitting"].value)
121 SR2 = np.ones(nmodes) * f[
"SR2"].value
122 SR = np.ones(nmodes) * f[
"SR"].value
124 return data, data2, SR, SR2
129 Return the Modes to Volt matrix
131 filename: (str): path to the ROKET file
133 f = h5py.File(filename,
'r')
139 Return the Volt to Modes matrix
141 filename: (str): path to the ROKET file
143 f = h5py.File(filename,
'r')
149 Return the variance of an error contributor
152 filename: (str): path to the ROKET file
153 contributor: (str): contributor name
155 v: (np.array[ndim=1, dtype=np.float32]): variance of the contributor
157 f = h5py.File(filename,
'r')
160 swap = np.arange(nmodes) - 2
161 swap[0:2] = [nmodes - 2, nmodes - 1]
163 return np.var(np.dot(P, f[contributor][:]), axis=1)
168 Return the sum of the specified contributors error buffers
171 filename: (str): path to the ROKET file
172 contributors: (list): list of contributors
174 err: (np.ndarray[ndim=2,dtype=np.float32]): Sum of the error buffers
176 f = h5py.File(filename,
'r')
178 err = f[
"noise"][:] * 0.
179 for c
in contributors:
188 Return the sum of all the error buffers
191 filename: (str): path to the ROKET file
193 err: (np.ndarray[ndim=2,dtype=np.float32]): Sum of the error buffers
196 f = h5py.File(filename,
'r')
199 err += f[
"aliasing"][:]
200 err += f[
"tomography"][:]
201 err += f[
"filtered modes"][:]
202 err += f[
"non linearity"][:]
203 err += f[
"bandwidth"][:]
211 Return the error covariance matrix considering statistical independence between contributors
214 filename: (str): path to the ROKET file
216 err: (np.ndarray[ndim=2,dtype=np.float32]): Covariance matrix
219 f = h5py.File(filename,
'r')
221 N = f[
"noise"][:].shape[1]
222 err = f[
"noise"][:].dot(f[
"noise"][:].T)
223 err += f[
"aliasing"][:].dot(f[
"aliasing"][:].T)
224 err += f[
"tomography"][:].dot(f[
"tomography"][:].T)
225 err += f[
"filtered modes"][:].dot(f[
"filtered modes"][:].T)
226 err += f[
"non linearity"][:].dot(f[
"non linearity"][:].T)
227 err += f[
"bandwidth"][:].dot(f[
"bandwidth"][:].T)
235 Return the error covariance matrix considering statistical independence between specified contributors
238 filename: (str): path to the ROKET file
239 contributors: (list): list of contributors
241 err: (np.ndarray[ndim=2,dtype=np.float32]): Covariance matrix
244 f = h5py.File(filename,
'r')
246 N = f[
"noise"][:].shape[1]
247 err = np.zeros((f[
"noise"][:].shape[0], f[
"noise"][:].shape[0]))
248 for c
in contributors:
249 err += f[c][:].dot(f[c][:].T)
258 Return the covariance matrix of the specified contributors
261 filename: (str): path to the ROKET file
262 contributor: (list): name of a contributor of the ROKET file
263 modal: (bool): if True (default), return the matrix expressed in the modal basis
265 covmat: (np.ndarray(ndim=2, dtype=np.float32)): covariance matrix
267 h5f = h5py.File(filename,
'r')
268 contrib = h5f[
"bandwidth"][:] * 0.
269 for c
in contributors:
271 covmat = contrib.dot(contrib.T) / contrib.shape[1]
274 covmat = P.dot(covmat).dot(P.T)
282 Return the pupil saved in a ROKET file
284 filename: (str): path to the ROKET file
286 spup: (np.ndarray[ndim=2,dtype=np.float32]): pupil
289 f = h5py.File(filename,
'r')
290 if (list(f.keys()).count(
"spup")):
293 indx_pup = f[
"indx_pup"][:]
294 pup = np.zeros((f[
"dm_dim"].value, f[
"dm_dim"].value))
295 pup_F = pup.flatten()
297 pup = pup_F.reshape(pup.shape)
298 spup = pup[np.where(pup)[0].min():np.where(pup)[0].max() + 1,
299 np.where(pup)[1].min():np.where(pup)[1].max() + 1]
307 Computes the error breakdown in nm rms from a ROKET file
310 filename: (str): path to the ROKET file
312 breakdown: (dict): dictionnary containing the error breakdown
314 f = h5py.File(filename,
'r')
316 noise = f[
"noise"][:]
317 trunc = f[
"non linearity"][:]
318 bp = f[
"bandwidth"][:]
319 tomo = f[
"tomography"][:]
320 aliasing = f[
"aliasing"][:]
321 filt = f[
"filtered modes"][:]
323 swap = np.arange(nmodes) - 2
324 swap[0:2] = [nmodes - 2, nmodes - 1]
325 N = np.var(P.dot(noise), axis=1)
326 S = np.var(P.dot(trunc), axis=1)
327 B = np.var(P.dot(bp), axis=1)
328 T = np.var(P.dot(tomo), axis=1)
329 A = np.var(P.dot(aliasing), axis=1)
330 F = np.var(P.dot(filt), axis=1)
331 C = np.var(P.dot(filt + noise + trunc + bp + tomo + aliasing), axis=1)
332 inde = N + S + B + T + A + F
334 if (list(f.attrs.keys()).count(
"_Param_target__Lambda")):
335 Lambda = f.attrs[
"_Param_target__Lambda"][0]
339 print(
"noise :", np.sqrt(np.sum(N)) * 1e3,
" nm rms")
340 print(
"trunc :", np.sqrt(np.sum(S)) * 1e3,
" nm rms")
341 print(
"bp :", np.sqrt(np.sum(B)) * 1e3,
" nm rms")
342 print(
"tomo :", np.sqrt(np.sum(T)) * 1e3,
" nm rms")
343 print(
"aliasing :", np.sqrt(np.sum(A)) * 1e3,
" nm rms")
344 print(
"filt :", np.sqrt(np.sum(F)) * 1e3,
" nm rms")
346 np.mean(np.sqrt(f[
"fitting"].value / ((2 * np.pi / Lambda)**2)) * 1e3),
348 print(
"cross-terms :", np.sqrt(np.abs(np.sum(C) - np.sum(inde))) * 1e3,
" nm rms")
351 np.sqrt(np.sum(N)) * 1e3,
353 np.sqrt(np.sum(S)) * 1e3,
355 np.sqrt(np.sum(B)) * 1e3,
357 np.sqrt(np.sum(T)) * 1e3,
359 np.sqrt(np.sum(A)) * 1e3,
361 np.sqrt(np.sum(F)) * 1e3,
364 np.sqrt(f[
"fitting"].value /
365 ((2 * np.pi / Lambda)**2)) * 1e3)
446 Displays the covariance and correlation matrix between the contributors
448 filename: (str): path to the ROKET file
449 maparico: (str): (optional) matplotlib colormap to use
452 f = h5py.File(filename,
'r')
456 labels = [
"noise",
"WF deviation",
"aliasing",
"filt. modes",
"bandwidth",
"aniso"]
457 if (maparico
is None):
460 plt.matshow(cov, cmap=maparico)
462 plt.xticks(x, labels, rotation=
"vertical")
463 plt.yticks(x, labels)
465 plt.matshow(cor, cmap=maparico)
467 plt.xticks(x, labels, rotation=
"vertical")
468 plt.yticks(x, labels)
469 print(
"Total variance : ", cov.sum(),
" microns^2")
474 Return the influence functions of the pzt and tt DM saved in a ROKET file
476 filename: (str): path to the ROKET file
478 IF: (csr_matrix): pzt influence function (sparse)
479 T: (np.ndarray[ndim=2,dtype=np.float32]): tip tilt influence function
481 f = h5py.File(filename,
'r')
482 IF = csr_matrix((f[
"IF.data"][:], f[
"IF.indices"][:], f[
"IF.indptr"][:]))
483 if (list(f.keys()).count(
"TT")):
486 T = IF[-2:, :].toarray()
489 return IF, T.T.astype(np.float32)
494 Return the #n mode of the Btt modal basis contains in a ROKET file
496 filename: (str): path to the ROKET file
497 n: (int): mode number
499 sc: (np.ndarray[ndim=2,dtype=np.float32]): mode #n of the Btt basis
501 f = h5py.File(filename,
'r')
504 dim = f[
"dm_dim"].value
505 indx = f[
"indx_pup"][:]
506 sc = np.zeros((dim, dim))
508 mode = IF.T.dot(Btt[:-2, n])
509 mode += TT.T.dot(Btt[-2:, n])
512 return sc.reshape((dim, dim))
517 Return the PSF computed by COMPASS saved in the ROKET file
520 filename: (str): path to the ROKET file
522 psf: (np.ndarray[ndim=2,dtype=np.float32]): PSF computed by COMPASS
524 f = h5py.File(filename,
"r")
531 def getMap(filename, covmat):
533 Return the spatial representation of a covariance matrix expressed in the DM space
535 filename: (str): path to the ROKET file
536 covmat: (np.ndarray[ndim=2,dtype=np.float32]): covariance matrix
538 Map: (np.ndarray[ndim=2,dtype=np.float32]): covariance map
540 f = h5py.File(filename,
'r')
542 xpos = f[
"dm.xpos"][:]
543 ypos = f[
"dm.ypos"][:]
544 pitch = xpos[1] - xpos[0]
545 nact = f.attrs[
"_Param_dm__nact"][0]
546 x = ((xpos - xpos.min()) / pitch).astype(np.int32)
547 y = ((ypos - ypos.min()) / pitch).astype(np.int32)
553 xx = np.tile(np.arange(nact), (nact, 1))
555 dx = np.zeros((x.size, x.size), dtype=np.int32)
557 for k
in range(x.size):
558 dx[k, :] = xx[nn][k] - xx[nn]
559 dy[k, :] = yy[nn][k] - yy[nn]
566 Map = np.zeros((nact * 2 - 1, nact * 2 - 1)).flatten()
568 ind = dy.flatten() + (nact * 2 - 1) * (dx.flatten())
569 Cf = covmat.flatten()
570 for k
in range(ind.size):
574 div[np.where(div == 0)] = 1
577 return Map.reshape((nact * 2 - 1, nact * 2 - 1))
580 def SlopesMap(covmat, filename=None, nssp=None, validint=None):
582 Return a part of the spatial representation of a covariance matrix expressed in the slopes space.
583 Need to be called 4 times to get the full map (XX, YY, XY, YX)
586 covmat: (np.ndarray[ndim=2,dtype=np.float32]): part of the covariance matrix
587 filename: (str): (optional) path to the ROKET file
588 nssp: (int): (optional) Number of ssp in the diameter
589 validint: (float): (optional) Central obstruction as a ratio of D
592 Map: (np.ndarray[ndim=2,dtype=np.float32]): covariance map
594 if filename
is not None:
595 f = h5py.File(filename,
'r')
596 nssp = f.attrs[
"_Param_wfs__nxsub"][0]
597 validint = f.attrs[
"_Param_tel__cobs"]
600 if nssp
is None or validint
is None:
601 raise ValueError(
"nssp and validint not defined")
603 nsub = covmat.shape[0]
604 x = np.linspace(-1, 1, nssp)
605 x, y = np.meshgrid(x, x)
606 r = np.sqrt(x * x + y * y)
608 rorder = np.sort(r.reshape(nssp * nssp))
609 ncentral = nssp * nssp - np.sum(r >= validint, dtype=np.int32)
610 validext = rorder[ncentral + nsub]
611 valid = (r < validext) & (r >= validint)
614 xx = np.tile(np.arange(nssp), (nssp, 1))
618 dx = np.zeros((xx.size, xx.size), dtype=np.int32)
620 for k
in range(xx.size):
621 dx[k, :] = xx[k] - xx
622 dy[k, :] = yy[k] - yy
629 Map = np.zeros((nssp * 2 - 1, nssp * 2 - 1)).flatten()
631 ind = dy.flatten() + (nssp * 2 - 1) * (dx.flatten())
632 Cf = covmat.flatten()
633 for k
in range(ind.size):
637 div[np.where(div == 0)] = 1
640 return Map.reshape((nssp * 2 - 1, nssp * 2 - 1))
643 def covFromMap(Map, nsub, filename=None, nssp=None, validint=None):
645 Return a part of the spatial representation of a covariance matrix expressed in the slopes space.
646 Need to be called 4 times to get the full map (XX, YY, XY, YX)
649 covmat: (np.ndarray[ndim=2,dtype=np.float32]): part of the covariance matrix
650 filename: (str): (optional) path to the ROKET file
651 nssp: (int): (optional) Number of ssp in the diameter
652 validint: (float): (optional) Central obstruction as a ratio of D
655 Map: (np.ndarray[ndim=2,dtype=np.float32]): covariance map
657 if filename
is not None:
658 f = h5py.File(filename,
'r')
659 nssp = f.attrs[
"_Param_wfs__nxsub"][0]
660 validint = f.attrs[
"_Param_tel__cobs"]
663 if nssp
is None or validint
is None:
664 raise ValueError(
"nssp and validint not defined")
666 x = np.linspace(-1, 1, nssp)
667 x, y = np.meshgrid(x, x)
668 r = np.sqrt(x * x + y * y)
670 rorder = np.sort(r.reshape(nssp * nssp))
671 ncentral = nssp * nssp - np.sum(r >= validint, dtype=np.int32)
672 validext = rorder[ncentral + nsub]
673 valid = (r < validext) & (r >= validint)
676 xx = np.tile(np.arange(nssp), (nssp, 1))
680 dx = np.zeros((xx.size, xx.size), dtype=np.int32)
682 for k
in range(xx.size):
683 dx[k, :] = xx[k] - xx
684 dy[k, :] = yy[k] - yy
691 covmat = np.zeros((nsub, nsub))
692 ind = dy.flatten() + (nssp * 2 - 1) * (dx.flatten())
693 Cf = covmat.flatten()
695 for k
in range(ind.size):
698 return Cf.reshape((nsub, nsub))
701 def getCovFromMap(Map, nsub, filename=None, nssp=None, validint=None):
703 Return the full spatial representation of a covariance matrix expressed in the slopes space.
706 covmat: (np.ndarray[ndim=2,dtype=np.float32]): part of the covariance matrix
707 filename: (str): (optional) path to the ROKET file
708 nssp: (int): (optional) Number of ssp in the diameter
709 validint: (float): (optional) Central obstruction as a ratio of D
711 Map: (np.ndarray[ndim=2,dtype=np.float32]): covariance map
713 if filename
is not None:
714 f = h5py.File(filename,
'r')
715 nssp = f.attrs[
"_Param_wfs__nxsub"][0]
717 mapSize = 2 * nssp - 1
718 covmat = np.zeros((nsub, nsub))
720 covmat[:nsub // 2, :nsub // 2] =
covFromMap(Map[:mapSize, :mapSize], nsub // 2,
722 covmat[nsub // 2:, nsub // 2:] =
covFromMap(Map[mapSize:, mapSize:], nsub // 2,
724 covmat[:nsub // 2, nsub // 2:] =
covFromMap(Map[:mapSize, mapSize:], nsub // 2,
726 covmat[nsub // 2:, :nsub // 2] =
covFromMap(Map[mapSize:, :mapSize], nsub // 2,
732 def get_slopessMap(covmat, filename=None, nssp=None, validint=None):
734 Return the full spatial representation of a covariance matrix expressed in the slopes space.
737 covmat: (np.ndarray[ndim=2,dtype=np.float32]): part of the covariance matrix
738 filename: (str): (optional) path to the ROKET file
739 nssp: (int): (optional) Number of ssp in the diameter
740 validint: (float): (optional) Central obstruction as a ratio of D
742 Map: (np.ndarray[ndim=2,dtype=np.float32]): covariance map
744 if filename
is not None:
745 f = h5py.File(filename,
'r')
746 nssp = f.attrs[
"_Param_wfs__nxsub"][0]
748 nsub = covmat.shape[0] // 2
749 mapSize = 2 * nssp - 1
750 Map = np.zeros((2 * mapSize, 2 * mapSize))
752 Map[:mapSize, :mapSize] =
SlopesMap(covmat[:nsub, :nsub], filename=filename,
753 nssp=nssp, validint=validint)
754 Map[:mapSize, mapSize:] =
SlopesMap(covmat[:nsub, nsub:], filename=filename,
755 nssp=nssp, validint=validint)
756 Map[mapSize:, :mapSize] =
SlopesMap(covmat[nsub:, :nsub], filename=filename,
757 nssp=nssp, validint=validint)
758 Map[mapSize:, mapSize:] =
SlopesMap(covmat[nsub:, nsub:], filename=filename,
759 nssp=nssp, validint=validint)
764 def ensquare_PSF(filename, psf, N, display=False, cmap="jet"):
766 Return the ensquared PSF
769 filename: (str): path to the ROKET file
770 psf: (np.ndarray[ndim=2,dtype=np.float32]): PSF to ensquare
771 N: (int): size of the square in units of Lambda/D
772 display: (bool): (optional) if True, displays also the ensquare PSF
773 cmat: (str): (optional) matplotlib colormap to use
775 psf: (np.ndarray[ndim=2,dtype=np.float32]): the ensquared psf
777 f = h5py.File(filename,
'r')
778 Lambda_tar = f.attrs[
"_Param_target__Lambda"][0]
779 RASC = 180 / np.pi * 3600.
780 pixsize = Lambda_tar * 1e-6 / (psf.shape[0] * f.attrs[
"_Param_tel__diam"] / f.attrs[
781 "_Param_geom__pupdiam"]) * RASC
782 x = (np.arange(psf.shape[0]) - psf.shape[0] / 2) * pixsize / (
783 Lambda_tar * 1e-6 / f.attrs[
"_Param_tel__diam"] * RASC)
784 w = int(N * (Lambda_tar * 1e-6 / f.attrs[
"_Param_tel__diam"] * RASC) / pixsize)
785 mid = psf.shape[0] // 2
786 psfe = np.abs(psf[mid - w:mid + w, mid - w:mid + w])
788 plt.matshow(np.log10(psfe), cmap=cmap)
790 xt = np.linspace(0, psfe.shape[0] - 1, 6).astype(np.int32)
791 yt = np.linspace(-N, N, 6).astype(np.int32)
796 return psf[mid - w:mid + w, mid - w:mid + w]
801 Return the ensquared energy in a box width of N * lambda/D
804 filename: (str): path to the ROKET file
805 N: (int): size of the square in units of Lambda/D
810 def cutsPSF(filename, psf, psfs):
812 Plots cuts of two PSF along X and Y axis for comparison
814 filename: (str): path to the ROKET file
815 psf: (np.ndarray[ndim=2,dtype=np.float32]): first PSF
816 psfs: (np.ndarray[ndim=2,dtype=np.float32]): second PSF
818 f = h5py.File(filename,
'r')
819 Lambda_tar = f.attrs[
"_Param_target__Lambda"][0]
820 RASC = 180 / np.pi * 3600.
821 pixsize = Lambda_tar * 1e-6 / (psf.shape[0] * f.attrs[
"_Param_tel__diam"] / f.attrs[
822 "_Param_geom__pupdiam"]) * RASC
823 x = (np.arange(psf.shape[0]) - psf.shape[0] / 2) * pixsize / (
824 Lambda_tar * 1e-6 / f.attrs[
"_Param_tel__diam"] * RASC)
827 plt.semilogy(x, psf[psf.shape[0] // 2, :], color=
"blue")
828 plt.semilogy(x, psfs[psf.shape[0] // 2, :], color=
"red")
830 np.abs(psf[psf.shape[0] // 2, :] - psfs[psf.shape[0] // 2, :]),
832 plt.xlabel(
"X-axis angular distance [units of lambda/D]")
833 plt.ylabel(
"Normalized intensity")
834 plt.legend([
"PSF exp",
"PSF model",
"Diff"])
838 plt.semilogy(x, psf[:, psf.shape[0] // 2], color=
"blue")
839 plt.semilogy(x, psfs[:, psf.shape[0] // 2], color=
"red")
841 np.abs(psf[:, psf.shape[0] // 2] - psfs[:, psf.shape[0] // 2]),
843 plt.xlabel(
"Y-axis angular distance [units of lambda/D]")
844 plt.ylabel(
"Normalized intensity")
845 plt.legend([
"PSF exp",
"PSF model",
"Diff"])
853 Compute d/dt(slopes)*slopes from ROKET buffer
855 filename: (str): (optional) path to the ROKET file
856 slopes: (np.ndarray[ndim=2,dtype=np.float32]: (optional) Buffer of slopes arranged as (nsub x niter)
857 dt: (int): (optionnal) dt in frames
858 dd: (bool): (optionnal) if True, computes d/dt(slopes)*d/dt(slopes)
860 dCmm: (np.ndarray[ndim=2,dtype=np.float32]: covariance matrix of slopes with their derivative
862 if filename
is not None:
863 f = h5py.File(filename,
'r')
864 slopes = f[
"slopes"][:]
866 if slopes
is not None:
868 dCmm = (slopes[:, dt:] - slopes[:, :-dt]).dot(
869 (slopes[:, dt:] - slopes[:, :-dt]).T / 2)
871 dCmm = slopes[:, :-dt].dot(slopes[:, dt:].T)
873 dCmm = (slopes[:, dt:] - slopes[:, :-dt]).dot(
874 (slopes[:, dt:] + slopes[:, :-dt]).T / 2)
876 return dCmm / slopes[:, dt:].shape[1]
881 Identify turbulent parameters (wind speed, direction and frac. of r0) from ROKET file
884 filename: (str): path to the ROKET file
885 nlayers: (int): number of turbulent layers (maybe deduced in the future ?)
887 f = h5py.File(filename,
"r")
888 dt = f.attrs[
"_Param_loop__ittime"]
889 dk = int(2 / 3 * f.attrs[
"_Param_tel__diam"] / 20 / dt)
890 pdiam = f.attrs[
"_Param_tel__diam"] / f.attrs[
"_Param_wfs__nxsub"]
893 size = mapC.shape[0] // 2
894 minimap = mapC[size:, size:] + mapC[size:, :size] + mapC[:size,
895 size:] + mapC[:size, :size]
897 ws = np.zeros(nlayers)
898 wd = np.zeros(nlayers)
899 frac = np.zeros(nlayers)
901 for k
in range(nlayers):
904 x, y = np.where(minimap == minimap.max())
907 print(
"max ", k,
": x=", x,
" ; y=", y)
908 frac[k] = minimap[x, y]
909 r = np.linalg.norm([x - size / 2, y - size / 2]) * pdiam
910 ws[k] = r / (dk * dt)
911 wd[k] = np.arctan2(x - size / 2, y - size / 2) * 180 / np.pi
914 minimap[x - 2:x + 3, y - 2:y + 3] = 0
915 minimap[(size - x - 1) - 2:(size - x - 1) + 3, (size - y - 1) - 2:
916 (size - y - 1) + 3] = 0
919 ind = np.argsort(f.attrs[
"_Param_atmos__frac"])[::-1]
920 print(
"Real wind speed: ", f.attrs[
"_Param_atmos__windspeed"][ind].tolist())
921 print(
"Estimated wind speed: ", ws.tolist())
922 print(
"-----------------------------")
923 print(
"Real wind direction: ", f.attrs[
"_Param_atmos__winddir"][ind].tolist())
924 print(
"Estimated wind direction: ", wd.tolist())
925 print(
"-----------------------------")
926 print(
"Real frac: ", f.attrs[
"_Param_atmos__frac"][ind].tolist())
927 print(
"Estimated frac: ", frac.tolist())
928 print(
"-----------------------------")