40 import scipy.ndimage.interpolation
as interp
42 from .
import hdf5_util
as h5u
43 from .
import utilities
as util
47 EELT_data = os.environ.get(
'SHESHA_ROOT') +
"/data/apertures/"
50 def make_pupil(dim, pupd, tel, xc=-1, yc=-1, real=0, halfSpider=False):
51 """Initialize the system pupil
55 dim: (long) : = p_geom.pupdiam
57 pupd: (long) : linear size of total pupil = p_geom.pupdiam
59 tel: (Param_tel) : Telescope structure
61 xc: (int) = p_geom.pupdiam / 2. - 0.5
63 yc: (int) = p_geom.pupdiam / 2. - 0.5
70 if tel.type_ap == ApertureType.EELT_NOMINAL:
73 elif (tel.type_ap == ApertureType.EELT):
75 tel.pupangle, D=tel.diam, halfSpider=halfSpider,
76 pitch=1.244683637214, nseg=33, inner_rad=4.1,
77 outer_rad=15.4, R=95.7853, nominalD=40,
78 half_seg=0.75, refl=tel.referr)
79 elif (tel.type_ap == ApertureType.KECK):
81 kpitch = seg_corner / 2 * np.sqrt(3)
89 dim, tel.t_spiders, xc, yc, tel.diam / dim, tel.gap, tel.pupangle,
90 D=tel.diam, cobs=tel.cobs, halfSpider=halfSpider, pitch=kpitch,
91 nseg=knseg, inner_rad=0.9, outer_rad=3.4, R=kR, nominalD=knominalD,
92 half_seg=0.9, refl=tel.referr, rotSpiderDegree=-30)
93 elif tel.type_ap == ApertureType.EELT_BP1:
94 print(
"ELT_pup_cobs = %5.3f" % 0.339)
97 elif tel.type_ap == ApertureType.EELT_BP3:
98 print(
"ELT_pup_cobs = %5.3f" % 0.503)
101 elif tel.type_ap == ApertureType.EELT_BP5:
102 print(
"ELT_pup_cobs = %5.3f" % 0.632)
105 elif tel.type_ap == ApertureType.EELT_CUSTOM:
106 print(
"EELT_CUSTOM not implemented. Falling back to EELT-Nominal")
107 tel.set_type_ap(ApertureType.EELT_NOMINAL)
109 elif tel.type_ap == ApertureType.VLT:
111 print(
"force_VLT_pup_cobs = %5.3f" % 0.14)
113 elif tel.type_ap == ApertureType.GENERIC:
117 raise NotImplementedError(
"Missing Pupil type.")
121 yc=0, real=0, cobs=0):
123 Initialize the system pupil
127 dim: (long) : linear size of ???
129 pupd: (long) : linear size of total pupil
131 t_spiders: (float) : secondary supports ratio.
133 spiders_type: (str) : secondary supports type: "four" or "six".
141 cobs: (float) : central obstruction ratio.
146 pup = util.dist(dim, xc, yc)
149 pup = np.exp(-(pup / (pupd * 0.5))**60.0)**0.69314
151 pup = (pup < (pupd + 1.) / 2).astype(np.float32)
155 pup -= np.exp(-(util.dist(dim, xc, yc) / (pupd * cobs * 0.5))**60.)**0.69314
157 pup -= (util.dist(dim, xc, yc) < (pupd * cobs + 1.) * 0.5).astype(np.float32)
160 first = 0.5 * (step - 1)
162 X = np.tile(np.arange(dim, dtype=np.float32) * step + first, (dim, 1))
166 t_spiders = t_spiders * pupd / dim
168 if (spiders_type ==
"four"):
170 s4_2 = 2 * np.sin(np.pi / 4)
171 t4 = np.tan(np.pi / 4)
173 spiders_map = ((X.T > (X + t_spiders / s4_2) * t4) +
174 (X.T < (X - t_spiders / s4_2) * t4)).astype(np.float32)
175 spiders_map *= ((X.T > (-X + t_spiders / s4_2) * t4) +
176 (X.T < (-X - t_spiders / s4_2) * t4)).astype(np.float32)
178 pup = pup * spiders_map
180 elif (spiders_type ==
"six"):
184 s2ma_2 = 2 * np.sin(np.pi / 2 - angle)
185 s6pa_2 = 2 * np.sin(np.pi / 6 + angle)
186 s6ma_2 = 2 * np.sin(np.pi / 6 - angle)
187 t2ma = np.tan(np.pi / 2 - angle)
188 t6pa = np.tan(np.pi / 6 + angle)
189 t6ma = np.tan(np.pi / 6 - angle)
191 spiders_map = ((X.T > (-X + t_spiders / s2ma_2) * t2ma) +
192 (X.T < (-X - t_spiders / s2ma_2) * t2ma))
193 spiders_map *= ((X.T > (X + t_spiders / s6pa_2) * t6pa) +
194 (X.T < (X - t_spiders / s6pa_2) * t6pa))
195 spiders_map *= ((X.T > (-X + t_spiders / s6ma_2) * t6ma) +
196 (X.T < (-X - t_spiders / s6ma_2) * t6ma))
197 pup = pup * spiders_map
199 print(
"Generic pupil created")
205 Initialize the VLT pupil
209 dim: (long) : linear size of ???
211 pupd: (long) : linear size of total pupil
213 tel: (Param_tel) : Telescope structure
216 if (tel.set_t_spiders == -1):
217 print(
"force t_spider =%5.3f" % (0.09 / 18.))
218 tel.set_t_spiders(0.09 / 18.)
219 angle = 50.5 * np.pi / 180.
221 Range = (0.5 * (1) - 0.25 / dim)
222 X = np.tile(np.linspace(-Range, Range, dim, dtype=np.float32), (dim, 1))
224 R = np.sqrt(X**2 + (X.T)**2)
226 pup = ((R < 0.5) & (R > (tel.cobs / 2))).astype(np.float32)
228 if (tel.set_t_spiders == -1):
233 (X - tel.cobs / 2 + tel.t_spiders / np.sin(angle)) * np.tan(angle)) +
234 (X.T < (X - tel.cobs / 2) * np.tan(angle))) * (X > 0) * (X.T > 0)
235 spiders_map += np.fliplr(spiders_map)
236 spiders_map += np.flipud(spiders_map)
237 spiders_map = interp.rotate(spiders_map, tel.pupangle, order=0, reshape=
False)
239 pup = pup * spiders_map
241 print(
"VLT pupil created")
247 Initialize the EELT pupil
251 dim: (long) : linear size of ???
253 pupd: (long) : linear size of total pupil
255 tel: (Param_tel) : Telescope structure
260 TODO : add force rescal pup elt
263 EELT_file = EELT_data +
"EELT-Custom_N" + str(dim) +
"_COBS" + str(
264 100 * tel.cobs) +
"_CLOCKED" + str(tel.pupangle) +
"_TSPIDERS" + str(
265 100 * tel.t_spiders) +
"_MS" + str(
266 tel.nbrmissing) +
"_REFERR" + str(
267 100 * tel.referr) +
".h5"
269 EELT_file = EELT_data + tel.type_ap.decode(
'UTF-8') +
"_N" + str(
270 dim) +
"_COBS" + str(100 * tel.cobs) +
"_CLOCKED" + str(
271 tel.pupangle) +
"_TSPIDERS" + str(
272 100 * tel.t_spiders) +
"_MS" + str(
273 tel.nbrmissing) +
"_REFERR" + str(
274 100 * tel.referr) +
".h5"
275 if (os.path.isfile(EELT_file)):
276 print(
"reading EELT pupil from file ", EELT_file)
277 pup = h5u.readHdf5SingleDataset(EELT_file)
279 print(
"creating EELT pupil...")
280 file = EELT_data +
"Coord_" + tel.type_ap.decode(
'UTF-8') +
".dat"
281 data = np.fromfile(file, sep=
"\n")
282 data = np.reshape(data, (data.size // 2, 2))
286 file = EELT_data +
"EELT_MISSING_" + tel.type_ap.decode(
'UTF-8') +
".dat"
287 k_seg = np.fromfile(file, sep=
"\n").astype(np.int32)
289 W = 1.45 * np.cos(np.pi / 6)
293 Range = (0.5 * (tel.diam * dim / pupd) - 0.25 / dim)
294 X = np.tile(np.linspace(-Range, Range, dim, dtype=np.float32), (dim, 1))
296 if (tel.t_spiders == -1):
297 print(
"force t_spider =%5.3f" % (0.014))
298 tel.set_t_spiders(0.014)
302 if (tel.nbrmissing > 0):
303 k_seg = np.sort(k_seg[:tel.nbrmissing])
305 file = EELT_data +
"EELT_REF_ERROR" +
".dat"
306 ref_err = np.fromfile(file, sep=
"\n")
311 std_ref = np.std(ref_err)
313 ref_err = ref_err * tel.referr / std_ref
315 if (tel.nbrmissing > 0):
318 pup = np.zeros((dim, dim))
320 t_3 = np.tan(np.pi / 3.)
323 vect_seg = tel.vect_seg
327 pup+=(1-ref_err[i])*(Yt<0.5*W)*(Yt>=-0.5*W)*(0.5*(Yt+t_3*Xt)<0.5*W) \
328 *(0.5*(Yt+t_3*Xt)>=-0.5*W)*(0.5*(Yt-t_3*Xt)<0.5*W) \
329 *(0.5*(Yt-t_3*Xt)>=-0.5*W)
332 for i
in range(N_seg):
335 pup+=(1-ref_err[i])*(Yt<0.5*W)*(Yt>=-0.5*W)*(0.5*(Yt+t_3*Xt)<0.5*W) \
336 *(0.5*(Yt+t_3*Xt)>=-0.5*W)*(0.5*(Yt-t_3*Xt)<0.5*W) \
337 *(0.5*(Yt-t_3*Xt)>=-0.5*W)
338 if (tel.t_spiders == 0):
341 t_spiders = tel.t_spiders * (tel.diam * dim / pupd)
343 s2_6 = 2 * np.sin(np.pi / 6)
344 t_6 = np.tan(np.pi / 6)
346 spiders_map = np.abs(X) > t_spiders / 2
349 (X + t_spiders / s2_6) * t_6) + (X.T <
350 (X - t_spiders / s2_6) * t_6))
353 (-X + t_spiders / s2_6) * t_6) + (X.T <
354 (-X - t_spiders / s2_6) * t_6))
356 pup = pup * spiders_map
358 if (tel.pupangle != 0):
359 pup = interp.rotate(pup, tel.pupangle, reshape=
False, order=0)
361 print(
"writing EELT pupil to file ", EELT_file)
362 h5u.writeHdf5SingleDataset(EELT_file, pup)
364 print(
"EELT pupil created")
368 def make_phase_ab(dim, pupd, tel, pup=None, xc=-1, yc=-1, real=0, halfSpider=False):
369 """Compute the EELT M1 phase aberration
373 dim: (long) : linear size of ???
375 pupd: (long) : linear size of total pupil
377 tel: (Param_tel) : Telescope structure
384 if ((tel.type_ap == ApertureType.GENERIC)
or (tel.type_ap == ApertureType.VLT)):
385 return np.zeros((dim, dim)).astype(np.float32)
387 elif tel.type_ap == ApertureType.EELT:
390 dim, 0, xc, yc, tel.diam / dim, tel.gap, tel.pupangle, D=tel.diam,
391 halfSpider=halfSpider, pitch=1.244683637214, nseg=33, inner_rad=4.1,
392 outer_rad=15.4, R=95.7853, nominalD=40, half_seg=0.75,
393 refl=[tel.std_piston, tel.std_tt, tel.std_tt])
394 elif (tel.type_ap == ApertureType.KECK):
396 kpitch = seg_corner / 2 * np.sqrt(3)
404 dim, 0, xc, yc, tel.diam / dim, tel.gap, tel.pupangle, D=tel.diam,
405 cobs=tel.cobs, halfSpider=halfSpider, pitch=kpitch, nseg=knseg,
406 inner_rad=0.9, outer_rad=3.4, R=kR, nominalD=knominalD, half_seg=0.9,
407 refl=[tel.std_piston, tel.std_tt, tel.std_tt])
409 ab_file = EELT_data +
"aberration_" + tel.type_ap.decode(
'UTF-8') + \
410 "_N" + str(dim) +
"_NPUP" + str(np.where(pup)[0].size) +
"_CLOCKED" + str(
411 tel.pupangle) +
"_TSPIDERS" + str(
412 100 * tel.t_spiders) +
"_MS" + str(tel.nbrmissing) +
"_REFERR" + str(
413 100 * tel.referr) +
"_PIS" + str(
414 tel.std_piston) +
"_TT" + str(tel.std_tt) +
".h5"
415 if (os.path.isfile(ab_file)):
416 print(
"reading aberration phase from file ", ab_file)
417 phase_error = h5u.readHdf5SingleDataset(ab_file)
419 print(
"computing M1 phase aberration...")
421 std_piston = tel.std_piston
424 W = 1.45 * np.cos(np.pi / 6)
426 file = EELT_data +
"EELT_Piston_" + tel.type_ap.decode(
'UTF-8') +
".dat"
427 p_seg = np.fromfile(file, sep=
"\n")
429 std_pis = np.std(p_seg)
430 p_seg = p_seg * std_piston / std_pis
433 file = EELT_data +
"EELT_TT_" + tel.type_ap.decode(
'UTF-8') +
".dat"
434 tt_seg = np.fromfile(file, sep=
"\n")
436 file = EELT_data +
"EELT_TT_DIRECTION_" + tel.type_ap.decode(
'UTF-8') +
".dat"
437 tt_phi_seg = np.fromfile(file, sep=
"\n")
439 phase_error = np.zeros((dim, dim))
440 phase_tt = np.zeros((dim, dim))
441 phase_defoc = np.zeros((dim, dim))
443 file = EELT_data +
"Coord_" + tel.type_ap.decode(
'UTF-8') +
".dat"
444 data = np.fromfile(file, sep=
"\n")
445 data = np.reshape(data, (data.size // 2, 2))
449 Range = (0.5 * (tel.diam * dim / pupd) - 0.25 / dim)
451 X = np.tile(np.linspace(-Range, Range, dim, dtype=np.float32), (dim, 1))
452 t_3 = np.tan(np.pi / 3.)
454 for i
in range(N_seg):
457 SEG=(Yt<0.5*W)*(Yt>=-0.5*W)*(0.5*(Yt+t_3*Xt)<0.5*W) \
458 *(0.5*(Yt+t_3*Xt)>=-0.5*W)*(0.5*(Yt-t_3*Xt)<0.5*W) \
459 *(0.5*(Yt-t_3*Xt)>=-0.5*W)
462 N_in_seg = np.sum(SEG)
463 Hex_diam = 2 * np.max(
464 np.sqrt(Xt[np.where(SEG)]**2 + Yt[np.where(SEG)]**2))
468 np.cos(tt_phi_seg[i]) * Xt + np.sin(tt_phi_seg[i]) * Yt)
469 mean_tt = np.sum(TT[np.where(SEG == 1)]) / N_in_seg
470 phase_tt += SEG * (TT - mean_tt)
474 phase_error += SEG * p_seg[i]
476 N_EELT = np.where(pup)[0].size
477 if (np.sum(phase_tt) != 0):
478 phase_tt *= std_tt / np.sqrt(
479 1. / N_EELT * np.sum(phase_tt[np.where(pup)]**2))
483 phase_error += phase_tt + phase_defoc
485 if (tel.pupangle != 0):
486 phase_error = interp.rotate(phase_error, tel.pupangle, reshape=
False,
489 print(
"phase aberration created")
490 print(
"writing aberration filel to file ", ab_file)
491 h5u.writeHdf5SingleDataset(ab_file, phase_error)
498 _____ _ _____ ____ ___ ____ ___
499 | ____| | |_ _| | _ \|_ _/ ___/ _ \
500 | _| | | | | | |_) || | | | | | |
501 | |___| |___| | | _ < | | |__| |_| |
502 |_____|_____|_| |_| \_\___\____\___/
509 centerMark=0, halfSpider=False, pitch=1.244683637214, nseg=33,
510 inner_rad=4.1, outer_rad=15.4, R=95.7853, nominalD=40,
511 half_seg=0.75, refl=None, rotSpiderDegree=None):
513 Generates a boolean pupil mask of the binary EELT pupil
514 on a map of size (npt, npt).
517 :returns: pupil image (npt, npt), boolean
518 :param int npt: size of the output array
519 :param float dspider: width of spiders in meters
520 :param float i0, j0: index of pixels where the pupil should be centred = p_geom.pupdiam / 2. - 0.5
521 Can be floating-point indexes.
522 :param float pixscale: size of a pixel of the image, in meters = ptel.diam/(p_geom.pupdiam / 2. - 0.5)
523 :param float gap: gap between 2 segments in metres
524 :param float rotdegree: rotation angle of the pupil, in degrees.
525 :param float D: diameter of the pupil. For the nominal EELT, D shall
527 :param int centerMark: when centerMark!=0, a pixel is added at the centre of
528 symmetry of the pupil in order to debug things using compass.
529 centerMark==1 draws a point
530 centerMark==2 draws 2 lines
531 :param bool halfSpider: half Spider computation flag
532 :param float pitch: segment pitch
533 :param int nseg: number of segments across the diameter
534 :param float inner_rad: Inner radius [meters]
535 :param float outter_rad: outter radius [meters]
536 :param float R: M1 curvature radius
537 :param float nominalD: diameter needed to get nominal aperture after projection
538 :param float half_seg: segment half size
539 :param float refl: std of the reflectivity of each segment
547 pixscale = D/(npt / 2. - 0.5)
550 pup = generateEeltPupilMask(npt, dspider, i0, j0, pixscale, gap, rotdegree)
553 rot = rotdegree * np.pi / 180
555 if rotSpiderDegree
is None:
558 rotSpider = rotSpiderDegree * np.pi / 180
565 outer_rad=outer_rad, R=R, nominalD=nominalD)
571 elif np.isscalar(refl):
572 referr = np.random.random(hx.size)
573 referr = referr * refl / np.std(referr)
574 refl = np.ones(hx.size) - referr
575 elif type(refl) == list:
577 refpist = np.random.random(hx.size)
578 refpist = refpist * refl[0] / np.std(refpist)
579 reftip = np.random.random(hx.size)
580 reftip = reftip * refl[1] / np.std(reftip)
581 reftilt = np.random.random(hx.size)
582 reftilt = reftilt * refl[2] / np.std(reftilt)
583 refl = np.array([refpist, reftip, reftilt])
586 "refl param must be None, scalar (reflectivity std error) or list of 3 elements (piston, tip and tilt std errors)"
590 nominalD=nominalD, pitch=pitch, half_seg=half_seg)
593 if (dspider > 0
and nspider > 0):
594 if (halfSpider
is True):
595 pup = pup *
fillHalfSpider(npt, nspider, dspider, i0, j0, pixscale,
598 pup = pup *
fillSpider(npt, nspider, dspider, i0, j0, pixscale, rotSpider)
606 obstru = (util.dist(pup.shape[0], pup.shape[0] // 2 + 0.5,
607 pup.shape[0] // 2 + 0.5) >=
608 (pup.shape[0] * cobs + 1.) * 0.5).astype(np.float32)
613 def fillPolygon(x, y, i0, j0, scale, gap, N, index=0):
615 From a list of points defined by their 2 coordinates list
616 x and y, creates a filled polygon with sides joining the points.
617 The polygon is created in an image of size (N, N).
618 The origin (x,y)=(0,0) is mapped at pixel i0, j0 (both can be
619 floating-point values).
620 Arrays x and y are supposed to be in unit U, and scale is the
621 pixel size in U units.
623 :returns: filled polygon (N, N), boolean
624 :param float x, y: list of points defining the polygon
625 :param float i0, j0: index of pixels where the pupil should be centred.
626 Can be floating-point indexes.
627 :param float scale: size of a pixel of the image, in same unit as x and y.
628 :param float N: size of output image.
631 x = np.array([1,-1,-1.5,0,1.1])
632 y = np.array([1,1.5,-0.2,-2,0])
638 pol = fillPolygon(x, y, i0, j0, scale, gap, N, index=2)
642 X = (np.arange(N) - i0) * scale
643 Y = (np.arange(N) - j0) * scale
644 X, Y = np.meshgrid(X, Y, indexing=
'ij')
652 T = (np.arctan2(Y - y0, X - x0) + 2 * np.pi) % (2 * np.pi)
653 t = (np.arctan2(y - y0, x - x0) + 2 * np.pi) % (2 * np.pi)
659 sens = np.median(np.diff(t))
669 x = np.roll(x, -imin)
670 y = np.roll(y, -imin)
671 t = np.roll(t, -imin)
678 indx, indy = (np.array([], dtype=np.int), np.array([], dtype=np.int))
679 distedge = np.array([], dtype=np.float)
684 sub = np.where((T >= t[-1]) | (T <= (t[0])))
686 sub = np.where((T >= t[i]) & (T <= t[j]))
690 vnorm = np.sqrt(dx**2 + dy**2)
694 crossprod = dx * (Y[sub] - y[i]) - dy * (X[sub] - x[i])
695 tmp = crossprod > gap
696 indx = np.append(indx, sub[0][tmp])
697 indy = np.append(indy, sub[1][tmp])
698 distedge = np.append(distedge, crossprod[tmp])
703 return (indx, indy, distedge)
706 a[indx, indy] = distedge
709 a = np.zeros((N, N), dtype=np.bool)
718 i0 = j0 = N / 2. - 0.5
723 Utilisee dans compass/shesha/shesha/supervisor/canapassSupervisor.py
724 pour le slaving des actus.
727 msk = np.zeros((6, pupNoSpiders.shape[0], pupNoSpiders.shape[1]))
728 mid = int(pupNoSpiders.shape[0] / 2)
729 tmp = np.zeros((pupNoSpiders.shape[0], pupNoSpiders.shape[1]))
730 for angle
in np.linspace(0, 60, npts):
731 tmp +=
compute1Spider(0, N, dspider, i0, j0, pixscale, angle * 2 * np.pi / 360)
732 msk[5, :, :] = tmp < npts * pupNoSpiders
734 msk[2, :, :] = tmp < npts * pupNoSpiders
737 for angle
in np.linspace(0, 60, npts):
738 tmp +=
compute1Spider(1, N, dspider, i0, j0, pixscale, angle * 2 * np.pi / 360)
739 msk[0, :, :] = tmp < npts * pupNoSpiders
741 msk[3, :, :] = tmp < npts * pupNoSpiders
744 for angle
in np.linspace(0, 60, npts):
745 tmp +=
compute1Spider(2, N, dspider, i0, j0, pixscale, angle * 2 * np.pi / 360)
746 msk[1, :, :] = tmp < npts * pupNoSpiders
748 msk[4, :, :] = tmp < npts * pupNoSpiders
755 Fonction de fab pour creer le slaving.
756 La fonction cree un tableau de booleens avec une seule spider.
757 Utilisee par la fonction compute6Segments()
759 a = np.ones((N, N), dtype=np.bool)
760 X = (np.arange(N) - i0) * scale
761 Y = (np.arange(N) - j0) * scale
762 X, Y = np.meshgrid(X, Y, indexing=
'ij')
765 nn = (abs(X * np.cos(i * w - rot) + Y * np.sin(i * w - rot)) < dspider / 2.)
770 def fillSpider(N, nspider, dspider, i0, j0, scale, rot):
772 Creates a boolean spider mask on a map of dimensions (N,N)
773 The spider is centred at floating-point coords (i0,j0).
775 :returns: spider image (boolean)
776 :param int N: size of output image
777 :param int nspider: number of spiders
778 :param float dspider: width of spiders
779 :param float i0: coord of spiders symmetry centre
780 :param float j0: coord of spiders symmetry centre
781 :param float scale: size of a pixel in same unit as dspider
782 :param float rot: rotation angle in radians
785 a = np.ones((N, N), dtype=np.bool)
786 X = (np.arange(N) - i0) * scale
787 Y = (np.arange(N) - j0) * scale
788 X, Y = np.meshgrid(X, Y, indexing=
'ij')
789 w = 2 * np.pi / nspider
790 for i
in range(nspider):
791 nn = (abs(X * np.cos(i * w - rot) + Y * np.sin(i * w - rot)) < dspider / 2.)
797 a = np.ones((N, N), dtype=np.bool)
798 b = np.ones((N, N), dtype=np.bool)
799 X = (np.arange(N) - i0) * scale
800 Y = (np.arange(N) - j0) * scale
801 X, Y = np.meshgrid(X, Y, indexing=
'ij')
802 w = 2 * np.pi / nspider
803 for i
in range(nspider):
804 right = (X * np.cos(i * w - rot) + Y * np.sin(i * w - rot) <
805 dspider / 2) * (X * np.cos(i * w - rot) + Y * np.sin(i * w - rot) > 0.)
806 left = (X * np.cos(i * w - rot) + Y * np.sin(i * w - rot) >
807 -dspider / 2) * (X * np.cos(i * w - rot) + Y * np.sin(i * w - rot) < 0.)
815 Cree une liste de coordonnees qui decrit un maillage hexagonal.
816 Retourne un tuple (x,y).
818 Le maillage est centre sur 0, l'un des points est (0,0).
819 Une des pointes de l'hexagone est dirigee selon l'axe Y, au sens ou le
820 tuple de sortie est (x,y).
822 :param float pitch: distance between 2 neighbour points
823 :param int supportSize: size of the support that need to be populated
827 nx = int(np.ceil((supportSize / 2.0) / pitch) + 1)
828 x = pitch * (np.arange(2 * nx + 1) - nx)
829 ny = int(np.ceil((supportSize / 2.0) / pitch / V3) + 1)
830 y = (V3 * pitch) * (np.arange(2 * ny + 1) - ny)
831 x, y = np.meshgrid(x, y, indexing=
'ij')
834 peak_axis = np.append(x, x + pitch / 2.)
835 flat_axis = np.append(y, y + pitch * V3 / 2.)
836 return flat_axis, peak_axis
840 outer_rad=15.4, R=95.7853, nominalD=40):
842 Computes the coordinates of the corners of all the hexagonal
844 Result is a tuple of arrays(6, 798).
847 -----------------------------------------
848 D: (float) : pupil diameter in meters (it must be set to 40.0 m for the ELT)
849 rot: (float) : pupil rotation angle in radians
850 pitch: (float): Segment pitch [meters]
851 nseg: (int) : number of segments across the diameter
852 inner_rad : (float): Inner radius [meters]
853 outer_rad : (float): Outer radius [meters]
854 R : (float): Curvature radius of the M1
855 nominalD: (float): diameter for nominal pupil
866 ll = np.sqrt(lx**2 + ly**2)
870 nn = (ll > inner_rad * pitch) & (ll < outer_rad * pitch)
874 ll = np.sqrt(lx**2 + ly**2)
880 th = np.linspace(0, 2 * np.pi, 7)[0:6]
881 hx = np.cos(th) * pitch / V3
882 hy = np.sin(th) * pitch / V3
887 x = (lx[
None, :] + hx[:,
None])
888 y = (ly[
None, :] + hy[:,
None])
889 r = np.sqrt(x**2 + y**2)
891 rrc = R / r * np.arctan(r / R)
901 mrot = np.array([[np.cos(rot), np.sin(rot)], [-np.sin(rot), np.cos(rot)]])
907 xyrot = np.dot(mrot, np.transpose(np.array([x, y]), (1, 0, 2)))
909 return xyrot[0], xyrot[1]
914 La fonction est appelee quand l'utilisateur a demande une pupille
915 ELT, et renseigne un diametre de telescope different de 40 metres.
917 Faut vraiment que je commente ou t'as compris ??
921 "\n\n\n\n",
"__ ___ ____ _ _ ___ _ _ ___ _",
922 "\ \ / / \ | _ \| \ | |_ _| \ | |/ ___|",
923 " \ \ /\ / / _ \ | |_) | \| || || \| | | _ ",
924 " \ V V / ___ \| _ <| |\ || || |\ | |_| |",
925 " \_/\_/_/ \_\_| \_\_| \_|___|_| \_|\____|",
" \n",
926 "Vous utilisez un telescope de type ELT. Ce telescope",
927 "est fait pour etre utilise avec un diametre de 40 m.",
" ",
928 "Or, vous utilisez un diametre different. Cela signifie",
929 "que le telescope que vous etes en train de creer a une",
930 "taille differente du veritable E-ELT de l'ESO.",
" ",
931 " * Soit vous savez exactement ce que vous faites, auquel",
932 "cas bonne route.",
" ",
933 " * Soit vous desirez creer LE vrai E-ELT et il faut changer",
935 " 1) le diametre telescope de votre fichier de parametres et",
936 " le renseigner a 40 metres.",
937 " p_tel.set_diam(40.0) # Nominal size for the real EELT",
938 " 2) le nombre d'actionneurs de M4 a 75",
939 " p_dm0.set_nact(75) # 75 actu in 40m for pitch=54.05cm",
940 " 3) option: tourner la pupille de 90 degres pour revenir au",
941 " cas initial de compass",
942 " p_tel.set_pupangle(90.) # ELT pup rotation in degrees"
951 Reorganisation des segments facon ESO.
953 ESO-193058 Standard Coordinate System and Basic Conventions
955 :param float x: tableau des centres X des segments
956 :param float y: idem Y
957 :return tuple (x,y): meme tuple que les arguments d'entree, mais tries.
965 t = (np.arctan2(y, x) + pi_6 - 1e-3) % (pix2)
970 sector = (t > k * pi_3) & (t < (k + 1) * pi_3)
972 distance = (A * np.cos(u) - np.sin(u)) * x[sector] + (
973 np.cos(u) + A * np.sin(u)) * y[sector]
974 indsort = np.argsort(distance)
975 X = np.append(X, x[sector][indsort])
976 Y = np.append(Y, y[sector][indsort])
982 Returns the data type of a numpy variable, either scalar value or array
984 if np.isscalar(truc):
987 return type(truc.flatten()[0])
990 def generateSegmentProperties(attribute, hx, hy, i0, j0, scale, gap, N, D, softGap=0,
991 nominalD=40, pitch=1.244683637214, half_seg=0.75):
993 Builds a 2D image of the pupil with some attributes for each of the
994 segments. Those segments are described from arguments hx and hy, that
995 are produced by the function generateCoordSegments(D, rot).
997 When attribute is a phase, then it must be a float array of dimension
998 [3, 798] with the dimension 3 being piston, tip, and tilt.
999 Units of phase is xxx rms, and the output of the procedure will be
1003 :returns: pupil image (N, N), with the same type of input argument attribute
1005 :param float/int/bool attribute: scalar value or 1D-array of the reflectivity of
1006 the segments or 2D array of phase
1007 If attribute is scalar, the value will be replicated for all segments.
1008 If attribute is a 1D array, then it shall contain the reflectivities
1010 If attribute is a 2D array then it shall contain the piston, tip
1011 and tilt of the segments. The array shall be of dimension
1012 [3, 798] that contains [piston, tip, tilt]
1013 On output, the data type of the pupil map will be the same as attribute.
1014 :param float hx, hy: arrays [6,:] describing the segment shapes. They are
1015 generated using generateCoordSegments()
1016 :param float dspider: width of spiders in meters
1017 :param float i0, j0: index of pixels where the pupil should be centred.
1018 Can be floating-point indexes.
1019 :param float scale: size of a pixel of the image, in meters.
1020 :param float gap: half-space between segments in meters
1021 :param int N: size of the output array (N,N)
1022 :param float D: diameter of the pupil. For the nominal EELT, D shall
1024 :param bool softGap: if False, the gap between segments is binary 0/1
1025 depending if the pixel is within the gap or not. If True, the gap
1026 is a smooth region of a fwhm of 2 pixels with a depth related to the
1028 :param float nominalD: diameter needed to get nominal pupil aperture
1029 :param float pitch: segment pitch
1030 :param float half_seg: segment half size
1034 attribute = np.ones(798)+np.random.randn(798)/20.
1048 if np.isscalar(attribute):
1049 attribute = np.array([attribute] * nseg)
1052 pupil = np.zeros((N, N), dtype=
getdatatype(attribute))
1055 x0 = np.mean(hx, axis=0)
1056 y0 = np.mean(hy, axis=0)
1058 x0 = x0 / scale + i0
1059 y0 = y0 / scale + j0
1061 hexrad = half_seg * D / nominalD / scale
1062 ix0 = np.floor(x0 - hexrad).astype(int) - 1
1063 iy0 = np.floor(y0 - hexrad).astype(int) - 1
1064 segdiam = np.ceil(hexrad * 2 + 1).astype(int) + 1
1066 n = attribute.shape[0]
1080 for i
in range(nseg):
1081 indx, indy, distedge =
fillPolygon(hx[:, i], hy[:, i], i0 - ix0[i],
1082 j0 - iy0[i], scale, gap * 0., segdiam,
1084 pupil[indx + ix0[i], indy + iy0[i]] = attribute[i] * (
1085 1. - (gap / scale / np.pi) / (1 + (distedge / scale)**2))
1088 for i
in range(nseg):
1089 indx, indy, distedge =
fillPolygon(hx[:, i], hy[:, i], i0 - ix0[i],
1090 j0 - iy0[i], scale, gap, segdiam,
1092 pupil[indx + ix0[i], indy + iy0[i]] = attribute[i]
1095 minimap = np.zeros((segdiam, segdiam))
1096 xmap = np.arange(segdiam) - segdiam / 2
1097 xmap, ymap = np.meshgrid(xmap, xmap, indexing=
'ij')
1099 diamseg = pitch * 2 / np.sqrt(3)
1100 diamfrizou = (pitch + diamseg) / 2 * D / nominalD
1105 factunit = 4 * scale / diamfrizou
1106 for i
in range(nseg):
1107 indx, indy, _ =
fillPolygon(hx[:, i], hy[:, i], i0 - ix0[i], j0 - iy0[i],
1108 scale, 0., segdiam, index=1)
1109 minimap = attribute[0, i] + (factunit * attribute[1, i]) * xmap + (
1110 factunit * attribute[2, i]) * ymap
1111 pupil[indx + ix0[i], indy + iy0[i]] = minimap[indx, indy]
1118 Renvoie une image de boolens (False) de taille (N,N) avec un point
1119 ou une croix (True) centree sur (i0, j0).
1120 :param int N: taille de l'image de sortie
1121 :param float i0, j0: position du marqueur de sortie
1122 :param int centerMark: 0 (pour rien), 1 (option point) ou 2 (option croix)
1126 X = (np.arange(N) - i0) * scale
1127 Y = (np.arange(N) - j0) * scale
1128 X, Y = np.meshgrid(X, Y, indexing=
'ij')
1130 res = (X**2 + Y**2) < 1
1132 res = (np.abs(X) < 0.9) | (np.abs(Y) < 0.9)