40 shesha_db = os.environ[
'SHESHA_DB_ROOT']
41 except KeyError
as err:
43 shesha_db = os.environ[
'SHESHA_ROOT'] +
"/data/"
46 shesha_savepath = shesha_db
55 import scipy.ndimage.interpolation
as sci
58 def make_lgs_prof1d(p_wfs: conf.Param_wfs, p_tel: conf.Param_tel, prof: np.ndarray,
59 h: np.ndarray, beam: float, center=
""):
60 """same as prep_lgs_prof but cpu only. original routine from rico
63 p_tel: (Param_tel) : telescope settings
65 prof: (np.ndarray[dtype=np.float32]) : Na profile intensity, in arbitrary units
67 h: (np.ndarray[dtype=np.float32]) : altitude, in meters. h MUST be an array with EQUALLY spaced elements.
69 beam: (float) : size in arcsec of the laser beam
71 center: (string) : either "image" or "fourier" depending on where the centre should be.
75 p_wfs._profcum = np.zeros(h.shape[0] + 1, dtype=np.float32)
76 p_wfs._profcum[1:] = prof.cumsum()
78 subapdiam = p_tel.diam / p_wfs.nxsub
80 xsubs = np.linspace((subapdiam - p_tel.diam) / 2, (p_tel.diam - subapdiam) / 2,
81 p_wfs.nxsub).astype(np.float32)
83 xsubs = np.zeros(1, dtype=np.float32)
84 ysubs = xsubs.copy().astype(np.float32)
87 hG = np.sum(h * prof) / np.sum(prof)
88 x = np.arange(p_wfs._Ntot).astype(np.float32) - p_wfs._Ntot / 2
90 x = x * p_wfs._qpixsize
95 dOffAxis = np.sqrt((xsubs[p_wfs._validsubsy // p_wfs.npix] - p_wfs.lltx)**2 +
96 (ysubs[p_wfs._validsubsx // p_wfs.npix] - p_wfs.llty)**2)
98 dOffAxis = np.sqrt((xsubs - p_wfs.lltx)**2 + (ysubs - p_wfs.llty)**2)
100 profi = np.zeros((p_wfs._Ntot, p_wfs._nvalid), dtype=np.float32)
102 subsdone = np.ones(p_wfs._nvalid, dtype=np.int32)
103 dif2do = np.zeros(p_wfs._nvalid, dtype=np.int32)
105 while (np.any(subsdone)):
106 tmp = dOffAxis[np.where(subsdone)][0]
107 inds = np.where(dOffAxis == tmp)[0]
109 zhc = (h - hG) * (206265. * tmp / hG**2)
110 dzhc = zhc[1] - zhc[0]
112 if (p_wfs._qpixsize > dzhc):
113 avg_zhc = np.zeros(zhc.size + 1, dtype=np.float32)
115 avg_zhc[avg_zhc.size - 1] = zhc[zhc.size - 1]
116 avg_zhc[1:-1] = 0.5 * (zhc[1:] + zhc[:-1])
117 avg_x = np.zeros(x.size + 1, dtype=np.float32)
119 avg_x[avg_x.size - 1] = x[x.size - 1]
120 avg_x[1:-1] = 0.5 * (x[1:] + x[:-1])
122 for i
in range(inds.size):
123 profi[:, inds[i]] = np.diff(np.interp(avg_x, avg_zhc,
124 p_wfs._profcum)).astype(np.float32)
127 for i
in range(inds.size):
128 profi[:, inds[i]] = np.interp(x, zhc, prof)
131 w = beam / 2.35482005
135 g = np.zeros(n, dtype=np.float32)
136 if (center ==
"image"):
143 if (center ==
"image"):
144 if ((p_wfs.npix * p_wfs._nrebin) % 2 != p_wfs._Nfft % 2):
145 g = np.exp(-(x + p_wfs._qpixsize)**2 / (2 * w**2.))
147 g = np.exp(-(x + p_wfs._qpixsize / 2)**2 / (2 * w**2.))
150 g = np.exp(-x**2 / (2 * w**2.))
152 p_wfs._ftbeam = np.fft.fft(g).astype(np.complex64)
153 p_wfs._beam = g.astype(np.float32)
156 g_extended = np.tile(g, (p_wfs._nvalid, 1)).T
159 np.fft.fft(profi, axis=0) * np.fft.fft(g_extended, axis=0),
160 axis=0).real.astype(np.float32)
161 p1d = p1d * p1d.shape[0]
162 p1d = np.roll(p1d, int(p_wfs._Ntot / 2. - 0.5), axis=0)
165 im = np.zeros((p1d.shape[1], p1d.shape[0], p1d.shape[0]), dtype=np.float32)
166 for i
in range(p1d.shape[0]):
167 im[:, i, :] = g[i] * p1d.T
170 azimuth = np.arctan2(ysubs[p_wfs._validsubsy // p_wfs.npix] - p_wfs.llty,
171 xsubs[p_wfs._validsubsx // p_wfs.npix] - p_wfs.lltx)
173 azimuth = np.arctan2(ysubs - p_wfs.llty, xsubs - p_wfs.lltx)
175 p_wfs._azimuth = azimuth
177 if (center ==
"image"):
178 xcent = p_wfs._Ntot / 2. - 0.5
181 xcent = p_wfs._Ntot / 2.
189 for k
in range(im.shape[0]):
190 img = im[k, :, :] / im[k, :, :].max()
191 im[k, :, :] = sci.rotate(img, azimuth[k] * 180 / np.pi, reshape=
False)
196 im = sci.rotate(img, azimuth * 180 / np.pi, reshape=
False)
198 p_wfs._lgskern = im.T
201 def prep_lgs_prof(p_wfs: conf.Param_wfs, nsensors: int, p_tel: conf.Param_tel,
202 sensors: Sensors, center=
"", imat=0):
203 """The function returns an image array(double,n,n) of a laser beacon elongated by perpective
204 effect. It is obtaind by convolution of a gaussian of width "lgsWidth" arcseconds, with the
205 line of the sodium profile "prof". The altitude of the profile is the array "h".
208 p_wfs: (Param_wfs) : WFS settings
210 nsensors: (int) : wfs index
212 p_tel: (Param_tel) : telescope settings
214 Sensors: (Sensors) : WFS object
216 center: (string) : either "image" or "fourier" depending on where the centre should be.
218 Computation of LGS spot from the sodium profile:
219 Everything is done here in 1D, because the Na profile is the result of the convolution of a function
220 P(x,y) = profile(x) . dirac(y)
221 by a gaussian function, for which variables x and y can be split :
222 exp(-(x^2+y^2)/2.s^2) = exp(-x^2/2.s^2) * exp(-y^2/2.s^2)
223 The convolution is (symbol $ denotes integral)
224 C(X,Y) = $$ exp(-x^2/2.s^2) * exp(-y^2/2.s^2) * profile(x-X) * dirac(y-Y) dx dy
225 First one performs the integration along y
226 C(X,Y) = exp(-Y^2/2.s^2) $ exp(-x^2/2.s^2) * profile(x-X) dx
227 which shows that the profile can be computed by
228 - convolving the 1-D profile
229 - multiplying it in the 2nd dimension by a gaussian function
231 If one has to undersample the inital profile, then some structures may be "lost". In this case,
232 it's better to try to "save" those structures by re-sampling the integral of the profile, and
233 then derivating it afterwards.
234 Now, if the initial profile is a coarse one, and that one has to oversample it, then a
235 simple re-sampling of the profile is adequate.
237 if (p_wfs.proftype
is None or p_wfs.proftype ==
""):
238 p_wfs.set_proftype(scons.ProfType.GAUSS1)
240 profilename = scons.ProfType.FILES[p_wfs.proftype]
242 profile_path = shesha_savepath + profilename
243 print(
"reading Na profile from", profile_path)
244 prof = np.load(profile_path)
245 make_lgs_prof1d(p_wfs, p_tel, np.mean(prof[1:, :], axis=0), prof[0, :],
246 p_wfs.beamsize, center=
"image")
247 p_wfs.set_altna(prof[0, :].astype(np.float32))
248 p_wfs.set_profna(np.mean(prof[1:, :], axis=0).astype(np.float32))
250 p_wfs._prof1d = p_wfs._profna
251 p_wfs._profcum = np.zeros(p_wfs._profna.size + 1, dtype=np.float32)
252 p_wfs._profcum[1:] = p_wfs._profna.cumsum()
253 subapdiam = p_tel.diam / p_wfs.nxsub
255 if (p_wfs.nxsub > 1):
256 xsubs = np.linspace((subapdiam - p_tel.diam) / 2, (p_tel.diam - subapdiam) / 2,
257 p_wfs.nxsub).astype(np.float32)
259 xsubs = np.zeros(1, dtype=np.float32)
260 ysubs = xsubs.copy().astype(np.float32)
263 hG = np.sum(p_wfs._altna * p_wfs._profna) / np.sum(p_wfs._profna)
264 x = np.arange(p_wfs._Ntot).astype(np.float32) - p_wfs._Ntot / 2
267 x = x * p_wfs._qpixsize
269 dh = p_wfs._altna[1] - p_wfs._altna[0]
271 if (p_wfs.nxsub > 1):
272 dOffAxis = np.sqrt((xsubs[p_wfs._validsubsx // p_wfs.npix] - p_wfs.lltx)**2 +
273 (ysubs[p_wfs._validsubsy // p_wfs.npix] - p_wfs.llty)**2)
275 dOffAxis = np.sqrt((xsubs - p_wfs.lltx)**2 + (ysubs - p_wfs.llty)**2)
280 w = p_wfs.beamsize / 2.35482005
284 g = np.zeros(n, dtype=np.float32)
285 if (center ==
"image"):
292 if (center ==
"image"):
293 g = np.exp(-(x + p_wfs._qpixsize / 2)**2 / (2 * w**2.))
295 g = np.exp(-x**2 / (2 * w**2.))
297 p_wfs._ftbeam = np.fft.fft(g, axis=0).astype(np.complex64)
302 azimuth = np.arctan2(ysubs[p_wfs._validsubsy // p_wfs.npix] - p_wfs.llty,
303 xsubs[p_wfs._validsubsx // p_wfs.npix] - p_wfs.lltx)
305 azimuth = np.arctan2(ysubs - p_wfs.llty, xsubs - p_wfs.lltx)
307 p_wfs._azimuth = azimuth
309 sensors.d_wfs[nsensors].d_gs.d_lgs.lgs_init(
310 p_wfs._prof1d.size, hG, p_wfs._altna[0], dh, p_wfs._qpixsize, dOffAxis,
311 p_wfs._prof1d, p_wfs._profcum, p_wfs._beam, p_wfs._ftbeam, p_wfs._azimuth)