diff --git a/divtel/array_old.py b/divtel/array_old.py deleted file mode 100644 index 072792f..0000000 --- a/divtel/array_old.py +++ /dev/null @@ -1,522 +0,0 @@ -import numpy as np -import astropy.units as u -import matplotlib.pyplot as plt -import ipywidgets as widgets - -from astropy.table import Table -import astropy.table as tab - -from . import utils as utils -from . import visualization as visual -from . import pointing - -from .cta import CTA_Info - -from astropy.coordinates import SkyCoord - -from shapely.ops import unary_union, polygonize -from shapely.geometry import LineString, Point - -import copy - - -class Array: - - """ - Class containing information on array (a set of telescopes) - - Parameters - ---------- - telescope_list: array, class.Telescope - Array containing class.Telescope - frame: class.CTA_Info - CTA Information - kwargs: dict - args for class.CTA_Info - """ - - def __init__(self, telescope_list, frame=None, pointing2src=False, **kwargs): - - self.telescopes = telescope_list - - self._div = 0 - self._pointing = {"az":0*u.deg, "alt":0*u.deg, "ra": 0*u.deg, "dec": 0*u.deg} - - if frame == None: - self._frame = CTA_Info(verbose=False, **kwargs) - else: - self._frame = copy.copy(frame) - if pointing2src and (self.frame.source is not None): - self.set_pointing_coord(ra = self.frame.source.icrs.ra.deg, - dec = self.frame.source.icrs.dec.deg) - - self.__make_table__() - - def __make_table__(self): - - """ - Merge rows from Telescope.table - - Returns - ------- - astropy.table - """ - - table = [] - - for tel in self.telescopes: - - table.append(tel.table) - - units = 'rad' - if hasattr(self, "_table"): - if hasattr(self._table, "units"): - units = self.table.units - - self._table = tab.vstack(table) - - self._table.add_column(self._dist2tel(), name="d_tel") - self._table["d_tel"].unit = u.m - self._table["d_tel"].info.format = "{:.2f}" - - self._table.units = units - - def __convert_units__(self, toDeg=True): - """ - Convert units in the table from deg to rad, and vice versa. - - Parameters - ---------- - toDeg: bool, optional - """ - - self._table = utils.deg2rad(self._table, toDeg) - if toDeg: - self._table.units = 'deg' - else: - self._table.units = 'rad' - - def _dist2tel(self): - """ - Distance to the telescope from the barycenter - - Returns - ------- - float - """ - - dist = np.zeros(self.size_of_array) - for i, axis in enumerate(["x", "y", "z"]): - dist += (self.table[axis] - self.barycenter[i])**2. - dist = np.sqrt(dist) - return dist - - @property - def table(self): - """ - Return astropy table containing all telescope information - An attribute, Array.table.units, determines the units of `az`, `alt`, `zn`, - `radius`, `fov` units (deg or rad) - - Returns - ------- - astropy.table - """ - if hasattr(self._table, "units"): - if (self._table.units == 'deg')*(self._table["az"].unit == u.rad): - self._table = utils.deg2rad(self._table, True) - elif (self._table.units == 'rad')*(self._table["az"].unit == u.deg): - self._table = utils.deg2rad(self._table, False) - return self._table - - @property - def size_of_array(self): - """ - Return the number of telescopes - - Returns - ------- - float - """ - return self._table.__len__() - - @property - def frame(self): - """ - Return a frame - - Returns - ------- - class.CTA_info - """ - return self._frame - - @property - def barycenter(self): - """ - Return a barycenter - - Returns - ------- - array - [b_x, b_y, b_z] - """ - return np.array(utils.calc_mean(self.table, ["x", "y", "z"])) - - @property - def div(self): - """ - Return a divergence parameter - - Returns - ------- - float - """ - return self._div - - @property - def pointing(self): - """ - Return pointing information - - Returns - ------- - dict - keys: `ra`, `dec`, `az`, and `alt` - """ - return self._pointing - - def calc_mean(self, params): - """ - Calculate the mean values of parameters - - Parameters - ---------- - params: str or array - see ArrayConfig.utils.calc_mean - - Returns - ------- - array - """ - return np.array(utils.calc_mean(self.table, params)) - - def hFoV(self, m_cut=0, return_multiplicity=False, full_output=False): - """ - Return a hyper field of view (hFoV) above a given multiplicity. - - Parameters - ---------- - m_cut: float, optional - the minimum multiplicity - return_multiplicity: bool, optional - return average and variance of multiplicity - full_output: bool, optional - return all parameters; multiplicity, overlaps, geoms - Returns - ------- - fov: float - hFoV - m_ave: float - average of multiplicity - m_var: float - variance of multiplicity - multiplicity: array - array containing multiplicity and corresponding hFoV - overlaps: array - array containing the number of overlaps for each patch - geoms: shapely.ops.polygonize - geometry of each patch - """ - if self.table.units == 'rad': - self.__convert_units__(toDeg=True) - - coord = self.get_pointing_coord(icrs=False) - - if max(self.table["az"])-min(self.table["az"]) > 180: - polygons = [] - for az, alt, r in zip(coord.az.degree, coord.alt.degree, self.table["radius"]): - if az < 180: - polygons.append(Point(az, alt).buffer(r)) - else: - polygons.append(Point(az-360, alt).buffer(r)) - else: - polygons = [Point(az, alt).buffer(r) for az, alt, r in zip(coord.az.degree, coord.alt.degree, self.table["radius"])] - - union = unary_union([LineString(list(pol.exterior.coords)) for pol in polygons]) - geoms = [geom for geom in polygonize(union)] - hfov = [geom.area for geom in geoms] - - count_overlaps = np.zeros(len(geoms)) - for i, geom in enumerate(geoms): - count_overlaps[i] = sum([1 for pol in polygons if abs(geom.difference(pol).area)<1e-5]) - - hfov = np.array(hfov) - - # multiplicity associated with each patch - overlaps = np.array(count_overlaps) - eff_overlaps=[] - eff_geoms=[] - for i in range(len(overlaps)): - if overlaps[i]>m_cut: - eff_overlaps.append(overlaps[i]) - eff_geoms.append(geoms[i]) - multiplicity = np.array([[i, hfov[overlaps==i].sum()] for i in set(overlaps)]) - - fov = sum(multiplicity[:,1][multiplicity[:,0]>=m_cut])*u.deg**2 - - if full_output: - return multiplicity, eff_overlaps, eff_geoms - elif return_multiplicity: - m_ave = np.average(multiplicity[:,0], weights=multiplicity[:,1]) - m_var = np.average((multiplicity[:,0]-m_ave)**2, weights=multiplicity[:,1]) - return fov, m_ave, m_var - else: - return fov - - - def update_frame(self, site=None, time=None, delta_t=None, verbose=False): - """ - Update class.CTA_Info parameters (site and/or observation time) - - Parameters - ---------- - site: str, optional - Updated site name - time: str, optional - Updated observation time (yyyy-MM-ddThh:mm:ss) - delta_t: astropy.Quantity, optional - Elapsed time from the original observation time. - e.g., CTA_Info.update(delta_t= -0.5*u.hour) - -> t_new = t_old - 0.5 hour - verbose: optional - """ - - self.frame.update(site=site, time=time, delta_t=delta_t, verbose=verbose) - self.divergent_pointing(self.div, ra=self.pointing["ra"], dec=self.pointing["dec"]) - - def get_pointing_coord(self, icrs=True): - """ - Return pointing coordinates - - Parameters - ---------- - icrs: bool, optional - If True, return (ra, dec). - If False, return (alt, az) - """ - return pointing.pointing_coord(self.table, self.frame, icrs=icrs) - - def set_pointing_coord(self, src=None, ra = None, dec = None, alt=None, az = None, units='deg'): - """ - Set pointing coordinates - - Parameters - ---------- - src: astropy.coordinates.SkyCoord, optional - ra: float, optional - Right ascension of a source - dec: float, optional - Declination of a source - alt: float, optional - Mean altitude of the array - az: float, optional - Mean azumith angle of the array - units: str - Units of RA and DEC; either deg (default) or rad - """ - if type(src) == SkyCoord: - ra = src.icrs.ra.deg - dec = src.icrs.dec.deg - units = 'deg' - - if ra is not None and dec is not None: - src = self.frame.set_source_loc(ra=ra, dec=dec, units=units) - self._pointing["ra"] = src.icrs.ra.value * u.deg - self._pointing["dec"] = src.icrs.dec.value * u.deg - self._pointing["alt"] = src.alt.value * u.deg - self._pointing["az"] = src.az.value * u.deg - elif alt is not None and az is not None: - if units == "deg": - self._pointing["alt"] = alt * u.deg - self._pointing["az"] = az * u.deg - else: - self._pointing["alt"] = alt * u.rad - self._pointing["alt"] = az * u.rad - - for tel in self.telescopes: - tel.__point_to_altaz__(self.pointing["alt"], self.pointing["az"]) - - self.__make_table__() - - def divergent_pointing(self, div, ra=None, dec = None, alt=None, az=None, units="deg"): - """ - Divergent pointing given a parameter div. - Update pointing of all telescopes of the array. - - Parameters - ---------- - div: float between 0 and 1 - ra: float, optioanl - source ra - dec: float, optional - source dec - alt: float, optional - mean alt pointing - az: float, optional - mean az pointing - units: string, optional - either 'deg' (default) or 'rad' - """ - - self.set_pointing_coord(ra=ra, dec = dec, alt=alt, az=az, units=units) - - self._div = div - - if np.abs(div) > 1: #or div < 0: - print("[Error] The div abs value should be lower and 1.") - elif div!=0: - G = pointing.pointG_position(self.barycenter, self.div, self.pointing["alt"], self.pointing["az"]) - for tel in self.telescopes: - alt_tel, az_tel = pointing.tel_div_pointing(tel.position, G) - - if div < 0: - az_tel=az_tel - np.pi - - tel.__point_to_altaz__(alt_tel*u.rad, az_tel*u.rad) - - - self.__make_table__() - - def group_by(self, group = None): - - if type(group) == dict: - groupping = np.zeros(self.size_of_array) - labels = [] - j = 1 - for key in group.keys(): - labels.append(key) - for i in group[key]: - groupping[i-1] = j - j+=1 - tel_group = self._table.group_by(np.asarray(groupping)) - elif group: - tel_group = self._table.group_by("radius") - labels = ["group_{}".format(i+1) for i in range(len(tel_group.groups))] - else: - tel_group = self._table.group_by(np.zeros(self.size_of_array)) - labels = ["_nolegend_"] - return (tel_group, labels) - - - def display(self, projection, ax=None, group=False, **kwargs): - """ - Display the CTA array - - Parameters - ---------- - projection: str - any combination of `x`, `y`, and `z` or `skymap` - ax: pyplot.axes, optional - group: bool or dict, optional - kwargs: args for `pyplot.scatter` - - Returns - ------- - ax: `matplotlib.pyplot.axes` - """ - tel_group, labels = self.group_by(group) - - if projection == 'skymap': - for i, [table, label] in enumerate(zip(tel_group.groups, labels)): - - ax = visual.display_skymap(table, self.frame, - label=labels[i], ax=ax) - return ax - else: - proj = [] - for axis in ["x", "y", "z"]: - if axis in projection: - proj.append(axis) - - if len(proj) == 1: - ax = visual.display_1d(tel_group, proj, ax=ax, labels=labels, **kwargs) - elif len(proj) == 2: - ax = visual.display_2d(tel_group, proj, ax=ax, labels=labels, **kwargs) - else: - ax = visual.display_3d(tel_group, proj, ax=ax, labels=labels, **kwargs) - - return ax - - def skymap_polar(self, group=None, fig=None, filename=None): - """ - Plot skymap - - Parameters - ---------- - group: bool or dict, optional - fig: pyplot.figure, optional - filemane: str, optional - - """ - return visual.skymap_polar(self, group=group, fig=fig, filename=filename) - - def multiplicity_plot(self,fig=None): - """ - Plot multiplicity - - Parameters - ---------- - fig: pyplot.figure, optional - """ - return visual.multiplicity_plot(self, fig=fig) - - def export_cfg(self, filename=None, outdir="./", verbose=False): - """ - Export cfg file. - - Parameters - ---------- - filename: str, optional - A default name is 'CTA-ULTRA6-LaPalma-divX-azX-altX.cfg' - - outdir: str, optional, - - verbose: bool, optional - """ - - if filename==None: - filename = 'CTA-ULTRA6-LaPalma-div{}-az{}-alt{}.cfg'.format( - str(self.div).replace(".", "_"), - str(self.pointing["az"].value).replace(".", "_"), - str(self.pointing["alt"].value).replace(".", "_")) - - with open(outdir+filename, 'w') as f: - f.write('#ifndef TELESCOPE\n') - f.write('# define TELESCOPE 0\n') - f.write('#endif\n') - f.write('#if TELESCOPE == 0\n') - f.write(' TELESCOPE_THETA={:.2f} \n'.format(90 - self.pointing["alt"].value)) - f.write(' TELESCOPE_PHI={:.2f} \n'.format(self.pointing["az"].value)) - f.write('\n% Global and default configuration for things missing in telescope-specific config.\n') - f.write('# include \n') - for n, tel in enumerate(self.table, 1): - zd = 90 - tel['alt'] - f.write('\n#elif TELESCOPE == {:d}\n'.format(n)) - if n <= 4: - f.write('# include \n') - else: - f.write('# include \n') - f.write(' TELESCOPE_THETA={:.2f}\n'.format(zd)) - f.write(' TELESCOPE_PHI={:.2f}\n'.format(360 - tel["az"])) - f.write('#else\n') - f.write(' Error Invalid telescope for CTA-ULTRA6 La Palma configuration.\n') - f.write('#endif\n') - f.write('trigger_telescopes = 2 % In contrast to Prod-3 South we apply loose stereo trigger immediately\n') - f.write('array_trigger = array_trigger_ultra6_diver-test.dat\n') - - if verbose: - with open(outdir+filename, 'r') as f: - for line in f.readlines(): - print(line) - diff --git a/divtel/visualization_old.py b/divtel/visualization_old.py deleted file mode 100644 index 7ca9e11..0000000 --- a/divtel/visualization_old.py +++ /dev/null @@ -1,353 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -from ipywidgets import widgets -from IPython.display import display - -import copy -import astropy.units as u -from astropy.coordinates import SkyCoord -from astroplan.plots import plot_sky -from astroplan import FixedTarget -from . import utils -from .const import COLORS -from . import pointing - -from matplotlib.transforms import Affine2D -from astropy.visualization.wcsaxes import SphericalCircle - -from shapely.geometry import mapping -from descartes import PolygonPatch - -def display_1d(table, proj, ax=None, labels=None, **kwargs): - xb = utils.calc_mean(table, proj[0]) - - ax = plt.figure().add_subplot(111) if ax is None else ax - - for i, [tels, label] in enumerate(zip(table.groups, labels)): - c = COLORS(i) - for val in tels[proj[0]]: - ax.axvline(val, label=label, color=c, **kwargs) - label='_nolegend_' - - ax.axvline(xb, color="r", label='barycenter', **kwargs) - ax.set_xlabel(f"{proj[0]} [m]") - ax.set_yticks([0, 1]) - ax.legend(frameon=False) - - return ax - -def display_2d(table, proj, ax=None, labels=None, **kwargs): - if ax is None: - ax = plt.figure().add_subplot(111) - - scale = 1 - - b_output = utils.calc_mean(table, [proj[0], proj[1], f"p_{proj[0]}", f"p_{proj[1]}"]) - - for i, [tels, label] in enumerate(zip(table.groups, labels)): - xx = tels[proj[0]] - yy = tels[proj[1]] - xv = tels[f"p_{proj[0]}"] - yv = tels[f"p_{proj[1]}"] - ids = tels["id"] - - s = ax.scatter(xx, yy, label=label, **kwargs) - ax.quiver(xx, yy, xv, yv, color=s.get_facecolor()) - - for i, x, y in zip(ids, xx, yy): - ax.annotate(i, (x, y)) - - ax.scatter(b_output[0], b_output[1], marker='+', label='barycenter', color="r") - ax.quiver(*b_output, color="r") - ax.set_xlabel(f"{proj[0]} [m]") - ax.set_ylabel(f"{proj[1]} [m]") - - ax.grid('on') - ax.axis('equal') - xlim = ax.get_xlim() - ylim = ax.get_ylim() - ax.set_xlim(xlim[0] - 0.25 * np.abs(xlim[0]), xlim[1] + 0.25 * np.abs(xlim[1])) - ax.set_ylim(ylim[0] - 0.25 * np.abs(ylim[0]), ylim[1] + 0.25 * np.abs(ylim[1])) - ax.legend(frameon=False) - - return ax - -def display_3d(table, proj, ax=None, labels=None, **kwargs): - ax = plt.figure().add_subplot(111, projection='3d') - - scale = 1 - - max_range = [] - for axis in ["x", "y", "z"]: - max_range.append(table[axis].max() - table[axis].min()) - max_range = max(max_range) - - for i, [tels, label] in enumerate(zip(table.groups, labels)): - xx = tels["x"] - yy = tels["y"] - zz = tels["z"] - c = COLORS(i) - ax.quiver(xx, yy, zz, - tels["p_x"], tels["p_y"], tels["p_z"], - length=max_range, - label=label, - color=c, - ) - - Xb = scale * max_range * np.mgrid[-1:2:2, -1:2:2, -1:2:2][0].flatten() + scale * (xx.max() + xx.min()) - Yb = scale * max_range * np.mgrid[-1:2:2, -1:2:2, -1:2:2][1].flatten() + scale * (yy.max() + yy.min()) - Zb = scale * max_range * np.mgrid[-0.01:2:2, -0.01:2:2, -0.01:2:2][2].flatten() + scale * (zz.max() + zz.min()) - - for xb, yb, zb in zip(Xb, Yb, Zb): - ax.plot([xb], [yb], [zb], 'w') - - xx = utils.calc_mean(table, proj[0]) - yy = utils.calc_mean(table, proj[1]) - zz = utils.calc_mean(table, proj[2]) - xbv = utils.calc_mean(table, f"p_{proj[0]}") - ybv = utils.calc_mean(table, f"p_{proj[1]}") - zbv = utils.calc_mean(table, f"p_{proj[2]}") - - ax.quiver(xx, yy, zz, - xbv, ybv, zbv, - color="r", - length=max_range, - label='barycenter', - ) - - ax.set_xlabel('x [m]') - ax.set_ylabel('y [m]') - ax.set_zlabel('z [m]') - ax.legend(frameon=False) - - return ax - -def display_barycenter(table, proj, ax=None, labels=None, fig=None, **kwargs): - if fig is None: - fig = plt.figure() - - if ax is None: - ax = fig.add_subplot(111) - - scale = 1 - - for i, (tab, label) in enumerate(zip(table.groups, labels)): - output = utils.calc_mean(tab, [proj[0], proj[1], f"p_{proj[0]}", f"p_{proj[1]}"]) - s = ax.scatter(output[0], output[1], color=COLORS(i),) - ax.quiver(*output, color=s.get_facecolor()) - - ax.annotate(label, (output[0], output[1])) - - ax.set_xlabel(f"{proj[0]} [m]") - ax.set_ylabel(f"{proj[1]} [m]") - - ax.grid('on') - ax.axis('equal') - xlim = ax.get_xlim() - ylim = ax.get_ylim() - ax.set_xlim(xlim[0] - 0.25 * np.abs(xlim[0]), xlim[1] + 0.25 * np.abs(xlim[1])) - ax.set_ylim(ylim[0] - 0.25 * np.abs(ylim[0]), ylim[1] + 0.25 * np.abs(ylim[1])) - ax.legend(frameon=False) - - return ax - -def interactive_barycenter(array, proj="xy", overwrite=True, group=False): - if overwrite: - new_array = array - else: - new_array = copy.deepcopy(array) - - fig = plt.figure() - - def update(div=0, az=0, alt=0): - new_array.divergent_pointing(div, az=az, alt=alt, units='deg') - new_array.__convert_units__(toDeg=True) - plt.cla() - grouped_table, labels = new_array.group_by(group) - display_barycenter(grouped_table, proj, labels=labels, fig=fig) - fig.canvas.draw_idle() - - div_s = widgets.FloatLogSlider(value=new_array.div, base=10, min=-4, max=0, step=0.2, description='Divergence') - az_s = widgets.FloatSlider(value=new_array.pointing["az"].value, min=0, max=360, step=0.01, description='Azimuth [deg]') - alt_s = widgets.FloatSlider(value=new_array.pointing["alt"].value, min=0, max=90, step=0.01, description='Altitude [deg]') - - ui = widgets.HBox([div_s, alt_s, az_s]) - out = widgets.interactive_output(update, {'div': div_s, 'az': az_s, 'alt': alt_s}) - display(ui, out) - - return new_array - -def display_skymap(table, frame, ax=None, **kwargs): - ax = plt.figure().add_subplot(111, projection='polar') if ax is None else ax - - radec = pointing.pointing_coord(table, frame, icrs=True) - point = SkyCoord(ra=radec.ra, dec=radec.dec) - target = FixedTarget(coord=point, name="source") - - plot_sky(target, frame.observer, frame.t_obs, ax=ax, style_kwargs=kwargs) - - return ax - -def skymap_polar(array, group=False, fig=None, filename=None): - if fig is None: - fig = plt.figure() - else: - ax = fig.gca() - ax.set_xticklabels([]) - ax.set_yticklabels([]) - - tr = Affine2D().scale(np.pi/180., 1.).translate(+np.pi/2., 0) + plt.polar().get_transform() - - n = 20 - extreme_finder = angle_helper.ExtremeFinderCycle(10, 10, - lon_cycle=360, - lat_cycle=None, - lon_minmax=None, - lat_minmax=(-90, 90), - ) - - grid_locator1 = angle_helper.LocatorDMS(12) - - tick_formatter1 = angle_helper.FormatterDMS() - - grid_helper = GridHelperCurveLinear(tr, - extreme_finder=extreme_finder, - grid_locator1=grid_locator1, - tick_formatter1=tick_formatter1 - ) - - ax1 = plt.subplot(projection=grid_helper) - - array.__convert_units__(toDeg=True) - tel_group, labels = array.group_by(group) - - for tel_table, label in zip(tel_group.groups, labels): - s = ax1.scatter(tel_table["az"], tel_table["alt"], label=label, - s=20, edgecolor="black", transform=ax1.transData, zorder=10) - - for tel in tel_table: - r = SphericalCircle((tel["az"] * u.deg, tel["alt"] * u.deg), tel["radius"] * u.deg, - color=s.get_facecolor()[0], alpha=0.1, transform=ax1.transData) - ax1.add_patch(r) - ax1.annotate(tel["id"], (tel["az"], tel["alt"]), fontsize=12, xytext=(4, 4), - color="black", textcoords='offset pixels', zorder=10) - - ax1.grid(True) - ax1.set_xlabel("Azimuth [deg]", fontsize=20) - ax1.set_ylabel("Altitude [deg]", fontsize=20) - ax1.legend(loc=1) - - if filename is not None: - plt.savefig(filename) - plt.show(block=False) - -def interactive_polar(array, overwrite=True, group=False): - if overwrite: - new_array = array - else: - new_array = copy.deepcopy(array) - - fig = plt.figure() - - def update(div=0, az=0, alt=0): - new_array.divergent_pointing(div, az=az, alt=alt, units='deg') - new_array.__convert_units__(toDeg=True) - plt.cla() - new_array.skymap_polar(group=group, fig=fig) - fig.canvas.draw_idle() - - div_s = widgets.FloatLogSlider(value=new_array.div, base=10, min=-4, max=0, step=0.2, description='Divergence') - az_s = widgets.FloatSlider(value=new_array.pointing["az"].value, min=0, max=360, step=0.01, description='Azimuth [deg]') - alt_s = widgets.FloatSlider(value=new_array.pointing["alt"].value, min=0, max=90, step=0.01, description='Altitude [deg]') - - ui = widgets.HBox([div_s, alt_s, az_s]) - out = widgets.interactive_output(update, {'div': div_s, 'az': az_s, 'alt': alt_s}) - display(ui, out) - - return new_array - - -def multiplicity_plot(array, m_cut=0, fig=None): - - m, overlaps, geoms = array.hFoV(m_cut=m_cut, full_output=True) - max_m = int(array.size_of_array) - ave_multi = np.average(m[:, 0], weights=m[:, 1]) - var_multi = np.average((m[:, 0] - ave_multi) ** 2, weights=m[:, 1]) - - if fig is None: - fig, (ax, ax_mul) = plt.subplots(1, 2, figsize=(10, 4)) - ax_cb = fig.add_axes([0.44, 0.15, 0.01, 0.7]) - - cmap = plt.cm.get_cmap('rainbow') - color_list = cmap(np.linspace(0, 1, max_m)) - bounds = np.arange(max_m + 1) + 1 - - ax_cb = fig.add_axes([0.44, 0.15, 0.01, 0.7]) - - minmax = [] - for i, pol in enumerate(geoms): - colore = int(overlaps[i]) - pol_map = mapping(pol) - ax.add_patch(PolygonPatch(pol_map, color=color_list[colore - 1])) - patch_az = np.asarray(pol_map['coordinates'])[0][:][0] - minmax.append([min(patch_az), max(patch_az)]) - minmax = np.asarray(minmax) - - norm = plt.colors.BoundaryNorm(bounds, cmap.N) - - cb1 = mpl.colorbar.ColorbarBase(ax_cb, - norm=norm, - cmap=cmap, - boundaries=bounds, - orientation='vertical', - label='Multiplicity') - cb1.set_ticks(np.arange(0, max_m + 1, step=4) + 1) - cb1.set_ticklabels(np.arange(0, max_m + 1, step=4)) - - ax.set_xlabel("Azimuth [deg]") - ax.set_ylabel("Altitude [deg]") - ax.set_xlim(np.min(array.table["az"]) - 5, np.max(array.table["az"]) + 5) - ax.set_ylim(np.min(array.table["alt"]) - 5, np.max(array.table["alt"]) + 5) - - ax.text(0.9, 0.9, r"Average: {:.1f} $\pm$ {:.1f}".format(ave_multi, np.sqrt(var_multi)), - ha="right", transform=ax.transAxes) - - ax_mul.bar(m[:, 0], m[:, 1]) - ax_mul.text(0.9, 0.9, "Total hFoV = {:.0f}".format(sum(m[:, 1][m[:, 0] >= m_cut])), - ha="right", transform=ax_mul.transAxes) - ax_mul.set_xticks(np.arange(0, max_m + 1, step=2)) - ax_mul.set_xlim(0.5, max_m + 0.5) - ax_mul.set_ylabel('HFOV') - ax_mul.set_xlabel('Multiplicity') - - plt.show() - -def interactive_multiplicity(array, overwrite=True): - - if overwrite: - new_array = array - else: - new_array = copy.deepcopy(array) - - fig, (ax, ax_mul) = plt.subplots(1, 2, figsize=(10, 4)) - ax_cb = fig.add_axes([0.44, 0.15, 0.01, 0.7]) - - def update(div=0, az=0, alt=0): - - new_array.divergent_pointing(div, az=az, alt=alt, units='deg') - new_array.__convert_units__(toDeg=True) - plt.cla() - new_array.multiplicity_plot(fig=(fig, (ax, ax_mul))) - fig.canvas.draw_idle() - - div_s = widgets.FloatLogSlider(value=new_array.div, base=10, min=-4, max=0, step=0.2, description='Divergence') - az_s = widgets.FloatSlider(value=new_array.pointing["az"].value, min=0, max=360, step=0.01, - description='Azimuth [deg]') - alt_s = widgets.FloatSlider(value=new_array.pointing["alt"].value, min=0, max=90, step=0.01, - description='Altitude [deg]') - - ui = widgets.HBox([div_s, alt_s, az_s]) - out = widgets.interactive_output(update, {'div': div_s, 'az': az_s, 'alt': alt_s}) - display(ui, out) - - return new_array