diff --git a/measureEccentricity/eccDefinition.py b/measureEccentricity/eccDefinition.py index 6ba6012b..457823ce 100644 --- a/measureEccentricity/eccDefinition.py +++ b/measureEccentricity/eccDefinition.py @@ -72,7 +72,7 @@ def __init__(self, dataDict, spline_kwargs=None, extra_kwargs=None): self.t[1] - self.t[0]) if "hlm_zeroecc" in dataDict: - self.compute_res_amp_and_omega() + self.compute_res_amp_and_omega22() # Sanity check various kwargs and set default values self.spline_kwargs = check_kwargs_and_set_defaults( @@ -105,7 +105,8 @@ def get_default_extra_kwargs(self): "num_orbits_to_exclude_before_merger": 1, "extrema_finding_kwargs": {}, # Gets overriden in methods like # eccDefinitionUsingAmplitude - "debug": True + "debug": True, + "omega22_averaging_method": "average_between_extrema" } return default_extra_kwargs @@ -165,46 +166,92 @@ def interp_extrema(self, extrema_type="maxima"): f"Sufficient number of {extrema_type} are not found." " Can not create an interpolator.") - def measure_ecc(self, tref_in): + def measure_ecc(self, tref_in=None, fref_in=None): """Measure eccentricity and mean anomaly at reference time. parameters: ---------- tref_in: - Input reference time at which to measure eccentricity and mean anomaly. - Can be a single float or an array. NOTE: eccentricity/mean_ano are + Input reference time at which to measure eccentricity and mean + anomaly. + Can be a single float or an array. + NOTE: eccentricity/mean_ano are returned on a different time array tref_out, described below. + fref_in: + Input reference frequency at which to measure the eccentricity and + mean anomaly. It can be a single float or an array. + NOTE: eccentricity/mean anomaly are returned on a different freq + array fref_out, described below. + + Given an fref_in, we find the corresponding tref_in such that, + omega22_average(tref_in) = 2 * pi * fref_in. + Here, omega22_average(t) is a monotonically increasing average + frequency that is computed from the instantaneous omega22(t). + Note that this is not a moving average; depending on which averaging + method is used (see the omega22_averaging_method option below), + it means slightly different things. + + Currently, following options are implemented to calculate the + omega22_average + - "average_between_extrema": Mean of the omega22 given by the + spline through the peaks and the spline through the troughs. + - "orbital_average_at_extrema": A spline through the orbital + averaged omega22 evaluated at all available extrema. + - "omega22_zeroecc": omega22 of the zero eccentricity waveform + The default is "average_between_extrema". A method could be passed + through the "extra_kwargs" option with the key + "omega22_averaging_method". + returns: -------- - tref_out: - Output reference time where eccentricity and mean anomaly are - measured. - This is set as tref_out = tref_in[tref_in >= tmin && tref_in < tmax], + tref_out/fref_out: + tref_out is the output reference time, while fref_out is the + output reference frequency, at which eccentricity and mean anomaly + are measured. + + NOTE: Only one of these is returned depending on whether tref_in or + fref_in is provided. If tref_in is provided then tref_out is + returned and if fref_in provided then fref_out is returned. + + tref_out is set as tref_out = tref_in[tref_in >= tmin && tref_in < tmax], where tmax = min(t_peaks[-1], t_troughs[-1]), and tmin = max(t_peaks[0], t_troughs[0]). This is necessary because - eccentricity is computed using interpolants of omega_peaks and - omega_troughs. The above cutoffs ensure that we are not extrapolating - in omega_peaks/omega_troughs. - In addition, if num_orbits_to_exclude_before_merger in extra_kwargs is - not None, only the data up to that many orbits before merger is + eccentricity is computed using interpolants of omega22_peaks and + omega22_troughs. The above cutoffs ensure that we are not + extrapolating in omega22_peaks/omega22_troughs. + In addition, if num_orbits_to_exclude_before_merger in extra_kwargs + is not None, only the data up to that many orbits before merger is included when finding the t_peaks/t_troughs. This helps avoid unphysical features like nonmonotonic eccentricity near the merger. + fref_out is set as fref_out = fref_in[fref_in >= fmin && fref_in < fmax]. + where fmin = omega22_average(tmin)/2/pi, and + fmax = omega22_average(tmax)/2/pi. tmin/tmax are defined above. + ecc_ref: - Measured eccentricity at tref_out. + Measured eccentricity at tref_out/fref_out. mean_ano_ref: - Measured mean anomaly at tref_out. + Measured mean anomaly at tref_out/fref_out. """ - tref_in = np.atleast_1d(tref_in) - omega_peaks_interp, self.peaks_location = self.interp_extrema("maxima") - omega_troughs_interp, self.troughs_location = self.interp_extrema("minima") + self.omega22_peaks_interp, self.peaks_location = self.interp_extrema("maxima") + self.omega22_troughs_interp, self.troughs_location = self.interp_extrema("minima") t_peaks = self.t[self.peaks_location] t_troughs = self.t[self.troughs_location] - t_max = min(t_peaks[-1], t_troughs[-1]) - t_min = max(t_peaks[0], t_troughs[0]) + self.t_max = min(t_peaks[-1], t_troughs[-1]) + self.t_min = max(t_peaks[0], t_troughs[0]) + # check that only one of tref_in or fref_in is provided + if (tref_in is not None) + (fref_in is not None) != 1: + raise KeyError("Exactly one of tref_in and fref_in" + " should be specified.") + elif tref_in is not None: + tref_in = np.atleast_1d(tref_in) + else: + fref_in = np.atleast_1d(fref_in) + # get the tref_in and fref_out from fref_in + tref_in, self.fref_out = self.compute_tref_in_and_fref_out_from_fref_in(fref_in) # We measure eccentricity and mean anomaly from t_min to t_max. # Note that here we do not include the t_max. This is because # the mean anomaly computation requires to looking @@ -212,18 +259,24 @@ def measure_ecc(self, tref_in): # period. # If ref time is t_max, which could be equal to the last peak, then # there is no next peak and that would cause a problem. - self.tref_out = tref_in[np.logical_and(tref_in < t_max, - tref_in >= t_min)] + self.tref_out = tref_in[np.logical_and(tref_in < self.t_max, + tref_in >= self.t_min)] # Sanity checks + # check that fref_out and tref_out are of the same length + if fref_in is not None: + if len(self.fref_out) != len(self.tref_out): + raise Exception(f"Length of fref_out {len(self.fref_out)}" + " is different from " + f"Length of tref_out {len(self.tref_out)}") # Check if tref_out is reasonable if len(self.tref_out) == 0: - if tref_in[-1] > t_max: - raise Exception(f"tref_in is later than t_max={t_max}, " + if tref_in[-1] > self.t_max: + raise Exception(f"tref_in is later than t_max={self.t_max}, " "which corresponds to min(last periastron " "time, last apastron time).") - if tref_in[0] < t_min: - raise Exception(f"tref_in is earlier than t_min={t_min}, " + if tref_in[0] < self.t_min: + raise Exception(f"tref_in is earlier than t_min={self.t_min}, " "which corresponds to max(first periastron " "time, first apastron time).") raise Exception("tref_out is empty. This can happen if the " @@ -244,15 +297,15 @@ def measure_ecc(self, tref_in): if self.tref_out[0] < t_peaks[0] or self.tref_out[-1] >= t_peaks[-1]: raise Exception("Reference time must be within two peaks.") - # compute eccentricty from the value of omega_peaks_interp - # and omega_troughs_interp at tref_out using the fromula in + # compute eccentricty from the value of omega22_peaks_interp + # and omega22_troughs_interp at tref_out using the fromula in # ref. arXiv:2101.11798 eq. 4 - self.omega_peak_at_tref_out = omega_peaks_interp(self.tref_out) - self.omega_trough_at_tref_out = omega_troughs_interp(self.tref_out) - self.ecc_ref = ((np.sqrt(self.omega_peak_at_tref_out) - - np.sqrt(self.omega_trough_at_tref_out)) - / (np.sqrt(self.omega_peak_at_tref_out) - + np.sqrt(self.omega_trough_at_tref_out))) + self.omega22_peak_at_tref_out = self.omega22_peaks_interp(self.tref_out) + self.omega22_trough_at_tref_out = self.omega22_troughs_interp(self.tref_out) + self.ecc_ref = ((np.sqrt(self.omega22_peak_at_tref_out) + - np.sqrt(self.omega22_trough_at_tref_out)) + / (np.sqrt(self.omega22_peak_at_tref_out) + + np.sqrt(self.omega22_trough_at_tref_out))) @np.vectorize def compute_mean_ano(time): @@ -285,8 +338,11 @@ def compute_mean_ano(time): self.mean_ano_ref = self.mean_ano_ref[0] self.ecc_ref = self.ecc_ref[0] self.tref_out = self.tref_out[0] + if fref_in is not None and len(self.fref_out) == 1: + self.fref_out = self.fref_out[0] - return self.tref_out, self.ecc_ref, self.mean_ano_ref + return_array = self.fref_out if fref_in is not None else self.tref_out + return return_array, self.ecc_ref, self.mean_ano_ref def check_extrema_separation(self, extrema_location, extrema_type="extrema", @@ -372,7 +428,7 @@ def check_monotonicity_and_convexity(self, tref_out, ecc_ref, if any(self.d2ecc_dt > 0): warnings.warn("Ecc(t) is concave.") - def compute_res_amp_and_omega(self): + def compute_res_amp_and_omega22(self): """Compute residual amp22 and omega22.""" self.hlm_zeroecc = self.dataDict["hlm_zeroecc"] self.t_zeroecc = self.dataDict["t_zeroecc"] @@ -410,6 +466,195 @@ def compute_res_amp_and_omega(self): self.res_omega22 = (self.omega22 - self.omega22_zeroecc_interp) + def compute_orbital_averaged_omega22_at_extrema(self, t): + """Compute reference frequency by orbital averaging at extrema. + + We compute the orbital average of omega22 at the periastrons + and the apastrons following: + omega22_avg((t[i]+ t[i+1])/2) = int_t[i]^t[i+1] omega22(t)dt + / (t[i+1] - t[i]) + where t[i] is the time of ith extrema. + We do this for peaks and troughs and combine the results + """ + extrema_locations = {"peaks": self.peaks_location, + "troughs": self.troughs_location} + @np.vectorize + def orbital_averaged_omega22_at_extrema(n, extrema_type="peaks"): + """Compute orbital averaged omega22 between n and n+1 extrema.""" + # integrate omega22 between n and n+1 extrema + # We do not need to do the integration here since + # we already have phase22 available to us which is + # nothing but the integration of omega22 over time. + # We want to integrate from nth extrema to n+1 extrema + # which is equivalent to phase difference between + # these two extrema + integ = (self.phase22[extrema_locations[extrema_type][n+1]] + - self.phase22[extrema_locations[extrema_type][n]]) + period = (self.t[extrema_locations[extrema_type][n+1]] + - self.t[extrema_locations[extrema_type][n]]) + return integ / period + # get the mid points between the peaks as avg time for peaks + t_average_peaks = (self.t[self.peaks_location][1:] + + self.t[self.peaks_location][:-1]) / 2 + omega22_average_peaks = orbital_averaged_omega22_at_extrema( + np.arange(len(self.peaks_location) - 1), "peaks") + # get the mid points between the troughs as avg time for troughs + t_average_troughs = (self.t[self.troughs_location][1:] + + self.t[self.troughs_location][:-1]) / 2 + omega22_average_troughs = orbital_averaged_omega22_at_extrema( + np.arange(len(self.troughs_location) - 1), "troughs") + # combine results from avergae at peaks and toughs + t_average = np.append(t_average_troughs, t_average_peaks) + # sort the times + sorted_idx = np.argsort(t_average) + t_average = t_average[sorted_idx] + # check if the average omega22 are monotonically increasing + if any(np.diff(omega22_average_peaks) <= 0): + raise Exception("Omega22 average at peaks are not strictly " + "monotonically increaing") + if any(np.diff(omega22_average_troughs) <= 0): + raise Exception("Omega22 average at troughs are not strictly " + "monotonically increaing") + omega22_average = np.append(omega22_average_troughs, + omega22_average_peaks) + # sort omega22 + omega22_average = omega22_average[sorted_idx] + return InterpolatedUnivariateSpline(t_average, omega22_average)(t) + + def compute_omega22_average_between_extrema(self, t): + """Find omega22 average between extrema". + + Take mean of omega22 using spline through omega22 peaks + and spline through omega22 troughs. + """ + return ((self.omega22_peaks_interp(t) + + self.omega22_troughs_interp(t)) / 2) + + def compute_omega22_zeroecc(self, t): + """Find omega22 from zeroecc data.""" + return InterpolatedUnivariateSpline( + self.t_zeroecc, self.omega22_zeroecc)(t) + + def get_availabe_omega22_averaging_methods(self): + """Return available omega22 averaging methods.""" + available_methods = { + "average_between_extrema": self.compute_omega22_average_between_extrema, + "orbital_average_at_extrema": self.compute_orbital_averaged_omega22_at_extrema, + "omega22_zeroecc": self.compute_omega22_zeroecc + } + return available_methods + + def compute_tref_in_and_fref_out_from_fref_in(self, fref_in): + """Compute tref_in and fref_out from fref_in. + + Using chosen omega22 average method we get the tref_in and fref_out + for the given fref_in. + + When the input is frequencies where eccentricity/mean anomaly is to be + measured, we internally want to map the input frequencies to a tref_in + and then we proceed to calculate the eccentricity and mean anomaly for + these tref_in in the same way as we do when the input array was time + instead of frequencies. + + We first compute omega22_average(t) using the instantaneous omega22(t), + which can be done in different ways as described below. Then, we keep + only the allowed frequencies in fref_in by doing + fref_out = fref_in[fref_in >= omega22_average(tmin) / (2 pi) && + fref_in < omega22_average(tmax) / (2 pi)] + Finally, we find the times where omega22_average(t) = 2 * pi * fref_out, + and set those to tref_in. + + omega22_average(t) could be calculated in the following ways + - Mean of the omega22 given by the spline through the peaks and the + spline through the troughs, we call this "average_between_extrema" + - Orbital average at the extrema, we call this "orbital_average_at_extrema" + - omega22 of the zero eccentricity waveform, called "omega22_zeroecc" + + User can provide a method through the "extra_kwargs" option with the key + "omega22_averaging_method". Default is "average_between_extrema" + + Once we get the reference frequencies, we create a spline to get time + as function of these reference frequencies. This should work if the + refrence frequency is monotonic which it should be. + + Finally we evaluate this spine on the fref_in to get the tref_in. + """ + self.available_averaging_methods = self.get_availabe_omega22_averaging_methods() + method = self.extra_kwargs["omega22_averaging_method"] + if method in self.available_averaging_methods: + # The fref_in array could have frequencies that is outside the range + # of frequencies in omega22 average. Therefore, we want to create + # a separate array of frequencies fref_out which is created by + # taking on those frequencies that falls within the omega22 average + # Then proceed to evaluate the tref_in based on these fref_out + fref_out = self.get_fref_out(fref_in, method) + # get omega22_average by evaluating the omega22_average(t) + # on t, from tmin to tmax + self.t_for_omega22_average = self.t[ + np.logical_and(self.t >= self.t_min, self.t < self.t_max)] + self.omega22_average = self.available_averaging_methods[ + method](self.t_for_omega22_average) + # check if average omega22 is monotonically increasing + if any(np.diff(self.omega22_average) <= 0): + warnings.warn(f"Omega22 average from method {method} is not " + "monotonically increasing.") + t_of_fref_out = InterpolatedUnivariateSpline( + self.omega22_average / (2 * np.pi), + self.t_for_omega22_average) + tref_in = t_of_fref_out(fref_out) + # check if tref_in is monotonically increasing + if any(np.diff(tref_in) <= 0): + warnings.warn(f"tref_in from fref_in using method {method} is" + " not monotonically increasing.") + return tref_in, fref_out + else: + raise KeyError(f"Omega22 averaging method {method} does not exist." + " Must be one of " + f"{list(self.available_averaging_methods.keys())}") + + def get_fref_out(self, fref_in, method): + """Get fref_out from fref_in that falls within the valid average f22 range. + + Parameters: + ---------- + fref_in: + Input 22 mode reference frequency array. + + method: + method for getting average omega22 + + Returns: + ------- + fref_out: + Slice of fref_in that satisfies + fref_in >= omega22_average(t_min) / 2 pi and + fref_in < omega22_average(t_max) / 2 pi + """ + # get min an max value f22_average from omega22_average + self.omega22_average_min = self.available_averaging_methods[ + method](self.t_min) + self.omega22_average_max = self.available_averaging_methods[ + method](self.t_max) + self.f22_average_min = self.omega22_average_min / (2 * np.pi) + self.f22_average_max = self.omega22_average_max / (2 * np.pi) + fref_out = fref_in[ + np.logical_and(fref_in >= self.f22_average_min, + fref_in < self.f22_average_max)] + if len(fref_out) == 0: + if fref_in[0] < self.f22_average_min: + raise Exception("fref_in is earlier than minimum available " + "frequency " + f"{self.f22_average_min}") + if fref_in[-1] > self.f22_average_max: + raise Exception("fref_in is later than maximum available " + "frequency " + f"{self.f22_average_max}") + else: + raise Exception("fref_out is empty. This can happen if the " + "waveform has insufficient identifiable " + "periastrons/apastrons.") + return fref_out + def make_diagnostic_plots(self, usetex=True, **kwargs): """Make dignostic plots for the eccDefinition method. @@ -453,10 +698,10 @@ def make_diagnostic_plots(self, usetex=True, **kwargs): self.plot_measured_ecc(fig, ax[0]) self.plot_decc_dt(fig, ax[1]) self.plot_mean_ano(fig, ax[2]) - self.plot_extrema_in_omega(fig, ax[3]) + self.plot_extrema_in_omega22(fig, ax[3]) self.plot_phase_diff_ratio_between_peaks(fig, ax[4]) if "hlm_zeroecc" in self.dataDict: - self.plot_residual_omega(fig, ax[5]) + self.plot_residual_omega22(fig, ax[5]) self.plot_residual_amp(fig, ax[6]) fig.tight_layout() return fig, ax @@ -518,7 +763,7 @@ def plot_mean_ano(self, fig=None, ax=None, **kwargs): else: return axNew - def plot_extrema_in_omega(self, fig=None, ax=None, **kwargs): + def plot_extrema_in_omega22(self, fig=None, ax=None, **kwargs): """Plot omega22, the locations of the apastrons and periastrons, and their corresponding interpolants. This would show if the method is missing any peaks/troughs or @@ -528,10 +773,10 @@ def plot_extrema_in_omega(self, fig=None, ax=None, **kwargs): figNew, axNew = plt.subplots() else: axNew = ax - axNew.plot(self.tref_out, self.omega_peak_at_tref_out, + axNew.plot(self.tref_out, self.omega22_peak_at_tref_out, c=colorsDict["periastron"], label=r"$\omega_{p}(t)$", **kwargs) - axNew.plot(self.tref_out, self.omega_trough_at_tref_out, + axNew.plot(self.tref_out, self.omega22_trough_at_tref_out, c=colorsDict["apastron"], label=r"$\omega_{a}(t)$", **kwargs) # plot only upto merger to make the plot readable @@ -585,10 +830,10 @@ def plot_phase_diff_ratio_between_peaks(self, fig=None, ax=None, **kwargs): else: return axNew - def plot_residual_omega(self, fig=None, ax=None, **kwargs): + def plot_residual_omega22(self, fig=None, ax=None, **kwargs): """Plot residual omega22, the locations of the apastrons and periastrons, and their corresponding interpolants. - Useful to look for bad omega data near merger. + Useful to look for bad omega22 data near merger. We also throw away post merger before since it makes the plot unreadble. """ diff --git a/measureEccentricity/measureEccentricity.py b/measureEccentricity/measureEccentricity.py index 0ec7a985..355ce632 100644 --- a/measureEccentricity/measureEccentricity.py +++ b/measureEccentricity/measureEccentricity.py @@ -47,7 +47,8 @@ def get_available_methods(): return models -def measure_eccentricity(tref_in, dataDict, method="Amplitude", +def measure_eccentricity(tref_in=None, fref_in=None, + dataDict=None, method="Amplitude", return_ecc_method=False, spline_kwargs=None, extra_kwargs=None): @@ -60,6 +61,31 @@ def measure_eccentricity(tref_in, dataDict, method="Amplitude", Can be a single float or an array. NOTE: eccentricity/mean_ano are returned on a different time array tref_out, described below. + fref_in: + Input reference frequency at which to measure the eccentricity and + mean anomaly. It can be a single float or an array. + NOTE: eccentricity/mean anomaly are returned on a different freq + array fref_out, described below. + + Given an fref_in, we find the corresponding tref_in such that, + omega22_average(tref_in) = 2 * pi * fref_in. + Here, omega22_average(t) is a monotonically increasing average + frequency that is computed from the instantaneous omega22(t). + Note that this is not a moving average; depending on which averaging + method is used (see the omega22_averaging_method option below), + it means slightly different things. + + Currently, following options are implemented to calculate the + omega22_average + - "average_between_extrema": Mean of the omega22 given by the + spline through the peaks and the spline through the troughs. + - "orbital_average_at_extrema": A spline through the orbital + averaged omega22 evaluated at all available extrema. + - "omega22_zeroecc": omega22 of the zero eccentricity waveform + The default is "average_between_extrema". A method could be passed + through the "extra_kwargs" option with the key + "omega22_averaging_method". + dataDict: Dictionary containing waveform modes dict, time etc. Should follow the format: @@ -97,23 +123,36 @@ def measure_eccentricity(tref_in, dataDict, method="Amplitude", debug: Run additional sanity checks if debug is True. Default: True. + omega22_averaging_method: + Methods for getting average omega22. Default is + "average_between_extrema". For more see fref_in. returns: -------- - tref_out: - Output reference time where eccentricity and mean anomaly are - measured. - This is set as tref_out = tref_in[tref_in >= tmin && tref_in <= tmax], + tref_out/fref_out: + tref_out is the output reference time where eccentricity and mean + anomaly are measured and fref_out is the output reference frequency + where eccentricity and mean anomaly are measured. + + NOTE: Only of these is returned depending on whether tref_in or + fref_in is provided. If tref_in is provided then tref_out is returned + and if fref_in provided then fref_out is returned. + + tref_out is set as tref_out = tref_in[tref_in >= tmin && tref_in < tmax], where tmax = min(t_peaks[-1], t_troughs[-1]), and tmin = max(t_peaks[0], t_troughs[0]). This is necessary because - eccentricity is computed using interpolants of omega_peaks and - omega_troughs. The above cutoffs ensure that we are not extrapolating - in omega_peaks/omega_troughs. - In addition, if num_orbits_to_exclude_before_merger in extra_kwargs is - not None, only the data up to that many orbits before merger is + eccentricity is computed using interpolants of omega22_peaks and + omega22_troughs. The above cutoffs ensure that we are not + extrapolating in omega22_peaks/omega22_troughs. + In addition, if num_orbits_to_exclude_before_merger in extra_kwargs + is not None, only the data up to that many orbits before merger is included when finding the t_peaks/t_troughs. This helps avoid unphysical features like nonmonotonic eccentricity near the merger. + fref_out is set as fref_out = fref_in[fref_in >= fmin && fref_in < fmax]. + where fmin is the frequency at tmin, and fmax is the frequency at tmax. + tmin/tmax are defined above. + ecc_ref: Measured eccentricity at t_ref. Same type as t_ref. @@ -131,7 +170,8 @@ def measure_eccentricity(tref_in, dataDict, method="Amplitude", spline_kwargs=spline_kwargs, extra_kwargs=extra_kwargs) - tref_out, ecc_ref, mean_ano_ref = ecc_method.measure_ecc(tref_in) + tref_out, ecc_ref, mean_ano_ref = ecc_method.measure_ecc( + tref_in=tref_in, fref_in=fref_in) if not return_ecc_method: return tref_out, ecc_ref, mean_ano_ref else: diff --git a/measureEccentricity/utils.py b/measureEccentricity/utils.py index b63fe6e5..4d4b19fc 100644 --- a/measureEccentricity/utils.py +++ b/measureEccentricity/utils.py @@ -53,7 +53,7 @@ def check_kwargs_and_set_defaults(user_kwargs=None, if kw not in default_kwargs: raise ValueError(f"Invalid key {kw} in {name}." " Should be one of " - f"{default_kwargs.keys()}") + f"{list(default_kwargs.keys())}") for kw in default_kwargs.keys(): if kw not in user_kwargs: diff --git a/notebook/measure_eccentricity_at_ref_frequencies.ipynb b/notebook/measure_eccentricity_at_ref_frequencies.ipynb new file mode 100644 index 00000000..ec5309e9 --- /dev/null +++ b/notebook/measure_eccentricity_at_ref_frequencies.ipynb @@ -0,0 +1,271 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "c48881b2-a62c-4af3-b465-1b31357bab2c", + "metadata": {}, + "source": [ + "### Measuring eccentricity at reference frequencies" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "064133b3-5baa-44a5-baef-5f2353e565fd", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home1/md.shaikh/miniconda3/envs/eccimrct/lib/python3.8/site-packages/gwtools/rotations.py:63: UserWarning: Could not import GWFrames, needed for rotations module\n", + " _warnings.warn(\"Could not import GWFrames, needed for rotations module\")\n", + "/home1/md.shaikh/miniconda3/envs/eccimrct/lib/python3.8/site-packages/gwtools/__init__.py:11: UserWarning: Could not import rotations, decompositions, or fitfuncs. These are not needed by GWSurrogate.\n", + " _warnings.warn(\"Could not import rotations, decompositions, or fitfuncs. These are not needed by GWSurrogate.\")\n" + ] + } + ], + "source": [ + "import sys\n", + "sys.path.append(\"../\")\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "from measureEccentricity.measureEccentricity import get_available_methods\n", + "from measureEccentricity import measure_eccentricity\n", + "from measureEccentricity.load_data import load_waveform\n", + "\n", + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "markdown", + "id": "48f5cd1b-31ba-4a79-b14b-819b2c4e5c2c", + "metadata": {}, + "source": [ + "### Load a PN waveform" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "3c393885-f7ef-404d-ae18-af59edc0578d", + "metadata": {}, + "outputs": [], + "source": [ + "lal_kwargs = {\"approximant\": \"EccentricTD\",\n", + " \"q\": 1.0,\n", + " \"chi1\": [0.0, 0.0, 0.0],\n", + " \"chi2\": [0.0, 0.0, 0.0],\n", + " \"Momega0\": 0.01,\n", + " \"ecc\": 0.1,\n", + " \"mean_ano\": 0,\n", + " \"include_zero_ecc\": True}\n", + "dataDict = load_waveform(**lal_kwargs)" + ] + }, + { + "cell_type": "markdown", + "id": "d2009af3-dd57-4f28-af6e-26685d83a095", + "metadata": { + "tags": [] + }, + "source": [ + "### measure eccentricity at a single fref " + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "81e1233b-be3a-4e97-a7f4-eac547e69185", + "metadata": {}, + "outputs": [], + "source": [ + "from measureEccentricity import eccDefinition" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "f70c7ee7-4bcd-40e3-9e07-87e649534bca", + "metadata": {}, + "outputs": [], + "source": [ + "eccDef = eccDefinition.eccDefinition(dataDict)\n", + "omega_averaging_methods = list(eccDef.get_availabe_omega22_averaging_methods().keys())" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "7bbd2daa-aa15-4ab0-b263-e616828fabba", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['average_between_extrema', 'orbital_average_at_extrema', 'omega22_zeroecc']" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "omega_averaging_methods" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "0ec98e0a-3446-45f6-8aeb-636b8286abfd", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Averaging method = average_between_extrema, f_ref = 0.00477464829275686, ecc = 0.08038266012258283, mean anomaly = 6.0765170522356735\n", + "Averaging method = orbital_average_at_extrema, f_ref = 0.00477464829275686, ecc = 0.07918386927241171, mean anomaly = 2.832303771223433\n", + "Averaging method = omega22_zeroecc, f_ref = 0.00477464829275686, ecc = 0.0795874907107213, mean anomaly = 1.817802067299933\n" + ] + } + ], + "source": [ + "fref_in = 0.03 / (2 * np.pi)\n", + "for averaging_method in omega_averaging_methods:\n", + " fref_out, ecc, mean_ano, eccMethod = measure_eccentricity(\n", + " fref_in=fref_in,\n", + " dataDict=dataDict,\n", + " method=\"ResidualAmplitude\", \n", + " return_ecc_method=True,\n", + " extra_kwargs={\"debug\": False, \"omega22_averaging_method\": averaging_method})\n", + " print(f\"Averaging method = {averaging_method}, f_ref = {fref_out}, ecc = {ecc}, mean anomaly = {mean_ano}\")" + ] + }, + { + "cell_type": "markdown", + "id": "c59453d7-21d8-4217-a212-f15b58b20b14", + "metadata": {}, + "source": [ + "### measure eccentricity at a reference frequency array" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "98b5575b-6a83-4034-bcd1-67a29f6b4419", + "metadata": {}, + "outputs": [], + "source": [ + "fref_in = np.arange(0.01, 0.1, 0.001) / (2 * np.pi)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "ecbb02af-ac51-4d7f-995d-a1e795a01557", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtsAAAF3CAYAAAB0RoegAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAB1nElEQVR4nO3dd3hURdvH8e9sekIIvZfQSyC0hA4G6b0oKoIdERXrY33tnUd9sKIIiihKEwVRUSkSmoAh9N5L6ISWQkib94/EGJDoRtlsyu9zXXtlz5k5u/fcWcKdkzlzjLUWERERERG58hzuDkBEREREpLBSsS0iIiIi4iIqtkVEREREXETFtoiIiIiIi6jYFhERERFxERXbIiIiIiIu4unuAFypTJkyNjg4OGs7ISGBgIAA9wVUgChXzlOunKdcOU+5cp5y5TzlynnKlfOUK4iOjj5prS17ubZCXWwHBwezevXqrO3IyEgiIiLcF1ABolw5T7lynnLlPOXKecqV85Qr5ylXzlOuwBizP6c2TSMREREREXERFdsiIiIiIi6iYltERERExEUK9ZxtERERyZ9SUlKIiYkhKSnJ3aFcVlBQEFu3bnV3GAVCUcqVr68vVapUwcvLy+ljVGyLiIhInouJiSEwMJDg4GCMMe4O50/i4uIIDAx0dxgFQlHJlbWW2NhYYmJiqFGjhtPHaRqJiIiI5LmkpCRKly6dLwttkcsxxlC6dOlc/zVGxbaIiIi4hQptKWj+yWc2T4ttY0wPY8x2Y8wuY8wTl2mvb4xZYYy5YIx55JK2icaY48aYTXkXsYiIiIjIP5dnxbYxxgMYC/QEGgJDjDENL+l2CrgfePMyLzEJ6OHKGEVERETyq8jISPr06ZOrY95++20SExNdFJFrrVu3jrlz57o7jH8tL89stwR2WWv3WGuTgWlA/+wdrLXHrbVRQMqlB1trl5BRjIuIiIjkO2lpae4O4U8Ka7Gdmpqax9H8c3m5Gkll4GC27Rig1ZV+E2PMCGAEQPny5YmMjMxqi4+Pv2hbcqZcOU+5cp5y5TzlynnKlfPyU66CgoKIi4tzawxDhgzh0KFDJCUlcffdd5OWlsb+/ft56aWXSEtLY9y4caxdu5Y333yTadOmMW7cOFJSUggLC2PMmDF4eHhQsWJF7r33XhYuXMirr77K4sWL+fHHH0lKSqJVq1a88847GGOIjo5m1KhR+Pv706ZNG+bPn8+qVatIS0vjueeeY+nSpSQnJ3PnnXdy++23XzbexMRETp8+Td++fdm5cyft2rVjzJgxOByOrPdPTk6mRo0afPDBB0yePJnDhw9z1VVXUbp0ae644w6ioqJ47bXX+OCDDxg3bhwbNmxgz549jBw5knnz5rF27Vr+7//+j4SEBEqVKsW4ceOoUKECe/bs4T//+Q+xsbH4+fnx3nvvUbduXUaOHEmxYsVYt24dx48f58UXX2TAgAE55vydd97hm2++ITk5mT59+vDUU0/x3XffMWHCBL799luOHTtGz549mTNnDs888wznz59nyZIlPPzww2zfvp2jR4+yf/9+Spcuzeuvv86DDz7IwYMZ5eV///tfWrduzauvvsr+/fs5evQou3fv5tVXXyUqKor58+dTsWJFZsyYgZeXF6NHj77s9+rvJCUl5e7fkbU2Tx7AYODjbNs3Ae/l0Pd54JHL7A8GNjn7ni1atLDZLVq0yOaF9PR0u+bYmjx5L1fJq1wVBsqV85Qr5ylXzlOunJefcrVly5aLtqs//v0Vf/yd2NhYa621iYmJNiQkxB49etTWqlXLWmvtuXPnbI8ePezSpUvtli1bbJ8+fWxycrK11tq7777bfvbZZ9ZaawE7ffr0P72mtdYOGzbMzpkzx1prbUhIiF2+fLm11trHH3/choSEWGut/eijj+xLL71krbU2KSnJtmjRwu7Zs+ey8S5atMj6+PjY3bt329TUVNulSxf71Vdf2RMnTtgOHTrY+Ph4a621o0ePti+88EJGXqtXtydOnLDWWnvkyBEbFhZmrbX2mmuusWFhYTYmJsZOmjTJPvHEEzY5Odm2adPGHj9+3Fpr7bRp0+xtt91mrbX26quvtjt27LDWWrty5UrbqVMna621t9xyix0wYIBNS0uzmzdvzsrf5fz888/2zjvvtOnp6TYtLc327t3bLl682Fpr7dChQ+17771ne/fubadMmWKttfbTTz+19957b9bxzz33nG3evLlNTEy01lo7ZMgQu3TpUmuttfv377f169fP6teuXTubnJxs161bZ/38/OzcuXOttdYOGDDAzpo16y+/V3/n0s+utdYCq20O9WhentmOAapm264CHM7D988Tp5JO8fLKl5m/fz5vRbxFl+pd3B2SiIhIvrdvdO88f893332XWbNmAXDw4EH27t1LzZo1WblyJRUqVGD79u20a9eOsWPHEh0dTXh4OADnz5+nXLlyAHh4eHDNNddkveaiRYt4/fXXSUxM5NSpU4SEhNChQwfi4uJo27YtADfeeCPff/89APPmzWPDhg3MnDkTgLNnz7Jz584c13Fu2bIlNWvWBDLOzC9btgxfX1+2bNlCu3btAEhOTqZNmzZ/OrZChQrEx8cTFxfHwYMHufHGG1myZAlLly5l0KBBbN++nU2bNtG1a1cgY1pMxYoViY+P59dff2Xw4MFZr3XhwoWs571798bhcNCwYUOOHTuWY77nzZvHvHnzaNasGZDxl5adO3fSsWNH3nvvPRo1akTr1q0ZMmRIjq/Rr18//Pz8AFiwYAFbtmzJajt37lzWX0t69uyJl5cXjRs3Ji0tjR49Mi77a9y4Mfv27QMu/73q27dvju/9T+VlsR0F1DHG1AAOATcAN+bh+7vcyiMreXzJ45xKypha/uKKF2larill/Mq4OTIRERHJLjIykgULFrBixQr8/f2JiIggKSmJ66+/nhkzZhAcHMzAgQMxxmCt5ZZbbuG111770+v4+vri4eEBZEwvuOeee1i9ejVVq1bl+eefJykp6fe/zl+WtZb33nuP7t27OxX3pdMcfo+va9euTJ069W+Pb9OmDZ9++in16tWjQ4cOTJw4kRUrVvC///2PAwcOEBISwooVKy465ty5c5QoUYJ169Zd9jV9fHwuGk9OrLU8+eST3HXXXX9qO3ToEA6Hg2PHjpGeno7DcfnLCgMCArKep6ens2LFiqzi+3IxORwOvLy8svLmcDhITU3N8XvlCnl2gaS1NhUYBfwMbAVmWGs3G2NGGmNGAhhjKhhjYoCHgaeNMTHGmOKZbVOBFUC9zP135FXszirpU5Jzyeeytk9fOM3o30a7MSIRERG5nLNnz1KyZEn8/f3Ztm0bK1euBGDQoEHMnj2bmTNncv311wPQuXNnZs6cyfHjxwE4deoU+/fv/9Nr/l6slSlThvj4+Kyz1SVLliQwMDDrPaZNm5Z1TPfu3fnwww9JSclYG2LHjh0kJCTkGPdvv/3G3r17SU9PZ/r06bRv357WrVuzfPlydu3aBWTM7d6xYwcAgYGBF82N79ixI2+++SYdO3akWbNmLFq0CB8fH4KCgqhXrx4nTpzIKrZTUlLYvHkzxYsXp0aNGnz11VdARtG8fv36XOX797FOnDiR+Ph4IKPAPn78OKmpqdx2221MmTKFBg0aMGbMmMvGfqlu3brx/vvvZ23n9MvA5eT0vXKFPL1du7V2LjD3kn3jsj0/Ssb0kssdm/PfFPKJeqXqcU+Te3h37bsAhJQO4Z4m97g5KhEREblUjx49GDduHKGhodSrV4/WrVsDGYVxw4YN2bx5My1btgSgYcOGvPzyy3Tr1o309HS8vLwYO3Ys1atXv+g1S5QowZ133knjxo0JDg7OmnYC8Mknn3DnnXcSEBBAREQEQUFBAAwfPpx9+/bRvHlzrLWULVuW2bNn5xh3mzZteOKJJ9i4cSMdO3Zk4MCBOBwOJk2axJAhQ7Kmd7z88svUrVuXESNG0LNnTypWrMiiRYvo0KEDBw8epGPHjnh4eFC1alXq168PgLe3NzNnzuT+++/n7NmzpKam8uCDDxISEsKXX37J3Xffzcsvv0xKSgo33HADTZo0yVXOu3XrxtatW7OmuBQrVowvvviCcePG0aFDBzp06EDTpk0JDw+nd+/edOrUidGjR9O0aVOefPLJP73eu+++y7333ktoaCipqal07NiRcePG/anf5fzV9+pKM391ur+gCwsLs6tXr87ajoyMJCIiwqXvmZqeyh0/30HrSq0Z3ng4Xg4vl76fq+RFrgoL5cp5ypXzlCvnKVfOy0+52rp1Kw0aNHB3GDmKi4sjMDDwir1efHw8xYoVA2D06NEcOXKEd95554q9vjtd6Vzld5f77Bpjoq21YZfrn6dntosCT4cnE2tch0fNTlBAC20RERG5sn744Qdee+01UlNTqV69OpMmTXJ3SJJHVGxfaYtew2PxaAgfDr3/5+5oREREJB+4/vrrs+aA/52NGzdy0003XbTPx8eHVatWuSK0K66gx3+lqdi+0ur3gqX/g6iPoXZXqKc7zIuIiIjzGjdunKuL/fKbgh7/lZaXt2svGio2gS7PZTz/9h6I+2O9yfjkeCZsmEBaev67nauIiIiIXHkqtl2h9b1QMwISY2H23ZCezqojqxg0ZxDvrn2XSZsnuTtCEREREckDKrZdweGAAePArxTsXsis+Q8yfN5wjiQcAWDsurFsP7U9q3thXhFGREREpChTse0qxStC/4yF1jvsXkkJnxJZTSnpKUzYOAGAsYt28fHSve6IUERERERcTMW2K9XvDQM/oswdv/B066cBMBhuaXgLr7R/BYB+TSrx4eLdbDp01p2RioiIyGX8vjb2pZ599lkWLFgAwNtvv01iYuLfvlZERATZ7/9R2L366qtX9PX27dvHlClTruhr5gUV267W5AbwKUb34O7cFnIbE7tP5JHwR/Dx8AGgail/nu3TkAemreV8si6cFBERyQ+staSnp+fY/uKLL9KlSxfA+WLbHdLS3Fdb5GWxnZqaekXf60pSsZ1XkhN4+MgBws7F/qlpQLPKNK4cxEs/bHFDYCIiIvnA80E5P1Z/+ke/1Z/+dd9cGDNmDI0aNaJRo0a8/fbb7Nu3jwYNGnDPPfdk3dYc4D//+Q/Nmzenc+fOnDhxAoBbb72VmTNn8u6773L48GE6depEp06dALj77rsJCwsjJCSE5557zul4Lnfcjz/+yHXXXZfVJzIykr59+wIwb9482rRpQ/PmzRk8eDDx8fEABAcH8+KLL9K+fXu++uorJkyYQHh4OE2aNOGaa67J+sVg9+7dtG7dmvDwcJ599tmLzuK/8cYbhIeHExoa+rdjGDJkCC1atCAkJITx48cD8MQTT3D+/HmaNm3K0KFDczz2iy++oGXLljRt2pS77rqLtLQ0oqKiCA0NJSkpiYSEBEJCQti0aRNPPPEES5cupWnTprz11ltMmjSJwYMH07dvX7p160ZCQgK333474eHhNGvWjG+//RaASZMmMWDAAPr27UuNGjV4//33GTNmDM2aNaN169acOnUKIMc8/VsqtvPKhukQ/SnMvng5wN+9OKARS3ee4OfNR90QnIiISNESHR3Np59+yqpVq1i5ciUTJkzg9OnTbN++nZtvvplly5ZRvXp1EhISaN68OWvWrOGqq67ihRdeuOh17r//fipVqsSiRYtYtGgRAK+88gqrV69mw4YNLF68mA0bNjgV0+WO69q1KytXriQhIQGA6dOnc/3113Py5ElefvllFixYwJo1awgLC2PMmDFZr+Xr68uyZcu44YYbGDRoEFFRUaxfv54GDRrwySefAPDAAw/wwAMPEBUVRaVKlbKOnTdvHjt37uS3335j3bp1REdHs2TJkhzjHjt2LNHR0axevZp3332X2NhYRo8ejZ+fH+vWrePLL7+87HFbt25l+vTpLF++nHXr1uHh4cGXX35JeHg4/fr14+mnn+axxx5j2LBhNGrUiNGjR9OhQwfWrVvHQw89BMCKFSv47LPP+OWXX3jllVe4+uqriYqKYtGiRTz66KNZedu0aRNTpkzht99+46mnnsLf35+1a9fSpk0bPv/8c4Ac8/Rv6aY2eaX5rbB5FuxdkrH+9tCZYExWc6CPJ29f34y7Jq+mSZUSVAjydV+sIiIiee15J69dCrst4/EvLVu2jIEDBxIQEABkFFpLly6levXqtG7dmri4OAAcDkfWnR+HDRvGoEGD/va1Z8yYwfjx40lNTeXIkSNs2bKF0NDQf3xcjx49+O6777j22mv54YcfeP3111m8eDFbtmyhXbt2ACQnJ9OmTZus18p+t8pNmzbx9NNPc+bMGeLj4+nevTuQUajOnj0bgBtvvJFHHnkEyCi2582bR7NmzQCIj49n586ddOzY8bJxjxs3jrlz5wJw8OBBdu7cSenSpf92vAsXLiQ6Oprw8HAAzp8/T7ly5YCMOfHh4eH4+vry7rvv5vgaXbt2pVSpUllxz5kzhzfffBOApKQkDhw4AECnTp0IDAwkMDCQoKCgrL8ONG7cOOuXoZzy9G+p2M4rDgcM/Ag+bAu7FsCqj6D1SABWH13NW9Fv8X7n97mpdTD/+Wodk29vhcNh/uZFRURE5J/Iadnd34vvnBjz1/837927lzfffJOoqChKlizJrbfeSlJS0t/G81fHXX/99YwdO5ZSpUoRHh5OYGAg1lq6du3K1KlT/3Yct956K7Nnz6ZJkyZMmjSJyMjIv4zFWsuTTz7JXXfd9bdxR0ZGEhkZyYoVK/D39yciIsKp8f7+Prfccguvvfban9pOnTpFfHw8KSkpJCUl5fh9yb7fWsvXX39NvXr1LuqzatUqfHx8srYdDkfWtsPhyJrvnds8OUvTSPJS8UrQ772M5/OfJe3oRj5a/xF3zLuDDSc38PTyp7k7ogYXUtKZsHSPe2MVEREpxDp27Mjs2bNJTEwkISGBWbNm0aFDhz/1S09PZ+bMmQBMmTKF9u3b/6lPYGBg1pnwc+fOERAQQFBQEMeOHePHH390Kp6/Oi4iIoI1a9YwYcKErDPWrVu3Zvny5ezatQuAxMREduzYcdnXjouLo2LFiqSkpFw0paN169Z8/fXXAEybNi1rf/fu3Zk4cWLWHPBDhw5x/Pjxy7722bNnKVGiBP7+/mzbto2VK1dmtXl5eZGSkpLjmDt37szMmTOzXvvUqVPs378fgBEjRvDSSy8xdOhQHn/8ceDiPF9O9+7dee+997J+kVq7dm2OfS8npzz9Wzqzndca9IXmt8Caz/hmzm2873Mhq2lJzBKmbPuCt66/jgFjl9OudhkaVc7dxR4iIiLy95o3b86tt95Ky5YtARg+fDglS5b8U7+AgAA2b95MixYtCAoKYvr06X/qM2LECHr27EnFihVZtGgRzZo1IyQkhJo1a2ZN8/g7TZo0yfE4Dw8P+vTpw6RJk/jss88AKFu2LJMmTWLIkCFcuJBRS7z88svUrVv3T6/90ksv0apVK6pXr07jxo2zCta3336bYcOG8b///Y/evXsTFJRRc3Tr1o2tW7dmTUspVqwYX3zxRdYUj+x69OjB+++/T2hoKPXq1aN169YX5SU0NJTmzZtftnht2LAhL7/8Mt26dSM9PR0vLy/Gjh3L4sWL8fT05MYbbyQtLY22bdvyyy+/0KFDBzw9PWnSpAm33nrrn75fzzzzDA8++CChoaFYawkODub77793Kv9/lad/yxTmuxeGhYXZ7OtZRkZGEhER4b6AfpecAB91JKViU271TWBD7OasprtC72JUs1F8u+4Q7yzYyff3t8ffO+9/J8o3uSoAlCvnKVfOU66cp1w5Lz/lauvWrTRo0MDdYeQoLi6OwMBAd4fhUomJifj5+WGMYdq0aUydOjVrBY/cKAq5yu5yn11jTLS1Nuxy/XVm2x28A2D4Arz8SvJG/GEGfzcYD+PBax1eo13ljN9k+zetzOLtJ3jp+y28NujvL6oQERERyY3o6GhGjRqFtZYSJUowceJEd4dUKKnYdhe/jD99VCpWiXc6/JdqPqUoV7bhRV1e6B9C73eX8ePGI/RsXNEdUYqIiMgV1qpVq6ypH7+bPHkyjRs3ztM4OnTowPr1653qGxsbS+fOnf+0f+HChXh7e//jY51ZtaSgU7Htbqf2EvbtI+ATCLf+AJ5/fGADfb14+4amjPh8NU2rlaBikJ8bAxUREZErYdWqVe4OIddKly7NunXrLtv2d3Ob/+rYokCrkbibTyBcOAcxv8H8Z/7U3LxaSW5uE8zD09eTll5459eLiIiIFEYqtt0toAxc9zk4vGDVONg486LmtcfXcnO78qSmpzN+iZYDFBERESlIVGznB1XCoEfmgu5z7ofj20i36YzfMJ5bf7qVF1c+z5jrmjBh6R42xJxxa6giIiIi4jwV2/lF+HBofB2kJHB2xjBG/jyc99a+R7pNZ/7++fx6/Hue7xfCA9PWkXAh1d3RioiIiIgTVGznF8ZA37ehbAN8Y3dy+uy+i5rfWP0Gbet606J6SZ6atTHH28yKiIhIwTF06FDq1atHo0aNuP3227PuuPjll18SGhpKaGgobdu2dXrVEMl/VGznJ94BcP0X+PR5hzd7TMTf0x+Akj4lebvT25TxK8NL/Rux7Wgck1fud3OwIiIiV07jzxpf9MjJVzu+uqjf878+n3dBusDQoUPZtm0bGzdu5Pz583z88ccA1KhRg8WLF7NhwwaeeeYZRowY4ZL3T0tLc8nryh9UbOc3ZWpDi1upHhTM822fp0X5FnzV9yvaV24PgJ+3B+OGteCdBTuJ3n/azcGKiIgUXGPGjKFRo0Y0atSIt99+m3379lG/fn2GDx9Oq1atGDp0KAsWLKBdu3bUqVOH3377DYCEhARuv/12wsPDadasWdZdFxMTE7nuuusIDQ3l+uuvp1WrVvx+J+u7776bsLAwQkJCeO6557Ji6NWrF8YYjDG0bNmSmJgYANq2bZt1O/LWrVtn7b+ccePG0bRpU5o2bUqNGjXo1KkTAPPmzaNNmzY0b96cwYMHEx8fD0BwcDAvvvgi7du356uvvmLq1Kk0btyYRo0a8fjjj2e9bk7HR0VF0bZtW5o0aULLli2Ji4sjLS2NRx55hMaNGxMaGsp77713Rb5HhYGK7XysZ7FaTDwWS/n0i/cHlwlg9DWhjJqyhhNxFy5/sIiIiOQoOjqaTz/9lFWrVrFy5UomTJjA6dOn2bVrFw888AArVqxg27ZtTJkyhWXLlvHmm2/y6quvAvDKK69w9dVXExUVxaJFi3j00UdJSEjggw8+oGTJkllno6Ojo7Pe75VXXmH16tVs2LAh64x1dikpKUyePJkePXr8KdZPPvmEnj175jiWkSNHsm7dOqKioqhSpQoPP/wwJ0+e5OWXX2bBggWsWbOGsLAwxowZk3WMr68vy5Yto2PHjjz++OP88ssvWa8xe/bsHI9PTk7m+uuv55133mH9+vUsWLAAPz8/xo8fz969e1m7di0bNmxg6NCh//ZbVGjopjb52byncexdAjNvh5u/BY8/vl1dG5Zn3cHT3Dd1DV/c0QpPD/3eJCIi4qxly5YxcOBAAgICABg0aBBLly6lRo0aNG7cmLi4OEJCQujcuTPGGBo3bsy+ffuAjDO+c+bM4c033wQgKSmJAwcOsGzZMh544AEAGjVqRGhoaNb7zZgxg/Hjx5OamsqRI0fYsmXLRe333HMPHTt2pEOHDhfFuWjRIj755BOWLVv2t2N64IEHuPrqq+nbty/ff/89W7ZsoV27dgAkJyfTpk2brL7XX389kHGWOiIigrJlywIZ01qWLFmCp6fnZY/fvn07FStWJDw8HIDixYsTFxfHggULGDlyJJ6eGbVKqVKlnPk2FAkqtvOzfu/CRx1h/zJY+AJ0e+mi5vs71+KOSWt4Y952nuzZwE1BioiI/Hsbb9noVL/BdQczuO7gf/1+OS004OPjk/Xc4XBkbTscDlJTU7OO/frrr6lXr55Tr7l3717efPNNoqKiKFmyJLfeeitJSUlZ7S+88AInTpzgo48+uui4DRs2MHz4cH788ce/va35pEmT2L9/P++//35WLF27dmXq1KmX7f/7Lxk5xZzT8Rs2bMAYc9n+l9sveTyNxBjTwxiz3RizyxjzxGXa6xtjVhhjLhhjHsnNsYVSYAUYPAmMB/z6LmyZA0BaehpjosfwQOR9vHVDKN+vP8JPm464N1YREZECpGPHjsyePZvExEQSEhKYNWvWn84q56R79+689957WYXq2rVrAWjfvj0zZswAYMuWLWzcmPELxLlz5wgICCAoKIhjx47x448/Zr3Wxx9/zM8//8zUqVNxOP4oyw4cOMCgQYOYPHkydevW/ct4oqOjefPNN/niiy+yXqN169YsX76cXbt2ARnzyXfs2PGnY1u1asXixYs5efIkaWlpTJ06lauuuirH4+vXr8/hw4eJiooCMm7VnpqaSrdu3Rg3blzWLySnTp1yKpdFQZ6d2TbGeABjga5ADBBljJljrd2Srdsp4H5gwD84tnCq3ha6vgjznoLZ93C2ZDUe3zyO5YeWAzB5+4d8MPR2bpsURZ3ygdQqW8zNAYuIiOR/zZs359Zbb6Vly5YADB8+POuCxL/zzDPP8OCDDxIaGoq1luDgYL7//nvuuecebrnlFkJDQ2nWrBmhoaEEBQVRp04dmjVrRkhICDVr1syamgEZ862rV6+eNcVj0KBBPPvss7z44ovExsZyzz33AODp6Zl1seWl3n//fU6dOpV1YWRYWBgff/wxkyZNYsiQIVy4kHF918svv/ynwr1ixYq89tprdOrUCWstvXr1on///gA5Hj99+nTuu+8+zp8/j5+fH7NmzWL48OHs2LGD0NBQvLy8uPPOOxk1apRT+SzsTF6t12yMaQM8b63tnrn9JIC19rXL9H0eiLfWvpnbY7MLCwuz2T+YkZGRREREXInh5C1rYcbNsHUOT1StyQ+eF9/U5r8d/suZE42Y9OteZt3TjgCff/87VIHNlRsoV85TrpynXDlPuXJefsrV1q1badAg/06BjIuLIzAwMFfHpKWlkZKSgq+vL7t376Zz587s2LEDb29vF0WZP/yTXBVkl/vsGmOirbVhl+ufl9NIKgMHs23HZO5z9bEFnzHQfyyUb8zDje6krF/ZrKayfmWpHFiZIS2rElqlBE98oxveiIiIuENiYiLt27enSZMmDBw4kA8//LDQF9ry9/LyzPZgoLu1dnjm9k1AS2vtfZfp+zwXn9nOzbEjgBEA5cuXbzFt2rSstvj4eIoVK8DTLGwaGA/2XdjHO0ffoYp3FYaXHU6QZxAAyWmWl1cm0aGyJ12Dvf7VWxX4XOUh5cp5ypXzlCvnKVfOy0+5CgoKonbt2u4OI0dpaWl4eHi4O4yLxMbG0q9fvz/tnzNnzt9eQOlK+TFXrrRr1y7Onj170b5OnTrleGY7L1cjiQGqZtuuAhy+0sdaa8cD4yFjGkn2P5flpz+f/VuNtxajsfHDu37vi/bXb5rIoA+XMzAilLDgf77sTmHKlaspV85TrpynXDlPuXJefsrV1q1b8/XUg/w4NSIwMPBP63PnB/kxV67k6+tLs2bNnO6fl9NIooA6xpgaxhhv4AZgTh4cW/jE7qbFN/fj/fVwOLL+oqZqpf1549omjJqyluNxSTm8gIiIiPtp2qMUNP/kM5tnxba1NhUYBfwMbAVmWGs3G2NGGmNGAhhjKhhjYoCHgaeNMTHGmOI5HZtXsec7pWpCw/6QkghTboC4oxc1d6pfjuvCqzJqylpS0tJzeBERERH38fX1JTY2VgW3FBjWWmJjY/H19c3VcXl6Uxtr7Vxg7iX7xmV7fpSMKSJOHVtkGQN934bTe+HACpg6BG6bC15+pKan8nb023Rq3Jn1Bz14/adtPNW7obsjFhERuUiVKlWIiYnhxIkT7g7lspKSknJdVBVVRSlXvr6+VKly2VI1R7qDZEHl6QPXfwETrobDa2D23ZzpM4ZHljzGqqOr+GHvD4zvP5nbPt5B4yol6NekkrsjFhERyeLl5UWNGjXcHUaOIiMjczUvtyhTrv5ant5BUq6wgDJw43TwDiRhy2yGfNObVUdXAXDy/EmeXfEoY4c25vk5m1lz4LSbgxUREREpelRsF3TlGsDgTwnw9KdH6aYXNZ1MOkmJ4om8cW0oIydHE3M60T0xioiIiBRRKrYLgzpd4cGNjOr2Ph0qdwCgebnmTOs9jerFq9O5QXlGdKzJHZNWE5eU4uZgRURERIoOFduFRUBpPBwe/Lfjf7mn5kA+bvUcpf3+WOD+jvY1aBFckvunriUtXVd+i4iIiOQFFduFTOChtdy9eDxe02+GC/FZ+40xvNAvhJQ0y8s/bHFjhCIiIiJFh4rtwqZ8IyheCY5thG9GQPof62x7eTgYO7Q5S3acYPLK/W4MUkRERKRoULFd2PiXghtngG8J2P4DLHw+qykpNYlPNr/H+0Mb8s6CnSzZkT/XNhUREREpLFRsF0ZlasN1n4PDE5a/A2u/4FTSKe6Ydwefbv6UtzY8zbtDQnlo+jp2Hotzd7QiIiIihZaK7cKq5lXQ600Ajsx9mJu+vYYNJzYAsOLICuYeeZcne9bnjs9WExt/wZ2RioiIiBRaKrYLs7DboN2DlDAeBHn6X9S0OXYzXRsXp2+Titw1OZoLqWluClJERESk8FKxXdh1eR6/u5bxXq/PqVKsCgCtKrbi856fE+QTxH+61qNsoA9Pfr0Ra7UkoIiIiMiVpGK7sDMGSteitF9pPuzyITdUuooPWz5LoHcgAA6HYcx1Tdl1Ip6xi3a5OVgRERGRwkXFdhESfHIPT62cjteUG+D8maz9ft4efHxzGFNWHeCHDUfcF6CIiIhIIaNiuyip2BRKBsPxzTDtRkhJymoqV9yXCbeE8ey3m1hz4LTbQhQREREpTFRsFyX+pWDYNxBYCfYvh6/vgPSMCyMTUhL4LuZDXhpUixGfR3MoPv1vXkxERERE/o6K7aKmRFUY9jX4BsG27+GHhzmRcJzbfrqNL7d+yaxDr/F4z1r8b3USh86cd3e0IiIiIgWaiu2iqHxDGDIdPH05uH4yQ2f1Y+uprQCsPLKSdec/pluwJzd/sopTCcluDlZERESk4FKxXVRVbwPXfkoZz0DK+pe9qGnn6Z1cVS2N7iEVuG1SFAkXUt0UpIiIiEjBpmK7KKvfC78H1vNer8+oGlgVgI5VOjKpxyT8HH482r0e9csHMvKLaJJTNYdbREREJLdUbBd1fiUo5VuKD7t8yK2VIninxmD8vTLuNmmM4ZWBjfDz8uA/X60nPV03vRERERHJDRXbAkD1s8f5z69f4jn9Jji6MWu/p4eDd4c04/i5JF74brPuMikiIiKSCyq2JUPl5tCgD1w4B5MH4p8Qk9Xk6+XB+Jtb8Nu+07z/i+4yKSIiIuIsFduSweEBgyZAzU6QcIIm65+G2N0AnE46zb2Lbue+3mnMXBPDl6v2uzlYERERkYJBxbb8wdMHbpgCwR3wST4Nn/Xl+JFobvvpNtafWM9zKx/l8QHevLtwJ3M36rbuIiIiIn9HxbZczNsfhkzjbPEGnI87zK0/3cbusxlnuJPSkngp6hFeuLY0z8zexK+7Tro5WBEREZH8TcW2/JlPMTaEPotfzU4MrXf9RU0hZULoWKMeY4c2576pa9kYc9ZNQYqIiIjkfyq25bLSPP3hplkMbfsU9zW7D4BOVSIY23ks/l7+tK5ZmlcHNeb2z6LYdTzezdGKiIiI5E+e7g5A8r87G99JtbPH6bz5R7xax0GADwDdQyoQl5TKsI9XMeXOVtQsW8zNkYqIiIjkLzqzLX/LpKfRY9NPeB3dBJ/3h8RTWW3XtqjCQ13rMPTjVew7meDGKEVERETynzwtto0xPYwx240xu4wxT1ym3Rhj3s1s32CMaZ6t7QFjzCZjzGZjzIN5GXeR5+EJw76G0rXh2CaYPADOn8lq7hEaxJ1XVWTox6s4EJvotjBFRERE8ps8K7aNMR7AWKAn0BAYYoxpeEm3nkCdzMcI4MPMYxsBdwItgSZAH2NMnTwKXQACy8Mt30HJGnBkPXwxCJLOcfL8SW77+Tbmn3qJW9uXY8iElRw8pYJbREREBPL2zHZLYJe1do+1NhmYBvS/pE9/4HObYSVQwhhTEWgArLTWJlprU4HFwMA8jF0AilfKKLhLVIND0Rz5ciC3/ngzO0/vZFPsJuaffpGhbctw48crOXTmvLujFREREXG7vCy2KwMHs23HZO5zps8moKMxprQxxh/oBVR1YaySkxJVMwru4pV5KzmG/XF/fLu2ntrKhYCF3NImmBsnrOTIWRXcIiIiUrQZa23evJExg4Hu1trhmds3AS2ttfdl6/MD8Jq1dlnm9kLgMWtttDHmDuBeIB7YApy31j50mfcZQcYUFMqXL99i2rRpWW3x8fEUK6YVM5zxd7nySzxMSsop3jy/gN0XMm5609ivMbeXvR1P48ncvcksPpjKEy19KelbuK/D1efKecqV85Qr5ylXzlOunKdcOU+5gk6dOkVba8Mu15aXS//FcPHZ6CrAYWf7WGs/AT4BMMa8mtn3T6y144HxAGFhYTYiIiKrLTIykuzbkjNnc9U85XZG/TIKvwsJvNX5HbwDygIQEQEfRO7ivegYpo1oTblAX9cG7Eb6XDlPuXKecuU85cp5ypXzlCvnKVd/LS9POUYBdYwxNYwx3sANwJxL+swBbs5claQ1cNZaewTAGFMu82s1YBAwNe9Cl5z4e/kztuoAxqybj/cX1160LOA9EbUZ2LQyN05Yxcn4C26MUkRERMQ98qzYzrywcRTwM7AVmGGt3WyMGWmMGZnZbS6wB9gFTADuyfYSXxtjtgDfAfdaa0/nVezy1/wqNcUnqBocWQeTekPcsay2+zrXoXfjigydsIpYFdwiIiJSxOTpZFpr7VxrbV1rbS1r7SuZ+8ZZa8dlPrfW2nsz2xtba1dnO7aDtbahtbaJtXZhXsYtf6NENbjtRyhTD45vgUm94OwhAI4lHGM7Y2hbz5OhH6/idEKym4MVERERyTuF+8o1yTvFK8Jtc6FCY4jdBZ/24Pjh1dwx7w6WH17OyvMvEV7bMOwTFdwiIiJSdKjYlisnoEzGsoCVwzh1LoY7fr6d/ef2AxATH8NvF16mRU0Prh+/guPnktwcrIiIiIjrqdiWK8uvJNw8m8BaXalTvsVFTQ1KNeCZXuH0b1qZwR+t0J0mRUREpNBTsS1Xnk8gXkNn8N8e4+ke3B2AiErteKPjG3g5vLi3U21uaxvM9R+tYPeJeDcHKyIiIuI6KrbFZbwcXozuMJpHqnTjf9Fz8YrJut6VW9vV4KGudblh/Eo2Hz7rxihFREREXEfFtriUp8OTW87G4Z10Fr4YBLsXZbUNDqvKC/1CuGXib0Tv10qOIiIiUvio2BbXG/AhNB0KKYkw5XrYNjerqUG1JG7qHMeIz1ezfNdJNwYpIiIicuWp2BbXc3hAv/chfDikXYDpQyF6EkcTjjJi/ggm7nyBG64+zP1T17Jgy7G/fz0RERGRAkLFtuQNhwN6vQlXPQ42nTM/PMRdcwZzNOEo6Tadz3e9yYBO23j8mw18u+6Qu6MVERERuSJUbEveMQY6/R/0eYv1gSXZn3LuouY0Ryxf3tGK1+ZuY8qqA24KUkREROTKUbEteS/sdq4aGc3bnd7Bx8MHgC7VOvN066epX7E400a05oPIXUxYssfNgYqIiIj8Oyq2xT38ShJRNYJxXcbROaguo3dvwvN8xhKAwWUC+GpkG6ZGHWDM/B1Ya90crIiIiMg/o2Jb3CqsTGPePrAHn5gomNgdTmfc3r1ikB8z7mrDwq3HeGr2JlLT0t0cqYiIiEjuqdgW9/L0gVvnQvlGELsTPukKRzYAUKaYD6/dUI69p49w1+RoEpNT3RysiIiISO6o2Bb3K14RbpsLwR0g/hh82gv2LCYmLoYHI+/mTNBb+PjFMmT8Sk7GX3B3tCIiIiJOU7Et+YNvEAz7GhpdA8lxnJwymBE/DOXE+RPExMewPv0l6gWf4JoPf2XvyQR3RysiIiLiFBXbkn94+sCgj6HNKCaVLMHBC6eyms4lnyO4yhHuvqoW1320gjUHdHt3ERERyf9UbEv+4nBA91d4YOgCBtQekLW7V3Av7m5yNze0rMbr14Qy/LPVzNt81H1xioiIiDjB090BiFyOV4lqvNj2RaoGVuXX7V/z0v7tmAvnwDeITvXLMem2cIZ/tpqj55K4uU2wu8MVERERuSyd2ZZ8yxjDiHo3MmH/Prz3LIKPu8KpjBvdhFYpwcyRbZm0fB+v/biV9HStxS0iIiL5j4ptyd98iuF1x89QtgGc3A4TOsO+5QBUK+3Ps9f68ev+7Tw4fR0XUtPcHKyIiIjIxVRsS/5XMhjumAd1usH5U/B5f1j7BXvO7OH/VjzE6RL/40TKVm6Z+Btnz6e4O1oRERGRLCq2pWDwLQ5DpkHreyE9hdjv7uOeucOIS47j7IUzbDdvElBqI4PH/crBU4nujlZEREQEULEtBYnDA3q8Cn3f4ZsyFTiUEpfVlJKeQpu6Dm4Ir8bAD35l1Z5YNwYqIiIikkGrkUjB0+JWhje6ltStk/lg/QcAXFOjN3c0vgNjDLXLFeOeL9fwaPd63NCympuDFRERkaJMZ7alQDI+xbi76d282v5VrgqoxlOrZmJiogDoWLcsM0a2YfySPbzw3WZS09LdHK2IiIgUVSq2pUDrW6M37yV64ZVwAib1hg0zAKhVthiz7mnHruPx3DYpirOJunBSRERE8p6KbSnYHA7M0K8g7A5IS4Zv7oT5z0JaKkH+XtwQEYtH8d8Y+MFydp+Id3e0IiIiUsTkutg2xgQYYzxcEYzIP+LhBX3GQM83wHjA8nfgi0FsO7SS51c8R3TieKrV/YnB45ayZMcJd0crIiIiRcjfFtvGGIcx5kZjzA/GmOPANuCIMWazMeYNY0wd14cp4oRWI+CWORBQltP7l/LAgntISksCIPr0D9QJncFDM9YxcdlerNUdJ0VERMT1nDmzvQioBTwJVLDWVrXWlgM6ACuB0caYYS6MUcR5we1hxGJ21GzHaY+LP943NOzH7HvaMWP1QZ78ZiPJqbpwUkRERFzLmWK7i7X2JWvtBmttVnVirT1lrf3aWnsNMN2ZNzPG9DDGbDfG7DLGPHGZdmOMeTezfYMxpnm2tocyz6ZvMsZMNcb4OvOeUgQFVabVsO+Z3OtLKherDMCQkk0YWKMXVUv5M/PutsQmJDPs41XExl9wc7AiIiJSmP1tsW2t/dtlHJzpkznPeyzQE2gIDDHGNLykW0+gTuZjBPBh5rGVgfuBMGttI8ADuOHv3lOKtnql6jG191RuC2rEo2u+g097wtlDFPPx5KNhLQivUZL+Y5ezMeasu0MVERGRQuqfXCDZ1RgzwRjTNHN7hJOHtgR2WWv3WGuTgWlA/0v69Ac+txlWAiWMMRUz2zwBP2OMJ+APHM5t7FL0lPQtycMtH8MrqCocioaPOsLepTgchke71+f/ejXg5knLmPbbAc3jFhERkSvO5LbAMMbMAm4DngbmAtdaa+9x4rhrgR7W2uGZ2zcBray1o7L1+R4Yba1dlrm9EHjcWrvaGPMA8ApwHphnrR2aw/uMIOOsOOXLl28xbdq0rLb4+HiKFSuWq/EWVYUtV17J52iw9U1KnV6PxcHuWrcSU6UfqxJ+47vTc0k5ciN1/YK5qaE33h4mV69d2HLlSsqV85Qr5ylXzlOunKdcOU+5gk6dOkVba8Mu1/ZPbtd+wlp7BnjEGDMaCHfyuMtVMJdW+pftY4wpScZZ7xrAGeArY8wwa+0Xf+ps7XhgPEBYWJiNiIjIaouMjCT7tuSsUOaqS2/45SXMsreovXsiSV6HmJGym+T0ZDwrjiPe3MBbm9rw0bAwqpX2d/plC2WuXES5cp5y5TzlynnKlfOUK+cpV3/tn9zU5offn1hrnwA+d/K4GKBqtu0q/HkqSE59ugB7rbUnMueHfwO0zWXcUtQ5PKDL83DdZGJ9A3kgaTvJ6ckApKancth8Q6+mfgz8YDkLthxzb6wiIiJSKOS62LbWfnvJ9ntOHhoF1DHG1DDGeJNxgeOcS/rMAW7OXJWkNXDWWnsEOAC0Nsb4G2MM0BnYmtvYRQBo2I/it/1El1p9L9r9QtsXuO+qcMbf3IJnvt3EGz9vIy1d87hFRETkn/snF0h2+ScXSFprU4FRwM9kFMozrLWbjTEjjTEjM7vNBfYAu4AJwD2Zx64CZgJrgI2ZcY/Pbewiv/Mq34gnO7zM/676HwFeAVxvStBj/RxITqBF9VJ8d1971uw/wy0Tf9PygCIiIvKP/ZM52/eSeYGkMaYU0NTZA621c8koqLPvG5ftuc18/csd+xzw3D+IVyRH3YK70cB6Uu7LGyF1A8REweBJlCkfwuQ7WvK/+Tvo+94y3h/anObVSro7XBERESlgnLld+zuZX/0yd52w1p6x1j4CdMP5CyRF8qWqNa7GZ8QiKFsfTu6ACVfD6k/xdBge616P8OYruWPK93y+Yp+WBxQREZFccWYaSefMr8syv/7TCyRF8q9yDeDORdDsJkhNgu8fhJm3MX3TJH45Og1H5XeYsHYGD05fR/yFVHdHKyIiIgWEM8X2T8aYFUAFY8ztQEz2W6Xn4gJJkfzN2x/6vw+DPgbvYmze+QOvr3kbgOT0C5wJmEyM/Y4+7y7VXSdFRETEKc7crv0RYCiQRsY6188AG40xm40x010cn0jeCx0Mdy1hYa1WpJCetdvf05//9bqZ/3Srx62f/saEJXtI12olIiIi8hecukDSWrvHGNPFWrvj933GmGJAI5dFJuJOpWtx//XfU3XnLF5d9SpJaUm84BNMsFdxgpuUomnVEtw/bS3Ldp3kmsoquEVEROTynLlA0gBkL7Qzt+OttSuz9xEpbAbWGciU3lN4wKsKPTb9BOM6wL5lVC3lz4y72hBSqTjP/nqeZTtPujtUERERyYecmbO9yBhznzGmWvadxhhvY8zVxpjPgFtcE56I+9UpWYfhPcdB5TA4FwOT+sC8Z/CyKTzWoz43NU7hwTlTGP3jNlLS0v/+BUVERKTIcKbY7kHGfO2pxpjDxpgtxpg9wE5gCPCWtXaSC2MUcb+S1eH2n6DjY2AM/PouTOiMPbqZ6PTpXCgzgXnHPuSacYs5eCrR3dGKiIhIPuHMBZJJ1toPrLXtgOpkLAXY3Fpb3Vp7p7V2nauDFMkXPLzg6qfg9p+hZA04tpGvpvZm4/mNAJwwv3C25Ov0G/cDc9YfdnOwIiIikh84fbt2Y8woIMBae8Rae8Z1IYnkc1VbwshlxDUbylulL76rZJWgUkwcdjVj5m3nsZnrSUzWmtwiIiJFmdPFNlABWG2MmWGM6aGLIqVI8ylGYP8P+KjHJMp4lgHA38OH0eWvolm10nx/fwdS0y193lvG+oNn3BuriIiIuI3Txba19mmgDvAJcCuw0xjzqjGmlotiE8n3Qss14fGKjzOgRm+ePHueqt8/CjPvoFh6HGOua8pDXepyx2dRjJm/QxdPioiIFEG5ObONtdYCRzMfqUBJYKYx5nUXxCZSIPg6fHmpw2sMaPUIePnDppnwYTvYE0nfJpX44f4OrIhZS78P5rPzWJy7wxUREZE8lJs52/cbY6KB14HlQGNr7d1AC+AaF8UnUjAYA2G3wchlmUsEHoLP+8NP/4eP4zQnA8ZxptSrDP58Ih8v1Z0nRUREiorcnNkuAwyy1na31n5lrU0BsNamA31cEp1IQVO6VsZqJRH/B8YDVo7ltZkDOHn+JPGpp0gr9zGTdrzBDRNWaIlAERGRIiA3xbaPtXZ/9h3GmP8CWGu3XtGoRAoyD0+IeByGz2dphTrMNRcX1dc3D+Hq+uXpP3Y5M1YfJGN2loiIiBRGuSm2u15mX88rFYhIoVO5Ba1vX8qI0BF4GA8A6vuUZmRQI0ZeVYspd7bi0+X7uPPzaE7EXXBzsCIiIuIKf1tsG2PuNsZsBOoZYzZke+wFNrg+RJGCy8vbj/ua3ccXvb6gfmB1Xt6zBa/P+8MP/6F+ScO397ajXoVi9HxnKT9tOuLucEVEROQKc+bM9hSgLzAn8+vvjxbW2mEujE2k0GhUphEz+n5NvTYPgcMToj6GD9rivX8xj3avz0c3teCFX6Zx37SVnE5Idne4IiIicoU4c7v2s9bafdbaIdba/dkep/IiQJHCwnj5QKcnYUQkVAiFswdg8gCYcx/pZj3xQROJSnmKzh9O5IcNRzSXW0REpBBwZhrJssyvccaYc9keccaYc64PUaSQqdAY7vwFrn4aPLxJXDuZZxY/CkBC+nFSyo7lpeXvctfkaI6dS3JzsCIiIvJvOHNmu33m10BrbfFsj0BrbXHXhyhSCHl4QcdH4a4lTAwO5VD6H0W1wfDegMHUrxBIr3eWMj3qgM5yi4iIFFC5uanNZ8aYEtm2SxpjJrokKpGiolwDbr/xJ26sf2PWrutNEOHHd/Bw17p8MbwVX6w8wNCPV3EgVutyi4iIFDS5Wfov1Fp75vcNa+1poNkVj0ikiPH38ufJVk8yqcckWgfW5IG9G2HWXfB5Pxp4HWfWPW2JqFeW/mOX8fHSPaTp7pMiIiIFRm6KbYcxpuTvG8aYUoDnlQ9JpGhqUb4FEwbOplj/D8CvFOxdAh+2xXPZm4xoW4Xpd4Uxaeer9PnoK3Yci3N3uCIiIuKE3BTL/wN+NcbMBCxwHfCKS6ISKaqMgaY3Qp3uMP8ZWPclLHoFNszg50adOeuxijj/1QyevoKbGwznvk4N8PbMze/MIiIikpec/l/aWvs5cC1wDDgBDLLWTnZVYCJFWkBpGPAB3PI9lK7DznP7mHjgJwDSSSM9aAHzDn1Bn/eW8ttercIpIiKSX+VqGoi1djOw2UWxiMilanSAu5fz08/3kHryt6zdZX1KMmPA4yzfeZ4Hpq2lXe0yPNmzPqWL+bgxWBEREbmU1tkWye88fRjV62PGRIyhjF8ZAJ48fICgacPoVTGB+Q9fRZCfF93eWsK03w6QrgsoRURE8g2n1tk2xhggROtsi7iHMYau1bvy7YBvebr29XRJdcC+pfBhG4otfYVnulXn8ztaMn31QQZ89ANbj+j3YBERkfzAqTnbNuOOGrP+7ZsZY3oYY7YbY3YZY564TLsxxryb2b7BGNM8c389Y8y6bI9zxpgH/208IgVNce/iXN/uacy9UdB0GKQlw7Ix8H44IacW8Np1pdnv/ww3fP0fnpmzivgLqe4OWUREpEjLzTIGK40x4f/0jYwxHsBYoCfQEBhijGl4SbeeQJ3MxwjgQwBr7XZrbVNrbVOgBZDIFSj+RQqsgNIwYCzcMR8qNoVzh7Azb+fVeSNIJ5X0wJX8cPpBIj54mx83HtEdKEVERNwkN8V2J2CFMWZ35lnnjcaYDbk4viWwy1q7x1qbDEwD+l/Spz/wuc2wEihhjKl4SZ/OwG5r7f5cvLdI4VS1Jdz5C/R5m+9KV2BtyumsphTiubN9XcbM38Htk6J0B0oRERE3MM6e8TLGVL/cfmeLXmPMtUAPa+3wzO2bgFbW2lHZ+nwPjLbW/n5R5kLgcWvt6mx9JgJrrLXv5/A+I8g4K0758uVbTJs2LastPj6eYsWKORNukadcOS+/5GpXwhamnvma46nHAWiVFsB/fDtzqPzV/Lw/jR/3ptC1uhc9a3jh7WHcEmN+yVVBoFw5T7lynnLlPOXKecoVdOrUKdpaG3a5ttws/XePtfbx7DuMMf8FHs+h/6Uu97/7pZX+X/YxxngD/YAnc3oTa+14YDxAWFiYjYiIyGqLjIwk+7bkTLlyXn7JVQQR3JJ2J59v+ZxPN0zguYM7qZq6lQYJK+nS63Ue8G/AKz9s5cXVZ3m6d326h1Qg49rnvJNfclUQKFfOU66cp1w5T7lynnL113IzjaTrZfb1zMXxMUDVbNtVgMO57NOTjLPax3LxviJFhpeHF3c0voP51/1C1X4fQrEKcGg1TOhMlaWP8+HA6rw2KIRnfn2S/p9MYvtR3fZdRETElZxZZ/tuY8xGoF7mXO3f52vvBTbm4r2igDrGmBqZZ6hvAOZc0mcOcHPmqiStgbPW2iPZ2ocAU3PxniJFkp+XP4ReB/ethrb3g8MD1nwO7zbn0O6XSPKJZq/XGAbPGsF/vlnImcRkd4csIiJSKDlzZnsK0JeMQrhv5qMP0MJaO9TZN7LWpgKjgJ+BrcAMa+1mY8xIY8zIzG5zgT3ALmACcM/vxxtj/Mk4u/6Ns+8pUuT5BEK3l+DuFVDras6mnGPsiRVZzdZ/E9uSZtFlzGImr9xPmm6IIyIickX97Zxta+1Z4Kwx5jZgEBD8+3HGGKy1Lzr7ZtbauWQU1Nn3jcv23AL35nBsIlDa2fcSkWzK1oVh37Dk19Gc2TUla7efw5tPWw/kZEQ4z3+3mS9X7uf5fiG0rql/aiIiIldCbuZszyZjab5UICHbQ0QKAmPo2+5Jvuj1BaFlQgG4I+485T4fRMOVjzL9+iqMuro2/5mxnnu/XEPMaS0VKCIi8m/lZjWSKtbaHi6LRETyRJOyTZjcazLzds0hYm8UxI6HDdMwW76lT9v76DxqFONWHKP3J2MZGNKKR69uS4BPbn5UiIiIyO9yc2b7V2NMY5dFIiJ5xmEc9KgzAN9ur8CoKAgZCKnnYcnr+I0L58bi8/GqMINvjj9A+wkP8cmvm0lNS3d32CIiIgVObort9sAaY8z2f3gHSRHJj0oGw+BJcPs8qNwC4o/xzoaPuJCWhCWV1MBIPtz+MN3fWcLCrcd063cREZFcyM3fhnOzpraIFDTVWsEdC9ge9QHfbfvooqYnQnpQolgDXp27lY+X7uWp3g1oVDnITYGKiIgUHLk5s30A6ADcknmLdguUd0lUIuIeDgd1Wt7Dq+1fpbx/xj/vmh4BDJj7PF02P8HPN1ehT5OK3DYpigenrdVFlCIiIn8jN8X2B0AbMm4sAxAHjL3iEYmIWzmMg761+vL9wO95qMVDPFKiGZ4ePrD5Gzw/bMXQE28TObI+1Ur503vcTP4z5zvOnk9xd9giIiL5Um6K7VbW2nuBJABr7WnA2yVRiYjb+Xr6cnuj2+nQ50O4fw00GwY2HVZPJGBcOA85ptG08ULmnf4/Ok66nTGLlpOcqosoRUREsstNsZ1ijPEgY/oIxpiygP5nFSkKgqpA/7EZd6Ks3wdSEvk1+gPWxK4BIN1/HZ/uv5eId6bzzZoY3YlSREQkU26K7XeBWUA5Y8wrwDLgVZdEJSL5U7n6cMOXpN/+M29VrXtRU6sStfhfv6v5YuV+er2zlHmbj2rlEhERKfKcLrattV8CjwGvAUeAAdbar1wVmIjkX45qrflPxOs0KNUga99DWxbT5ue+fN3xGI91r8OY+TsY+MGv/LrrpBsjFRERca9c3RbOWrsN2OaiWESkAGlTqQ2tKrbix70/sn3PPEJO/wInt2Nm3kbnciF06vYE311ozhOz1uJfejUvXn2bu0MWERHJc06f2TbGfGaMKZFtu6QxZqJLohKRAsFhHPSu2ZuHu7wDo6Khz1tQvDIc34xjxk30XzWEe1ot55DHFIb/cg0vblnIliOx7g5bREQkz+RmznaotfbM7xuZq5E0u+IRiUjB5OkNYbfD/Wuh5xtQrAIXjm1gwv4fALAe5zgRMJshX73IwzPWcfCU1ugWEZHCLzfFtsMYU/L3DWNMKXI5DUVEigBPH2g1Ah5Yx8zwIRxPT8pq8sKDbyOuokqQL33fX8bTszdy+Mx5NwYrIiLiWrkptv8H/GqMeckY8yLwK/CGa8ISkQLPy4/eV4/m9ka34+vhC8Dg+GSCv7mNhw8/zOLB3gT4eNLznaU8M3uTim4RESmUcrMayefANcAx4AQwKHOfiMhllfAtwUMtHuLHa37k5vo30qNYO/AtAfuXEzS9P08ee5Qlgz3w93LQ49MXuGXGWGJOx7k7bBERkSsmVxdIAoette9ba98DjugCSRFxRhm/Mjza6knOVh8KD26AiP8Dn+KwbylBMwZxx9H78SozjzXnx9Hz677cPOM9Fd0iIlIo6AJJEclbvkEQ8Tg8uBE6PQV+Jfk05SjJ6SkZ7V6xbL/wNb3eXcqz327iyFlNLxERkYJLF0iKiHv4lYCrHiNx1G/MKuZ3UdPDKSks6x2Pvyf0eFtFt4iIFFy5vUByxSUXSL7umrBEpKjwDyjH1wNmc23da/F0eFLOw5+BBzZR/Ps7eWLvbSzrcQx/T0uPt5fyf99sYPeJ0+4OWURExGm5vUByEBdfIDnZVYGJSNFRuVhlnmvzHD8M/IHRncbg3ectKFENTu4g8MdRPLFrGMu6HiDFsZUB3/Vi4Jcvsv7QYXeHLSIi8rdyc4GkAZoDpTIvkIw3xrR0WWQiUuRUKlaJ8MrtMm6Oc98a6P8BlKoFp/cROO8/xMW/Bx7x7Er9imHz+jPo83fZdOisu8MWERHJUW6mkXwAtAGGZG7HAWOveEQiIgAeXtBsKIyKgms+YUv5uixPy1ZYO5JoVaY0wz9bzS0Tf+O3vafcF6uIiEgOclNst7LW3gskQdZqJN4uiUpE5HcOD2h8LXu7PUeAV0DW7uYmgMejH2J5yLdcW+MCj3y1nuvGrSBy+3GstW4MWERE5A+5KbZTjDEegAUwxpQF0l0SlYjIJXrX6sPP1/zMfc3uo6RPSe7wKAtpF/BY+xl9F/clsvqn3FvvLK/O3UqXceP4YMV80tNVdIuIiHvlZum+d4FZQDljzCvAtcDTLolKROQygnyCGBE6gpsa3pRxC/iTO+HXd2D9dBxbv+Wqrd/SNrg9Xb0S+HBHLBM216Rf8BAe6zCIAB/9IU5ERPJeblYj+RJ4DHgNOAIMsNZ+5arARERy4ufphzEGytaF/mMz7krZ9j7wDuTn2PXEpsYCkOq1h9mH3qDj/+by9oIdxMZfcHPkIiJS1ORmGgnW2m3W2rGZt2zf6qqgRERypXgl6PYy9sGNfFK1wUVNg60385ptIe7UUTq9GclTszay50S8mwIVEZGiJlfFtohIvuZXgv/r+CrtKrXL2jUsZgelV/2XZ3YMZlXod9QxMQwet4IRn6/m523bdDGliIi4VJ7ebt0Y0wN4B/AAPrbWjr6k3WS29wISgVuttWsy20oAHwONyLhI83Zr7Yq8i15E8jtjDOEVwgmvEM72U9uJOvobwe2rwMoPYOc8/DZM5lYmc1PNq5ldvDuPrPwQv19rc33dYdzfph9eHh7uHoKIiBQyeVZsZ65kMhboCsQAUcaYOdbaLdm69QTqZD5aAR9mfoWMIvwna+21xhhvwD+vYheRgqdeqXrUK1UvY6NWp4yLKVd+COun4rHnFw6W3w7+Hpz32MGk3c/y5Zbp3BfyGte2qEKgr5d7gxcRkUIjL6eRtAR2WWv3WGuTgWlA/0v69Ac+txlWAiWMMRWNMcWBjsAnANbaZGvtmTyMXUQKujJ1oM8YeGgziVc/xYxiF/++PrK4YdeuHbT/7yJe+G4z+2MT3BSoiIgUJnlZbFcGDmbbjsnc50yfmsAJ4FNjzFpjzMfGmABERHLLvxSnmt1Ig3JNsnYFOXwYtv1HXtk3hJW1P6d24gYGjl3O8M+iWL7zBOnpuqWAiIj8MyavLg4yxgwGultrh2du3wS0tNbel63PD8Br1tplmdsLyVhu0AArgXbW2lXGmHeAc9baZy7zPiOAEQDly5dvMW3atKy2+Ph4ihUr5qohFirKlfOUK+flt1wdTD7IonOLqJRmuP/kEcqeWIHJvFfXuYBgIv178N/E0sSXmE8zn45cW7Elxbx88iS2/Jar/Ey5cp5y5TzlynnKFXTq1CnaWht2uba8LLbbAM9ba7tnbj8JYK19LVufj4BIa+3UzO3tQAQZF0SutNYGZ+7vADxhre39V+8ZFhZmV69enbUdGRlJRETElRtUIaZcOU+5cl6+z9W5w7B6IkRPgoQTANxTK4Sl6XEZ7Wn+dCg1gmc63UjFID+XhpLvc5WPKFfOU66cp1w5T7kCY0yOxXZeTiOJAuoYY2pkXuB4AzDnkj5zgJtNhtbAWWvtEWvtUeCgMSbzaic6A1sQEbmSileCq5+GhzbDwPHsq9zkj0IbwCORanHH6PHWEkZNWcPqfae0dKCIiPylPFuNxFqbaowZBfxMxtJ/E621m40xIzPbxwFzyVj2bxcZS//dlu0l7gO+zCzU91zSJiJy5Xj6QJPrWePvhWPFC6TbjKkljRwBPL77vzxaui6rzDU8+1Uo1juQW9pUp3/Tyvh5a+lAERG5WJ6us22tnUtGQZ1937hszy1wbw7HrgMue3peRMQVBtUZROuKrZm+fTozd8xkaIlmmONn8IjdQdvY1/jBuxhHKvZj/PoIXlu+kRoVz/NQq1u5qkZjd4cuIiL5RJ4W2yIiBU2lYpV4qMVDjGwyEk+HJ3QDtn0Pq8ZjDvxKpV1TeJopLKndgN0XEhi1ZB7FFtfj/tAnua5pczwcxt1DEBERN9Lt2kVEnODn6YeXwws8vCBkINz+I9y9AlqOYHFQaU6k/bEu9wWzhwXLNhPx5iI+Wryb0wnJboxcRETcSWe2RUT+qfINodcb/OiXBgfmZ+3ul5LO86fuI65SB2bu7M7Vi2rROaQyQ1tVo2nVEhijs90iIkWFim0RkX/ptY7/pcvBbkzbNo3oY9FcVzIUjh4h8PBSbmMptxSrwJqEfrwyNZxDJZfSonIwj7W/iWpB5dwduoiIuJiKbRGRf8nLw4sewT3oEdyD/ef2U714dTh/GtZNhdUTccTuJCx+PO97TqS7TxUWn4xk8azJVPQK45V2LxEeXNHdQxARERdRsS0icgVVL14944lfSWhzD7S+G/YthdUTmX3iN9Iz71CJScOm7eS/U34huXh1bmxZjb5NKhHgox/LIiKFiX6qi4i4kjFQoyM2uAPfz+oDcQeymm47tZ8bk+/hlEdbvlrdmavmNqBH02rc2LK6GwMWEZErScW2iEgeMMYwudcXzNk9h5k7ZnI04Qh9qraE+O8pdexX7uJXhvuVYt2JnrzwaUt2lplDh/Nr+E+bG6lWooK7wxcRkX9IxbaISB4p6VuSW0Ju4eaGN7Pv3D6KB9XImNu94StY8xkexzbR4tCXPOD7NXf5luOX4zv5ZfZnVPJqxX87/JcmWslERKTAUbEtIpLHjDHUCKqRseFXElqNgJZ3wqE1sGYSM89Eg43L7JxOcPI+xn85lT2+IdzQshoDmlWmhL+3+wYgIiJOU7EtIpIfGANVWpBQvj7LZnSC1D+a7jy2jrALqzjvXYP5a7sycH4LmtSvx/Xh1QgLLo6Xh5f74hYRkb+kYltEJB8J8Argp2t+4u15b7PRbiQl5Twtmg+EDdPwi9tLv7jx9DUeHDzZjg9nNeOu0vNoWLwDo8JvpF3VpppmIiKSz6jYFhHJZ0r5luLq4lfzwlUvEJsUi/ErA52fhV0LYN0XmO0/Uu3kEmqX3op1+LE5/ifuXvQTpWxrHmv+Ip0blMPXy8PdwxAREVRsi4jkW8YYyviVydjw8IR6PTIeCSex66cxe980SDmb1X/o+fWcWfAm/WaFER7aiGtaVKGZLqoUEXErFdsiIgVNQBmOhl7DuYNfQUrGLg8M15zYSen07dzIJ8TsD2PqttY859WGbs1rU63yIXrX6YinQz/2RUTykn7qiogUQBWLVWTBtQtYErOEb3Z9gweG0i17wIbpmB0/U/VMFI8RxSOpHzFlS1OePniU51cEElbmau4Ju4FmFRq5ewgiIkWCim0RkQLKy8OLztU707l6Z9LS08DhAQ36wPkzsHUObJiBY98yNgWcASDVxLEy9lu2z9xGq3JP0adpFdrWKo2Xh8Ot4xARKcxUbIuIFAIejmwXRPqVgOY3Q/ObSYzdxS9zr4f0P5r/F7+EBnEbWbC3DXdeCKNyow70aVKFljVK4eHQ/G4RkStJxbaISCGW6F+CnrX6MG/fPOJT4qngXYIWPhbH2YMMYBYDmEX89vIs2N6KW6nP3irLubpaV24J7U+tkjXdHb6ISIGnYltEpBAr41eGF9q+wP+1+j8iD0aSnJaMo0ZvOLgqY6rJlm8pdu4QA5jDidIrWZfuy6x9E5m1byI1fCJ4uf1LNK4cpBVNRET+IRXbIiJFgI+HD92Du/+xo3qbjEe3V+DwGtgym/knF0PKmawuA44tZO+k43zq2YbKTbvQo3FVQioVV+EtIpILKrZFRIoyhwOqhHGiVHV2zPw+a7fB0O/cEcqkHaJ/6k8krh5D5OowJjtaERgSQWLxxQxp1ItGZUJUfIuI/AUV2yIiQln/svxy3S8sPLCQefvmYa2lTLcvYMu3sGU2/rG76MVCerGQBTs+4KFyxfnuwGT8TBm6Ve3P8x3vx1OrmoiI/ImKbRERATJuEz+47mAG1x1ManoqODyhYihc/TSc2AZbv4dt37HIKxG4AMB5e5ITW79m9JLD2Lo9ad88hHa1yuDtqcJbRARUbIuIyGVcdKdJY6BcAyjXgJQODxI5PQKSL2Q1335uO63sBuz28Wzf04D3UpqTWLMnYc1a4Be0m9aVwvD19M37QYiI5AMqtkVExHkWnmz1JAv3L2TZoWX4evjQousjsH0uZvci6qdsoT5bYM8X/HaoOndUsjjwoX7xltzQsDcD6/V29whERPKUim0REXGal4cXfWr2oU/NPiSmJLL33F48S4dAi1vgQhzsWpAx3WTnPLYEJAD+pHOBLeeW8u6izWyeH0P1Zp3p2qgKFYJ0tltECj8V2yIi8o/4e/kTUjrkjx0+gRAyMOOReoEFPwyBMzuzmocn7WfomceI/yWIBfObs7l4B8o17UGX0GDKFrcU8y7mhlGIiLiWim0REbni0j28aFSxJUeSz3I88TgAnetdAzt/odipPQxgEQPOLeLCsjeIXBrKwGqnKOVTna7VuzCscR+qFq/q5hGIiFwZKrZFROSKcxgHT7R8gsfCH2PTyU1sOLGBCg2HQQ+bbWWT7/E5sg5/v/WkOcpxImUHU3btYMb2T7kp6HXaNw+lRfWSeDi0jreIFFx5WmwbY3oA7wAewMfW2tGXtJvM9l5AInCrtXZNZts+IA5IA1KttWF5GLqIiPwDDuMgtGwooWVDM3ZkW9mEqx6FMwdZuPT/4NS6rGP6pibz8Mb+bN1Sj7E2nOQ6vWjRPJy2tUrj4+nhnoGIiPxDeVZsG2M8gLFAVyAGiDLGzLHWbsnWrSdQJ/PRCvgw8+vvOllrT+ZRyCIi4molqnLSPwhO/bGrR0AV8DxCg9TtNGA7bP+CA7uqMymtBbNqnKdG2XrcFNqX1pWb4zBaz1tE8re8PLPdEthlrd0DYIyZBvQHshfb/YHPrbUWWGmMKWGMqWitPZKHcYqISB56v/P7HE04ysIDC1l+aDnhnd6BtOQ/VjbZ8TPVLuynn0cM71OZQyc2s2zhN/jYIB5q8DFBydbdQxARyVFeFtuVgYPZtmO4+Kx1Tn0qA0cAC8wzxljgI2vteBfGKiIieahCQAWGNhjK0AZDM3Z4eEHD/hmP1GTYt5QF6z6EhO1Zx9SyyfT9uTs/pzbj9e0bqdy8B10aV6N8cS0pKCL5R14W25e7wuXS0xF/1aedtfawMaYcMN8Ys81au+RPb2LMCGAEQPny5YmMjMxqi4+Pv2hbcqZcOU+5cp5y5Tzl6lIezE67+L+MzokXKJ5+hsGORRC7iMT5rxL5c1NW+7RiZ0VL3ZKlaV2iIV7Gy00x5z/6XDlPuXKecvXX8rLYjgGyr+VUBTjsbB9r7e9fjxtjZpExLeVPxXbmGe/xAGFhYTYiIiKrLTIykuzbkjPlynnKlfOUK+cpV3/WMqUlyw4tY8H+BSw5tITeN8+ExLPs/fE9apzfhP+xjfRyrKBrygraO2qx7lwKX531IdivBU+1+T9aVq1GxnX4RZc+V85TrpynXP21vCy2o4A6xpgawCHgBuDGS/rMAUZlzuduBZy11h4xxgQADmttXObzbsCLeRi7iIi4mb+XP92Cu9EtuBspaSl4eXhBYBX2B99AjYgIOLUXtn3P6m1fk2hiAbDmAsfOR2EmDeFdzzZ4NOxL+xZNCa0chENLCopIHsizYttam2qMGQX8TMbSfxOttZuNMSMz28cBc8lY9m8XGUv/3ZZ5eHlgVuYZCU9girX2p7yKXURE8hcvj8tMDSlVA9rexwJzBnbMyNrdPd1By/TNtEzeDOs+ZtuGOnxkWpJWrzcV65fnqto1KeNfKu+CF5EiJU/X2bbWziWjoM6+b1y25xa49zLH7QGauDxAEREp8PrW6ovDOFh4YCEnzp+gS8RLEB8HW+fAroXUT9lJfXbCli/pdbYBz685T2mPBlxdrQv3hF9DGf+S7h6CiBQiuoOkiIgUKk3LNaVpuaY82epJNpzYQMPSDcHDG0Kvg+RE2P0LbPueYzt/4qB3AgCx6Zv5at9myiyZS1qlPtQJ60rHBhXx99Z/kyLy7+iniIiIFEoO46BpuaYX7/T2hwZ9oEEfFm6ZDFGvZzU1dhTjnrSFcHAh5w4F8fM3zThcoTOVw3rRvkElyhTzz9sBiEihoGJbRESKpGI+QdQuUZtdZ3YB0KVmL6jQFbZ+T/FTuxloIuFYJElzX+HWlQ3YHehLizIduaVJH1pXaVzkVzYREeeo2BYRkSKpX61+9KvVj71n97Jg/wK61+wFxSpDlxfgxDbY9j1s/R6vI+s4VCKepPRYlsdOZfkvUxl4tg31G9xOhxahVC2lM94ikjMV2yIiUqTVCKrBnaF3/rHDGCjXIOPR8VHW7JrLmeWPZzUHOrx55tR0vJZPZ+Py2iz0b4dngz6Eh7WmboVAnfEWkYuo2BYREfkLO1PO4GE8SLNpAESUDsXLqz7sWkjj1F00TtwF0Z9xYE1FHgqqx64KFelVqwvDGnenuG+gm6MXEXdTsS0iIvIXbmxwI71q9CIyJjJjukndwVD1qotWNmHHT1Q7f4TjJYPYf2EPH25ZzoebX2aQ7UyXlvfRql41vD0d7h6KiLiBim0REZG/UcK3BANqD2BA7QF/7My2sglpqcTuns+mFX9MN8GkMfLgp5TZP4mVNOZI+QhKtxhA6yaNCPDRf78iRYX+tYuIiPxbHp4scyRjsVm7GgRUpmKlinDwNzqwFo6thblvsemHWswtW58NlSvRt24netXpSKC3ppuIFFYqtkVERK6APjX7UK14NRYdWMQvB3+hU43e0PRuiD8BO3+G7T/C7l9olLKbL7zSWH9uI+tX/8zLUQ76mFYMbjaKxiGN8PLQdBORwkTFtoiIyBXg4fCgWblmNCvXjIdaPESqTc1oKFYWmg3LeKScx+6OJGrNy5ByLqPdpNP/yGya753O7m8qs79EG7zrdaFB656ULlnCbeMRkStDxbaIiMgVZozBy3j9ucHLj13la3P890Ib8PPwoXm1q2HfUmolH6LWmZmwaiYXVnnxnV99/lfZnxYV29GlRmvCK7agjF+ZPByJiPxbKrZFRETyULXi1fiwy4csP7ScXw//SuVilfHu8gGkpcDB32D3Qti1EJ8j6zjpE0NsWiDzYmYyL2YmTS+UZZR3K0o37EhwaAc8fYu5ezgi8jdUbIuIiOQhHw8f2lduT/vK7QG4kHYho8HDC4LbZTw6PwsJJ1mx4G44sy3r2F4Je2l1OBr2fUDKDx7s8qlLQvkwitfryMri8dQu15hGZRrh5+nnjqGJyGWo2BYREXEjHw+fy+5P8Q1iXfz+i/a1aP8onNgLB1bgeXQTtZO3wsGtnDv0Ba9WqwoGHHhQJ6AWn4c9gX/lxuDlmxfDEJEcqNgWERHJh7w8vJh/7XxWH13NmuNr2H5qO7XDRoLJWK3EJJ2Fg1FwYAXrDywCcwKAdNKIP30I/0ldSMWD4z7VOV+6IT5VmhJw3off9qdxDku9kvWoHFgZh9HqJyKupGJbREQknwryCaJz9c50rt75z42+QVCnC9Tpwto1QbBxQlZT69KVsRfS8YjdSaULe+DwHjj8PVWAO45X5zf/jPXAfT38ea3mtXSpfw2UqA4OjzwamUjRoWJbRESkgGtarin9avVj7fG1HIw7SP2GgzGDboDkBDi2BY5ugKMbOLtjOfv9ATLmiSelJVIx8r8w72WSjTfn/KuTUrIOPhUbcL5iNT6M20yt0g2pGVST2iVrU7lYZbeOU6QgUrEtIiJSwHWs0pGOVToCcCrpFJ6OzP/evQOganjGA1j+y1yOHfzjlvIexkGNyq1JObYD78RjlEnYCQk7IWYui/18+bZCOdj3EwD1PSvycWA4gVUb4ihbD0rXBr8SAKTbdE1HEcmBim0REZFCpJRvqRzb0mwawxoMY9upbWw/vZ3y/uXx7z8rozHpLJzcCSe2w8nt7DyyHOzxrGMrnz1N0M73Yc0fr5foVZKk4jV5pIwnh709qFWiJjWKBzOk7nVUDqruqiGKFCgqtkVERIqIQI9AHm+ZcWbbWsu55D9uroNvEFQJy3gA+5Y9Dbu/zWoOC+kODa8h5egWUo/vxOvMHvxTTuMfG80O/5qcTU7lUHwMS1hC3/lvcM6rImmlauNXoR6+FetDmXqcLV6BHUnHqVWi1l/+UiBSmKjYFhERKYKMMQT5BOXYfl2966hfqj67z+5mz5k91KvXHyqE4wV4AaSnw7lDJBzbyNmVf0xNcWAITk7GJ3kfJOyDgwuy2tb6+XFfhbIABHmXoHeZJjwZ3B/K1IGgauDQVBQpfFRsi4iIyJ+Elg0ltGxozh0cDihRlYPpCRftrlq8Gj5PrYLY3RC7E3tiBxeObif1+Ha22yNZ/c4mn+HQtl9h+WQAUh0+JAbWwJSti3/lEGb7gKNsPeqWqkvtErVzXI9cJL9TsS0iIiL/WP1S9Vl14yr2n9vPnrN7sFjw8oMKjaBCIwzw+2119i99EvZ8n3VsyyqNOOddHs9Tu/C/cJziZ7fB2W2waw7vVK3J6Z2pADhw8IVnfUKqhOMo3wDKNYASwToTLgWCim0RERH5V/y9/GlQugENSjf4y36VilWmfqn67D27lwtpF6jX+l6KV2yV0Zh0LuMCzZPbOXV4PadPzM06zlpL3V0/4dj5U9a+VIcviSVqs61sVSKr1adhxTAalm5IcPFgrYwi+YqKbREREckTo5qNYlSzUaSlp3Eo/hBl/cv+0ehbHKq0gCot2FuxHvz0R7Fdq3g1fPo+woUjm0k6tBmv2O0ZZ8JPbWJD2kEmJ2+GXV8D0C2tJM/Z0vhXa4pnxcZQoTGUqqkb9ojbqNgWERGRPOXh8KBa8Wo5tpf1K8tdoXex4/QOdpzeQd0yIdD8ZnyArJnb50/D8W3sWP8unN6UdWzDs0cofnY9HPwla1+Kw5fzJevzWc0Qguv3JLRsKNUCq2GMcc0ARbJRsS0iIiL5SrXi1RjVbFTWdlp62p87+ZWE6m3Yue71i3Y37/MWJKeQengjiQfW4XliE/5Jx0g+tZ7xgafg5FIA/PFlzmkPildpjl9wGFRqBuUagqcuxJQrS8W2iIiI5GsefzEFZFTTUWw8uZEtsVvYfno79YKvBi9/PEMGUvz3Tomn2Lp9Fmx4O+u4YtZQ/sxWOLMVNn0JQJrxJKFEfWIqN+BQYARp6Wl/+d4izlCxLSIiIgVW5+qd6Vy9M5BxIeVlp4b4l2KLt9dFu1rVuAo6vYM9vIbEfdHYw2sJiNtL8dObGO9xhoUJq3h98tuU967DfxJOE1E+FJ8araBKOARVBU1BESep2BYREZFC4a/mYDcr14zhjYez6eQmNp7cSGj5FlC9DaZ6GwLaZHa6EAdHNrDvt+fg/DHSucCR5E0ExB7HJ2YFRH8EwHnv0iRXbEGxmq3wCB0MJXVreslZnhbbxpgewDuAB/CxtXb0Je0ms70XkAjcaq1dk63dA1gNHLLW9smzwEVERKRAa1y2MY3LNgYy5oCn2cvMA/cJ5Fylxuw+fyxrl8HQpP8npB3ZzPm9K/E6Eo1fcix+++ex4/Aibt49i7L+zWhRvjkD/AxNvD0wVVtn3BVTZ7+FPCy2MwvlsUBXIAaIMsbMsdZuydatJ1An89EK+DDz6+8eALbCH9OwRERERHLDw+GBB5efix2XHEdElQhWH15NfHo89UrVI7B+H6jfh2KdAGvh1B6IWc1v22eScH4bCRcWs+/AYrYnejP12C4AkrxKcKFiOMXqtMOjepuMCzB18WWRlJdntlsCu6y1ewCMMdOA/kD2Yrs/8Lm11gIrjTEljDEVrbVHjDFVgN7AK8DDeRi3iIiIFBGVi1Xmvc7vsWjRImq3qM3Z5LMXdzAGSteC0rXYcDYa9m7LauoS3IrEYvVwxKzC98JJfA/MhwPzAfiyfGs82rxKj7qNKeGVDskJEFA6L4cmbmIy6to8eCNjrgV6WGuHZ27fBLSy1o7K1ud7YLS1dlnm9kLgcWvtamPMTOA1IBB4JKdpJMaYEcAIgPLly7eYNm1aVlt8fDzFihVzyfgKG+XKecqV85Qr5ylXzlOunKdcOc+ZXL1w6AVOpp7M2n6o/EPU9K0J1uKbdIygs1vxP70VnzNb6FXRgxSTjk0NpOKFUnx7YhmJnpU4V6IBSaVDOFuiERd8y7l6WC6hzxV06tQp2lobdrm2vDyzfbmJS5dW+pftY4zpAxy31kYbYyL+6k2steOB8QBhYWE2IuKP7pGRkWTflpwpV85TrpynXDlPuXKecuU85cp5zuSqVWorNp3cRPSxaNafWM+wTsPw9vD+U7+1x9eS8uPNABjPOJK90/A66U2ZlEOUOXEITiwAIMG3ImnV2hB43UcYzz+/Tn6lz9Vfy8tiOwaomm27CnDYyT7XAv2MMb0AX6C4MeYLa+0wF8YrIiIikiM/Tz/CK4QTXiH8L/tFHY26aLt19Qi8hiyFoxuw+38lcecSvGJWEpB0hLFHVzFh3EPUD2pO5xptuObIR5SoFoIjuD2UrQ8OhyuHJC6Ql8V2FFDHGFMDOATcANx4SZ85wKjM+dytgLPW2iPAk5kPMs9sP6JCW0RERAqCaoHVaFupLWuPr+V86vmM4tzTG6qEYaqEEdDufkhPg+NbWLr8CdLil7A5fQlbdr9LyNGjtN2cBGRedFmpFYH1O+Go2RHKNlDxXQDkWbFtrU01xowCfiZj6b+J1trNxpiRme3jgLlkLPu3i4yl/27Lq/hEREREXKFHjR70qNGDlLQUNsVuolpgtT93cnhwtmQ1tsQfzNplSafeVc+QtC8au285fknH8N3/M+z/mXRgXtOxBLfuR91ygThS4sG7mJYbzIfydJ1ta+1cMgrq7PvGZXtugXv/5jUigUgXhCciIiLiMl4eXjQr1yzH9lVHVmGzXc5Wr2Q9SrceBa3JWHLw9F7Yt5yknZH8dHQlr537hLRZy0hLqMtYFtAsZTtp1doRUK8TpkZHKFVTxXc+oDtIioiIiOQDbSq14e1Ob7Py8EpWHllJ64qt/2g0JqN4LlUT3+Y3Eb38GRJ3zQb/k+C/nMWJDtoeOwk7v814AIm+5bHBHQgIHwa1OrlnUKJiW0RERCQ/CPQOpHO1znSu1hmA1PTUy/az1vLroV8v2td5wHjwLIHds5jEHZF4HFiOf9Ix7LaZjNnrxakGZWhfuwztSpwh8NxOqNEB/Eq6fEyiYltEREQkX/J0XL5MO5xwmNMXTmdt+3n60axcc/DwwpStR0CrEZCeDse3sGvrLCYd+Iqy54+xKqoBK4/s5vnkn0jHQWLpEHzrXY1nrU5QrTV4+eXV0IoUFdsiIiIiBUjlYpVZdsMyVh9bza+HfyXdpuPl4XVxJ4cDKjTi19ho7AHL8eRd4NjF0Xo1OXu8FcWOR1MsdiP8uhF+fYdUhzfxNXtT/MZPcTg0z/tKUrEtIiIiUsD4e/nTsUpHOlbp+Jf9lh5aetF2h4bXEXTN0IzbxR9YwYUdi0ja8QuBZ7YyJ+Y477w6jw61y9E52Itue0bjU7dzxnzvksEuHE3hpmJbREREpBBKt+lcSL1w0b4OlTtkPPEOgNpd8KndBZ9ekBR3hHdn98HX82WO+jQjckMavQ5/Bzu+A+B8YHW863bGo/bVENwB/Erk8WgKLhXbIiIiIoWQwziY3GsyJ8+fZPmh5WyJ3UK14pdZ4xtYfW43F9KTuZCczObkRZwsUQbb5E1Ob1mAX8xy/OL2Q/REiJ5ImvHkwB0bCa5cEWNMxrKEkiMV2yIiIiKFWBm/MvSv3Z/+tfvn2GfZoWUXbXeoGoFnqzsp2epOSEuFw2s5v30BiVvns/1CIg99vh5PxzY61inDwzvuJ+VIY7zqdIbanTXl5BIqtkVERESKOA/jQaBXIHEpcQB0rJxtLriHJ1QNx69qOH5dnuSNbweSFvcM9Us1J+1CFfySD+K14wDs+AGA5OLBeNXriql1NdS8KmPKShGmYltERESkiHs0/FEebPEg646vY+mhpbSq2Oqy/Y7EH2HXmV0ArDmxAodx0CX8ba4qnc7ZjT8ReHg5fuf2QdQEiJrAyo6Tadi2J8V9vSAhNmNtb4cjD0fmfiq2RURERAQvhxfhFcIJrxCeY58lMUsu2m5StgnefsH4tIqgXKvbIS0Veyia0xt/4sDuxXywN4DoyIWEVA7ijaQXqHx+Bx51OmNqd4FaV0Oxsq4eltup2BYRERERpxhjqBBQgaMJRwG4qspVEJutg4cnplorSlVrxejFj7HxwCO0atuCKj7NufDbITyTYmHjjIwHkFahCR51ukCja6F8QzeMyPVUbIuIiIiIU66rdx2D6w5m15ldLI5ZTJfqXdgbu/dP/VLSU1h2eBmp6alEHVtFFKsYdOvnWONH7IYfOb9lHuVPrcb76Ho4up7lJ3yp1KU6NcoEwNmYjBVOSlR1wwivPBXbIiIiIuI0Ywx1StahTsk6AOzlz8X2uuPriEuOy9ou4VOCxmVDMQ4PynSpD10egpTznN+1lF2rZ7EwvR4/jF+Bv7cnr/pPpc3xaaSXqYujdleo0wWqtQUv3zwb45WkYltERERErqiElASCiwez79w+IONmOh4Oj4s7efnh16AbK1L28c2652jevDn1g1oRvymFROOH/8kdcHIHrByL9fTD1OgIoddB42vzfkD/goptEREREbmiIqpGEFE1gn1n97EkZgkNSjfIse/ig4tJs2lEHYsi6lgUNa96Dv9aH3Fu53KORn+H/4FFVLmwG3b+zLr44jhKdqFRpSAc52Ph8FoIbg9efnk4utxRsS0iIiIiLhEcFExwUHCO7SfPn2TjyY0X7etYpSN4eFG8fgTF60cAkHLmEPtXfUf02XJ8OX0dcUmpPFYuisGHRmM9fTHB7aFON2g5Aoxx4YhyT8W2iIiIiLjFscRj1Ayqye6zuwFoWLoh5fzL/amfV4nK1O4+ktrAHcDekwnsjdzFnuN1qJmyE3Yt4OTxIzga3UapAO+8HcTfULEtIiIiIm4RUjqE2QNmc/DcQZYcWkKQT5BTx9UoE0CNa0cCI4mLPcSeFXNYdzSZ9onJKrZFRERERLKrWrwqQ4sP/UfHBpauTJM+d9PkCsd0pRSt+2WKiIiIiOQhFdsiIiIiIi6iYltERERExEVUbIuIiIiIuIiKbRERERERF1GxLSIiIiLiIiq2RURERERcRMW2iIiIiIiLqNgWEREREXERFdsiIiIiIi6iYltERERExEVUbIuIiIiIuIiKbRERERERFzHWWnfH4DLGmBPA/my7ygAn3RROQaNcOU+5cp5y5TzlynnKlfOUK+cpV85TrqC6tbbs5RoKdbF9KWPMamttmLvjKAiUK+cpV85TrpynXDlPuXKecuU85cp5ytVf0zQSEREREREXUbEtIiIiIuIiRa3YHu/uAAoQ5cp5ypXzlCvnKVfOU66cp1w5T7lynnL1F4rUnG0RERERkbxU1M5si4iIiIjkmQJdbBtjehhjthtjdhljnrhMuzHGvJvZvsEY0/zvjjXGvJTZd50xZp4xplJejceVXJGrbO2PGGOsMaaMq8eRF1z0uXreGHMo83O1zhjTK6/G40qu+lwZY+7LbNtsjHk9L8biai76XE3P9pnaZ4xZl0fDcSkX5aqpMWZlZq5WG2Na5tV4XMlFuWpijFlhjNlojPnOGFM8r8bjSv8yVxONMceNMZsuOaaUMWa+MWZn5teSeTEWV3NRrgZn/kxPN8YUrZVLrLUF8gF4ALuBmoA3sB5oeEmfXsCPgAFaA6v+7ligeLbj7wfGuXus+TVXme1VgZ/JWM+8jLvHml9zBTwPPOLu8RWQXHUCFgA+mdvl3D3W/JqrS47/H/Csu8eaX3MFzAN6Zjs+0t1jzce5igKuynx+O/CSu8fqzlxltnUEmgObLjnmdeCJzOdPAP9191jzca4aAPWASCDM3ePMy0dBPrPdEthlrd1jrU0GpgH9L+nTH/jcZlgJlDDGVPyrY62157IdHwAUhkntLslVpreAxygceQLX5qqwcVWu7gZGW2svAFhrj+fFYFzMpZ8rY4wBrgOmunogecBVubLA72dog4DDrh5IHnBVruoBSzKfzweucfVA8sC/yRXW2iXAqcu8bn/gs8znnwEDXBF8HnNJrqy1W621210ce75UkIvtysDBbNsxmfuc6fOXxxpjXjHGHASGAs9ewZjdxSW5Msb0Aw5Za9df6YDdyGWfK2BU5p/bJhaSPzW6Kld1gQ7GmFXGmMXGmPArGrV7uPJzBdABOGat3XlFonUvV+XqQeCNzJ/tbwJPXrmQ3cZVudoE9Mt8PpiMv2AWdP8mV3+lvLX2CEDm13L/Ms78wFW5KrIKcrFtLrPv0rOrOfX5y2OttU9Za6sCXwKj/nGE+ccVz5Uxxh94isLxy0h2rvpcfQjUApoCR8j4k39B56pceQIlyfjT5KPAjMwztwWZy35eZRpC4TirDa7L1d3AQ5k/2x8CPvnHEeYfrsrV7cC9xphoIBBI/scR5h//JldFjXJ1hRXkYjuGi3/brsKf/yyYUx9njgWYQuH485krclULqAGsN8bsy9y/xhhT4YpGnvdc8rmy1h6z1qZZa9OBCWT8ma6gc9W/wRjgm8w/T/4GpAMF/eJbl/28MsZ4AoOA6VcwXndyVa5uAb7JfP4V+jf4Vz+vtllru1lrW5DxS9zuKxy3O/ybXP2VY79Pn8j8WhimvbkqV0VXXk0Ov9IPMs5+7SGj4Pt9An/IJX16c/EE/t/+7ligTrbj7wNmunus+TVXlxy/j8JxgaSrPlcVsx3/EDDN3WPNx7kaCbyY+bwuGX+qNO4eb37MVWZ7D2Cxu8eY33MFbAUiMp93BqLdPdZ8nKtymV8dwOfA7e4eqztzla09mD9f9PcGF18g+bq7x5pfc5WtLZIidoGk2wP4lx+IXsAOMn7rfipz30hgZOZzA4zNbN+Y/Zt7uWMz939Nxny1DcB3QGV3jzO/5uqS199HISi2Xfi5mpzZdwMwh2zFd0F+uChX3sAXmf8O1wBXu3uc+TVXmW2Tfn+NwvJw0eeqPRBNRuGwCmjh7nHm41w9kLl/BzCaAv7L7hXK1VQypgCmkHFW947M/aWBhcDOzK+l3D3OfJyrgZnbF4BjwM/uHmdePXQHSRERERERFynIc7ZFRERERPI1FdsiIiIiIi6iYltERERExEVUbIuIiIiIuIiKbRERERERF1GxLSIiIiLiIiq2RURERERcRMW2iEgRYYy5yxhjjTFXZds3KnNfl8zt/xljthhj3nNfpCIihYenuwMQEZE8E0rGXUwbAIuNMf7AHcAJYKMxpibQzlrb0I0xiogUKjqzLSJSdDQm41bK9TO37we+AtKBEsBioLoxZq0xJsAtEYqIFDIqtkVEio4GwAygvjEmCLge+BXYZK3dDnwGPGOtbWatTXBjnCIihYamkYiIFAHGmKpArLV2jzGmHPAY8B5Ql4ypJZBx5vvbzP4DgN5AOWAs4J9921o7L08HICJSQKnYFhEpGkKBjZnP44AeQEvgbWBN5v4QYDOAtXY2MNsYUxJ401p7R/ZtQMW2iIgTVGyLiBQNjfmj2H6DjLPcacaYxsAkY0wgkGKtTbzkuKfJOLOd07aIiPwFFdsiIkVDY+BrAGvt99n2NwS2AE2BTb/vNMYYYDTwo7V2zaXbeRW0iEhBZ6y17o5BRETyGWPM/cAtQBSwDvDOvm2tHee+6ERECg4V2yIiIiIiLqKl/0REREREXETFtoiIiIiIi6jYFhERERFxERXbIiIiIiIuomJbRERERMRFVGyLiIiIiLiIim0RERERERdRsS0iIiIi4iIqtkVEREREXOT/AXrnq4wCo1VvAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots(figsize=(12, 6))\n", + "lstyles = [\"-\", \"--\", \":\"]\n", + "lwidths = [1, 2, 4]\n", + "for idx, averaging_method in enumerate(omega_averaging_methods):\n", + " fref_out, ecc, mean_ano, eccMethod = measure_eccentricity(\n", + " fref_in=fref_in,\n", + " dataDict=dataDict,\n", + " method=\"ResidualAmplitude\", \n", + " return_ecc_method=True,\n", + " extra_kwargs={\"debug\": False, \"omega22_averaging_method\": averaging_method})\n", + " ax.plot(fref_out, ecc, ls=lstyles[idx], lw=lwidths[idx], label=f\"{averaging_method}\")\n", + "ax.legend()\n", + "ax.set_xlabel(r\"$Mf_{22}$\")\n", + "ax.set_ylabel(\"eccentricity($f_{22}$)\")\n", + "ax.grid()" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "1d59dcbb-51ee-422f-94ac-e2a6beabdc6d", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuEAAAF3CAYAAAALl5VRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABmR0lEQVR4nO3dd3hc1bn+/e8aaTTqXXKTu417xRVsY9NNIJSQmHZCSUIIgZDkJCfkTSUHEpIAvwAh4VAMhFADoUMwBgswuMq99ypZvU3XSOv9Y0ZC7paRNJJ8f65rLs3s2XvPsx8Jc2tp7b2NtRYREREREWk/jmgXICIiIiJyqlEIFxERERFpZwrhIiIiIiLtTCFcRERERKSdKYSLiIiIiLQzhXARERERkXYWG+0CoiE7O9v269cv2mVEhcfjISkpKdpldBnqZ+tSP1ufetq61M/WpX62LvWz9bVGTwsKCsqstTmHLj8lQ3i/fv1Yvnx5tMuIivz8fGbOnBntMroM9bN1qZ+tTz1tXepn61I/W5f62fpao6fGmN1HWq7pKCIiIiIi7UwhXERERESknSmEi4iIiIi0s1NyTviR1NXVsW/fPvx+f7RLaVNpaWls3Lgx2mV0GdHuZ3x8PHl5eTidzqjVICIiIi2nEB6xb98+UlJS6NevH8aYaJfTZmpra0lJSYl2GV1GNPtpraW8vJx9+/bRv3//qNQgIiIiJ0fTUSL8fj9ZWVldOoBL12KMISsrq8v/9UZERKQrUghvRgFcOhv9zIqIiHROCuEiIiIiIu1MIVy+lPz8fC6++OIWbfOXv/wFr9fbRhW1rVWrVvHuu+9GuwwRERHp5BTCu7j6+vpol3CYrhrCQ6FQO1cjIiIinZWujnIEd721ng2FNa26z+E9U/nNJSOOu95ll13G3r178fv93HHHHdTX17Nz507+9Kc/AfD0009TUFDAww8/zD//+U8eeughgsEgkydP5m9/+xsxMTEkJyfz4x//mPfff5/777+fjz76iLfeegufz8fEiROZO3cuxhiWLVvGt771LZKSkpg2bRrvvfce69ato76+njvvvJP8/HwCgQDf//73+e53v3vUmmtqarj88svZvHkzM2bM4G9/+xsOh4N58+bxm9/8hkAgwMCBA3nqqaeYO3cuhYWFzJo1i+zsbL73ve+xePFiHnjgAR588EEefPBBduzYwfbt27n++utZuHAhBQUF/PjHP8btdpOdnc3TTz9Njx492L59O9///vcpLS0lMTGRxx9/nKFDh3LDDTeQmprK8uXLOXDgAH/605+48sorj1r/n//8Z15++WUCgQCXX345d911F6+99hqPPPIIH3zwAQcOHOCss85i/vz5/PrXv8bn87Fw4UJ+/vOfs2rVKsrLy9m1axfZ2dk8+OCD3HLLLezZswcI/8Jx5pln8tvf/padO3dSVFTEli1beOCBB1i8eDHvvfcevXr14q233sLpdPK73/2u6Xt1xhln8H//93+a9y0iItIFaSS8g5k7dy4FBQUsX76chx56iCuuuIJ///vfTe+/9NJLzJkzh40bN/LSSy/x2WefsWrVKmJiYnjuuecA8Hg8jBw5kiVLljBt2jRuu+02li1bxrp16/D5fLz99tsA3HjjjTz66KMsWrSImJiYps948sknSUtLY9myZSxbtozHH3+cnTt3HrXmpUuXcv/997N27Vq2b9/Ov//9b8rKyrj77ruZP38+K1asYMKECTzwwAP84Ac/oGfPnixYsIAFCxYwY8YMPv30UwA+/fRTsrKy2L9/PwsXLmT69OnU1dVx++2388orr1BQUMBNN93EL37xCwBuvvlmHn74YQoKCrjvvvu49dZbm2oqKipi4cKFvP3229x5551HrX3evHls3bqVpUuXsmrVKgoKCvjkk0+4/PLL6d69O4888gjf+c53uOuuu+jTpw+/+93vmDNnDqtWrWLOnDkAFBQU8MYbb/D8889zxx138KMf/Yhly5bx6quv8u1vf7vps7Zv384777zDG2+8wXXXXcesWbNYu3YtCQkJvPPOOwBH/V6JiIjISfKUQ8Ad7SoOo5HwIziREeu28tBDD/Haa68BsHfvXnbu3MmAAQNYvHgxgwcPZvPmzZx55pk88sgjFBQUMHHiRAB8Ph+5ubkAxMTE8LWvfa1pnwsWLOBPf/oTXq+X8vJyxo4dy/Tp06mtreWMM84A4JprrmkKfPPmzWPNmjW88sorAFRXV7N169ajXot60qRJDBgwAICrr76ahQsXEh8fz4YNGzjzzDMBCAaDTJ069bBtu3fvjtvtpra2lr1793LNNdfwySef8Omnn3LFFVewefNm1q1bx3nnnQeEp9f06NEDt9vN559/zte//vWmfQUCgabnl112GQ6Hg+HDh1NcXHzUfs+bN4958+Yxbtw4ANxuN1u3bmXGjBk8/PDDjBw5kilTpnD11VcfdR9f/epXSUhIAGD+/Pls2LCh6b2amhpqa2sBmD17Nk6nk1GjRlFfX8+FF14IwKhRo9i1axdw8PeqoqKCESNGcMkllxz1s0VEROTYat66k4S9n+D86eZol3IQhfAOJD8/n/nz57No0SISExOZOXMmfr+fOXPm8PLLLzN06FAuv/xyjDFYa7n++uv5wx/+cNh+4uPjm0a2/X4/t956K8uXL6d37978/Oc/x+/3Y609ah3WWh5++GEuuOCCE6r70OkSjfWdd955vPDCC8fdfurUqTz11FMMGTKE6dOnM3fuXBYtWsT999/Pnj17GDFiBIsWLTpom5qaGtLT01m1atUR9+lyuQ46nqOx1vLzn//8iNNt9u/fj8PhoLi4mIaGBhyOI//hKCkpqel5Q0MDixYtagrlR6rJ4XDgdDqb+uZwOAiFQod9r37729/qGuAiIiJfUsnuzbgDmYyNdiGH0HSUDqS6upqMjAwSExPZtGkTixcvBuCKK67g9ddf54UXXmiaAnHOOefwyiuvUFJSAkBFRQW7d+8+bJ+NIS47Oxu3280bb7wBQEZGBikpKU2f8eKLLzZtc8EFF/D3v/+duro6ALZs2YLH4zlq3UuXLmXnzp00NDTw0ksvMW3aNKZMmcJnn33Gtm3bAPB6vWzZsgWAlJSUptFhgBkzZnDfffcxY8YMxo0bx4IFC3C5XKSlpTFkyBBKS0ubQnhdXR3r168nNTWV/v37869//QsIh+nVq1e3qN+Nxzp37lzc7vCfqfbv309JSQmhUIgbb7yR559/nmHDhvHAAw8csfZDnX/++fz1r39ten20XxKO5NDvVeNfIkREROTkpQf2U+PqGe0yDqMQ3oFceOGFhEIhRo8eza9+9SumTJkChAPz8OHD2b17N5MmTQJg+PDh3H333Zx//vmMHj2a8847j6KiosP2mZ6ezne+8x1GjRrFZZddxvjx45vee/LJJ7n55puZOnUq1lrS0tIA+Pa3v83w4cMZP348I0eO5Lvf/e4xr/wxdepU7rzzTkaOHEn//v25/PLLycnJ4emnn+bqq69m9OjRTJkyhU2bNgHhudyzZ89m1qxZAEyfPp29e/cyY8YMYmJi6N27N9OmTQMgLi6OV155hZ/97GeMGTOGsWPH8vnnnwPw3HPP8eSTTzJmzBhGjBjR9AtGS5x//vlcc801TJ06lVGjRnHllVdSW1vL73//e6ZPn8706dN54IEHeOKJJ9i4cSOzZs1iw4YNjB07lpdeeumw/T300EMsX76c0aNHM3z4cB599NETruXQ71XjVCMRERE5SaEAmQ3lBFP7RLuSw5hj/am+q5owYYJdvnz5Qcs2btzIsGHDolRR+6mtrSUlJQUIz39OTk4G4N5776WoqIgHH3wwmuV1Os37GS1d6Wc3Pz+fmTNnRruMLkU9bV3qZ+tSP1uX+nk4T+Emkh6bzIdDf8c5V93R4u1bo6fGmAJr7YRDl2tO+CnsnXfe4Q9/+AOhUIi+ffvy9NNPR7skERERkVZTuncLG+JdbEr2c060izmEQvgpbM6cOU1zzI9n7dq1/Nd//ddBy1wuF0uWLGmL0lpdZ69fREREWq62aBsvpiSz1vMh3+dn0S7nIArhckJGjRrVopMMO5rOXr+IiIi0XF3ZTvY6nfRN7xftUg6jEzNFREREpEuKqdnFnlgnfdM63omZCuEiIiIi0jX59+CJMfRJUQgXEREREWl71uI3ZQD0Te0b5WIOpxAuIiIiIl2Ou2wvB5zhS3H36YDXCVcI72Qar+t9qF//+tfMnz8fgL/85S94vd7j7mvmzJkcer30ruz3v/99q+5v165dvPzyy626TxEREWkdJbs3stsZi8GQl5wX7XIOoxDeSVhraWhoOOr7v/vd7zj33HOBEw/h0VBfXx+1z26LEP6vf/3riO8d6w6jIiIi0vZq929iT2wsua4cnDHOaJdzGF2i8EjeuxMOrG3dfXYfBbPvPeYqDzzwAHPnzgXCt46/7LLLmm7vvmjRIl5//XUA/vu//5sFCxaQkZHBiy++SE5ODjfccAMXX3wxhYWFFBYWMmvWLLKzs1mwYAHf+973WLZsGT6fj0suuYR77z12HY2ab3fllVdy11138d577/HUU081jQDn5+dz//3389ZbbzFv3jx+85vfEAgEGDhwIE899RTJycn069ePm266iXnz5nHbbbdRW1vLY489RjAYZNCgQTz77LMkJiayfft2rr32Wurr65k9ezYPPPAAbrcbgD//+c+8/PLLBAIBLr/8cu66666j1n3ZZZexd+9e/H4/d9xxBzfffDN33nknPp+PsWPHMmLECJ577rkjbvvPf/6Thx56iGAwyOTJk/nb3/7GihUr+Na3vsXSpUupr69n0qRJvPTSS9x5551s3LiRsWPHcv3115ORkcE777yD3+/H4/Hw1ltvcfvtt7N27VpCoRC//e1vufTSS3n66ad5/fXXqa+vZ926dfz3f/83wWCQZ599FpfLxbvvvktmZiaPP/74EfskIiIixxcq3c4up5P+GQOjXcoRaSS8gygoKOCpp55iyZIlLF68mMcff5zKyko2b97MN7/5TVauXEnfvn3xeDyMHz+eFStWcNZZZx0WRn/wgx/Qs2dPFixYwIIFCwC45557WL58OWvWrOGzzz5jzZo1J1RT8+0+/vhj1qxZw3nnncfixYvxeDwAvPTSS8yZM4eysjLuvvtu5s+fz4oVK5gwYQIPPPBA077i4+NZuHAhV111FVdccQXLli1j9erVDBs2jCeffBKAO+64gzvuuINly5bRs2fPpm3nzZvH1q1bWbp0KatWraKgoIBPPvnkqHXPnTuXgoICli9fzkMPPUR5eTn33nsvCQkJrFq16qgBfOPGjbz00kt89tlnrFq1ipiYGJ577jkmTpzIV7/6VX75y1/yP//zP1x33XWMHDmSe++9l6lTp7Jq1Sp+9KMfAbBo0SKeeeYZPvroI+655x7OPvtsli1bxoIFC/jpT3/a1Ld169bx/PPPs3TpUn7xi1+QmJjIypUrmTp1Kv/4xz8AjtonEREROb7Y6h3scTrp3wGvEQ4aCT+y44xYt4WFCxdy+eWXk5SUBIQD2Keffkrfvn2ZMmVK03oOh6PpLpfXXXcdV1xxxXH3/fLLL/PYY48RCoUoLCxkw4YNjB49ukXbFRUVNW134YUX8tZbb3HllVfyzjvv8Kc//YmPP/6YDRs2cOaZZwIQDAaZOnVq076a35lz3bp1/PKXv6Sqqgq3280FF1wAcNBo/zXXXMNPfvITIBzC582bx7hx4wBwu91s3bqVGTNmHLHuhx56iNdeew2AvXv3snXrVrKyso57vB9++CEFBQVMnDgRAJ/PR25uLhCecz9x4kTi4+N56KGHjrqP8847j8zMzKa633zzTe677z4A/H4/e/bsAWDWrFmkpKSQkpJCWloal1xyCRC+qVDjL0lH65OIiIgcX4N/Dz6H6ZBXRgGF8A7DWnvE5Y2h/GiMMcd8f+fOndx3330sW7aMjIwMrr32Wvx+/3HrOXS7G264oWm7OXPm8Mgjj5CZmcnEiRNJSUnBWst5553HCy+8cNzjuOGGG3j99dcZM2YMTz/9NPn5+cesxVrLz3/+c7773e8et+78/Hzmz5/PokWLSExMZObMmSd0vI2fc/311/OHP/zhsPcqKipwu93U1dXh9/uP+n1pvtxay6uvvsqQIUMOWmfJkiW4XK6m1w6Ho+m1w+Fomk/e0j6JiIhIhLUEHGVAVoe8RjhoOkqHMWPGDF5//XW8Xi8ej4fXXnuN6dOnH7ZeQ0MDr7zyCgDPP/8806ZNO2ydlJQUamtrAaipqSEpKYm0tDSKi4v54IMPTqieQ7d77733mt6bOXMmK1as4PHHH28a4Z4yZQqfffYZ27ZtA8Dr9bJly5Yj7ru2tpYePXpQV1d30NSQKVOm8OqrrwLw4osvNi2/4IILmDt3btP88P3791NSUnLEfVdXV5ORkUFiYiKbNm1i8eLFTe85nU7q6uqOesznnHMOr7zyStO+Kyoq2L17NwA333wz//u//8u1117Lz372MyDc58aajuSCCy7g4YcfbvoFa+XKlUdd90iO1icRERE5ttqyvRRFzsXUSLgc0/jx47nhhhuYNGkSED4xMyMj47D1kpKSWL9+PaeffjppaWm89NJLh61z8803M3v2bHr06MGCBQsYN24cI0aMYMCAAQdNbTmWMWPGHLRd4zQTgJiYGC6++GKefvppnnnmGQBycnJ4+umnufrqqwkEAgDcfffdnHbaaYft+3//93+ZPHkyffv2ZdSoUU2/MPzlL3/huuuu4/777+crX/kKaWlpAJx//vls3LixaXpLcnIy//znP5umijR34YUX8uijjzJ69GiGDBly0PHefPPNjB49mvHjxx8x1A4fPpy7776b888/n4aGBpxOJ4888ggff/wxsbGxXHPNNdTX13PGGWfw0UcfMX36dGJjYxkzZgw33HDDYd+vX/3qV/zwhz9k9OjRWGvp168fb7/99gn1/1h9EhERkWMr3rmePc5YYnDQM7nn8TeIAnO0aRBd2YQJE+yh18feuHEjw4YNi1JF7ae2tpaUlJRol3FEXq+XhIQEjDG8+OKLvPDCC7zxxhvRLuuYOkI/u9LPbn5+PjNnzox2GV2Ketq61M/WpX62LvXzCwX//gvPHXiY9VkDeX/OvJPeT2v01BhTYK2dcOhyjYRLh1FQUMBtt92GtZb09PSmyzWKiIiItERd6TZ2OJ0Myhoc7VKOSiFcmDx5ctMUkkbPPvsso0aNatc6pk+fzurVq09o3fLycs4555zDln/44YfHvRLKl9lWREREOj5n9VZ25zo5K2NQtEs5KoVwYcmSJdEuocWysrJYtWpVu28rIiIiHV9d3W5CxsnA9I55ox7Q1VFEREREpAsJ+r3UxlYDMCBtQJSrOTqFcBERERHpMgp3rGNnXHiyR/+0/lGu5ugUwkVERESky6jYtZYdcU5ynJkkOY9908NoUggXERERkS4jWLSR7U4nAzM77pVRQCFcDnHttdcyZMgQRo4cyU033dR0h8nnnnuO0aNHM3r0aM4444wTvoqJiIiISHuKrdrKTqeTwVmH3zCwI1EIl4Nce+21bNq0ibVr1+Lz+XjiiScA6N+/Px9//DFr1qzhV7/6FTfffHObfH59fX2b7FdERERODXWBnQQchoFpHffKKKBLFB7RH5f+kU0Vm1p1n0Mzh/KzST875joPPPBA0w1qvv3tb3PZZZdx4YUXMm3aNBYvXsyYMWO48cYb+c1vfkNJSQnPPfcckyZNwuPxcPvtt7N27VpCoRC//e1vufTSS/F6vdxwww1s2rSJYcOGsWvXLv70pz9x1lln8b3vfY9ly5bh8/m48sorueuuuwC46KKLmuqZNGkS+/btA+CMM85oWj5lypSm5Ufy6KOP8uijjwJQXV1Nv379WLBgAfPmzeM3v/kNgUCAgQMH8tRTT5GcnEy/fv246aabmDdvXtPNen7/+99jreUrX/kKf/zjHwGOuv2yZcu444478Hg8uFwuPvzwQxITE/nZz37G+++/jzGG73znO9x+++0n8Z0TERGRzqI+FMLrKAcyGZDeca+MAhoJ7zAKCgp46qmnWLJkCYsXL+bxxx+nsrKSbdu2cccdd7BmzRo2bdrE888/z8KFC7nvvvv4/e9/D8A999zD2WefzbJly1iwYAE//elP8Xg8/O1vfyMjI6Np9LqgoKDp8+655x6WL1/OmjVrmka4m6urq+PZZ5/lwgsvPKzWJ598ktmzZx/1WG655RZWrVrFsmXLyMvL48c//jFlZWXcfffdzJ8/nxUrVjBhwgQeeOCBpm3i4+NZuHAhM2bM4Gc/+xkfffRR0z5ef/31o24fDAaZM2cODz74IKtXr2b+/PkkJCTw2GOPsXPnTlauXMmaNWu49tprv+y3SERERDq44t2b2eMyQMe+PCFoJPyIjjdi3RYWLlzI5ZdfTlJS+CzeK664gk8//ZT+/fs33blyxIgRnHPOORhjGDVqFLt27QLCI8Rvvvkm9913HwB+v589e/awcOFC7rjjDgBGjhzJ6NGjmz7v5Zdf5rHHHiMUClFUVMSGDRsOev/WW29lxowZTJ8+/aA6FyxYwJNPPsnChQuPe0x33HEHZ599Npdccglvv/02GzZs4MwzzwQgGAwyderUpnXnzJkDwLJly5g5cyY5OTlAeHrMJ598Qmxs7BG337x5Mz169GDixIkApKamAjB//nxuueUWYmPDP+KZmZnHrVdEREQ6t+Kda9jmdJIem0qaKy3a5RyTQngHYa094nKXy9X03OFwNL12OByEQqGmbV999VWGDBlyQvvcuXMn9913H8uWLSMjI4MbbrgBv9/f9P5dd91FaWkp//d//3fQdmvWrOHb3/4277333nFv7/7000+ze/du/vrXvzbVct555/HCCy8ccf3GXz6OVvPRtl+zZg3GmCOuf6TlIiIi0nV5961jS1wcQzKHHH/lKNN0lA5ixowZvP7663i9XjweD6+99tpho9BHc8EFF/Dwww83BdiVK1cCMG3aNF5++WUANmzYwNq1awGoqakhKSmJtLQ0iouLee+995r29cQTT/D+++/zwgsv4HB88eOxZ88errjiCp599llOO+3YZxsXFBRw33338c9//rNpH1OmTOGzzz5j27ZtAHi9XrZs2XLYtpMnT+bjjz+mrKyM+vp6XnjhBc4666yjbj906FAKCwtZtmwZALW1tYRCIc4//3weffTRpl9UKioqTqiXIiIi0nk5yjewPc7J8JyR0S7luDQS3kGMHz+eG264gUmTJgHhEzMzMjJOaNtf/epX/PCHP2T06NFYa+nXrx9vv/02t956K9dffz2jR49m3LhxjB49mtTUVMaMGcO4ceMYMWIEAwYMaJriAeH53H379m2aKnLFFVfw61//mt/97neUl5dz6623AhAbG8vy5cuPWM9f//pXKioqmDVrFgATJkzgiSee4Omnn+bqq68mEAgAcPfddx8W6Hv06MEf/vAHZs2ahbWWiy66iEsvvRTgqNu/9NJL3H777fh8PhISEpg/fz7f/va32bJlC6NHj8bpdPKd73yH22677YT6KSIiIp2TP7iNOmM4LbNjX54QwBztz/9d2YQJE+yhAXLjxo0MGzYsShW1jfr6eurq6oiPj2f79u2cc845FBQUHHcqiZy42tpaUlJSolpDV/rZzc/PZ+bMmdEuo0tRT1uX+tm61M/Wdar30+/z8v5fB/PL3Ez+/dV/Mzjjy9+spzV6aowpsNZOOHS5RsK7MK/Xy6xZs6irq8Nay9///nfi4uKiXZaIiIhIq9u7ZRXb42KJJYZ+af2iXc5xKYR3YSkpKYdNGamtrW21/ZeXl3POOecctvzDDz/UaLuIiIi0q4qdK9gS56RPUh5OhzPa5RyXQrictKysLFatWhXtMkRERERoKFrHprg4zug2+vgrdwC6Okozp+L8eOnc9DMrIiISVl+zkfLYmE5xeUJQCG8SHx9PeXm5Qo10GtZaysvLiY+Pj3YpIiIiUWWtxWP3AHSaEK7pKBF5eXns27eP0tLSaJfSpvx+v0JbK4p2P+Pj48nLy4va54uIiHQEZcV7ORBXB8CQDIXwTsXpdNK/f/9ol9Hm8vPzGTduXLTL6DLUTxERkegr3LScja44MmNSyYg/sfusRJumo4iIiIhIp+bZXcD6uDhGdII7ZTZSCBcRERGRTq2+dDW74pyM6T4+2qWcMIVwEREREenUakNbARiRPSLKlZw4hXARERER6bQqSosoifMAMCJLIVxEREREpM3tXb+I9a44cmLTO81JmaAQLiIiIiKdmHfXcta74hiZ2znulNlIIVxEREREOq1Q2Ur2OjvXSZnQziHcGHOhMWazMWabMebOI7xvjDEPRd5fY4wZ3+y9ucaYEmPMukO2yTTGfGCM2Rr52nn+DiEiIiIiX0p1aDvQuU7KhHYM4caYGOARYDYwHLjaGDP8kNVmA4Mjj5uBvzd772ngwiPs+k7gQ2vtYODDyGsRERER6eIqSvZzwOUDYFjmsChX0zLtORI+Cdhmrd1hrQ0CLwKXHrLOpcA/bNhiIN0Y0wPAWvsJUHGE/V4KPBN5/gxwWVsULyIiIiIdy971i1jriqO7M4s0V1q0y2mR9rxtfS9gb7PX+4DJJ7BOL6DoGPvtZq0tArDWFhljco+0kjHmZsKj63Tr1o38/PwWFd9VuN3uU/bY24L62brUz9annrYu9bN1qZ+t61Tsp2fFe6xKdzEwpm+bHHtb9rQ9Q7g5wjJ7EuucFGvtY8BjABMmTLAzZ85sjd12Ovn5+Zyqx94W1M/WpX62PvW0damfrUv9bF2nYj/nL/01FbExfH/sV5g5ZGar778te9qe01H2Ab2bvc4DCk9inUMVN05ZiXwt+ZJ1ioiIiEgHVx8KUWl2AzA2d2x0izkJ7RnClwGDjTH9jTFxwFXAm4es8ybwzchVUqYA1Y1TTY7hTeD6yPPrgTdas2gRERER6Xj2bC5gU7whgTgGpg2Mdjkt1m4h3FobAm4D3gc2Ai9ba9cbY24xxtwSWe1dYAewDXgcuLVxe2PMC8AiYIgxZp8x5luRt+4FzjPGbAXOi7wWERERkS6sdMOnrHK5GJY+lBhHTLTLabH2nBOOtfZdwkG7+bJHmz23wPePsu3VR1leDpzTimWKiIiISAcX2LeIrRlObulzZrRLOSm6Y6aIiIiIdDoVwY1YYxjbbWy0SzkpCuEiIiIi0qlUlxVS6KrFAKOzR0e7nJOiEC4iIiIincquVfmscrnIi+tOclxytMs5KQrhIiIiItKp1G77hJXxLqb2mRbtUk6aQriIiIiIdCpVNUvxORxMyTsj2qWcNIVwEREREek0aqvKKXSWAnB6t9OjXM3JUwgXERERkU5j54r5FMS76B3Xjcz4zGiXc9IUwkVERESk06jZuoAV8S6m9pke7VK+FIVwEREREek0qqqX4nc4mJrXOW/S00ghXEREREQ6BXd1Bfu7wHxwUAgXERERkU5iR8EHFCSE54NnxGdEu5wvRSFcRERERDqF6s0fsdLV+eeDg0K4iIiIiHQSpTWL8TscTOs9I9qlfGkK4SIiIiLS4RUX7WF3fCUxGCb1mBTtcr40hXARERER6fB2LnmHzxMSGJY0iCRnUrTL+dIUwkVERESkw6vZOY9NrjjOHnxBtEtpFQrhIiIiItKh1dc3UFq/FoAz8qZFuZrWoRAuIiIiIh3a1nVLWJvQQIqJZ1jmsGiX0yoUwkVERESkQyte+S6LEuKZ3G0yDtM14mvXOAoRERER6bIqSz6iIiaGWYPOj3YprUYhXEREREQ6rANlFex37gfgjJ5nRLma1qMQLiIiIiId1sbP32JhUhxDE3qTnZAd7XJajUK4iIiIiHRYtdveZJ3LxXmnXRLtUlqVQriIiIiIdEhuf5DyhtUAzOp7bpSraV0K4SIiIiLSIa1Zkk9BEnSLSWVQ+qBol9OqFMJFREREpEMqX/s6i+NdnN3vXIwx0S6nVSmEi4iIiEiHEww1UOX9lIDDwdkDZ0e7nFanEC4iIiIiHc6y1avZmugmESendzs92uW0OoVwEREREelw9i/9N/mJCUzvPgmnwxntclqdQriIiIiIdCj+uno81e9TGRPD7KFfj3Y5bUIhXEREREQ6lCVrN7ErqZQEYjiz15nRLqdNKISLiIiISIdStOQlPkpK4Kzuk4mPjY92OW1CIVxEREREOgx/XT3e6v9QFRPDhV10KgoohIuIiIhIB/LZqo3NpqJMi3Y5bUYhXEREREQ6jL2LX2JBUgJndZ/UZaeigEK4iIiIiHQQ5e4Ade73I1NR5kS7nDalEC4iIiIiHcIHS9ewPaWMJBPLtLyuOxUFFMJFREREpIOoKPgnHyUlMDtvFq4YV7TLaVMK4SIiIiISdVuKazH2I3wOB5eMuC7a5bQ5hXARERERibqPF37MilQfPWNTGJc7LtrltDmFcBERERGJqlB9A8FNz7EkPp6vDr4cY0y0S2pzCuEiIiIiElX5Gw8QSijAGsMlXfyqKI1io12AiIiIiJzaVn76Bp+nGEYn5dEntU+0y2kXGgkXERERkajZX+XDVf4q2+LiuGzEf0W7nHajEC4iIiIiUfPG5+s4kL6TBGK4aNCl0S6n3SiEi4iIiEhUhOobqF7xDz5IiueivBkkOZOiXVK7UQgXERERkahYsKmExPiP8DscfH3Md6NdTrvSiZkiIiIiEhVLPnmXlakhhsbnMTxreLTLaVcaCRcRERGRdrerzENK+YtsccXx9ZHXnxLXBm9OIVxERERE2t2Ln6ylOG078Ti4aPAV0S6n3bU4hBtjkowxMW1RjIiIiIh0fbX+Onxrn2FecjwX9zqL5LjkaJfU7o4bwo0xDmPMNcaYd4wxJcAmoMgYs94Y82djzOC2L1NEREREuoqXl+0hMSWfgMPBtaf/INrlRMWJjIQvAAYCPwe6W2t7W2tzgenAYuBeY8x1bVijiIiIiHQR9Q2WTQv/xbw0mJoygEEZg6JdUlScyNVRzrXW1h260FpbAbwKvGqMcbZ6ZSIiIiLS5Xy0qYQejn9TEhvLb06/I9rlRM1xQ/iRAvjJrCMiIiIi8v6CfIrSq+gTm820PjOjXU7UnMyJmecZYx43xoyNvL651asSERERkS5nzb4qulU8w5p4F9eMvBGHOXUv1HcyN+u5FbgR+KUxJhMY26oViYiIiEiX9PT8VdjMLSSTzGXDr4l2OVF1Mr9+lFprq6y1PwHOBya2ck0iIiIi0sVsK3GTvPsJPkqKZ07/r5DkTIp2SVF1MiH8ncYn1to7gX+0XjkiIiIi0hXNXbCeYOZSYo2D6yb+KNrlRF2LQ7i19o1DXj/ceuWIiIiISFdTWOWjYcNTvJ8Sx+U9Z5CdkB3tkqLuZE7MPFcnZoqIiIjIiZr7yVZcmflYDDdMuTPa5XQIJ3Ni5vfRiZkiIiIicgJKawOUr/gHn/WN4cKcceSl5EW7pA7hRG5b/2Dka0JkkU7MFBEREZET8mj+VjLT3sPncHDTlF9Eu5wO40Smo5wT+bow8vWkT8w0xlxojNlsjNlmjDnsbxEm7KHI+2uMMeOPt60xZqwxZrExZpUxZrkxZtKJ1iMiIiIibae4xs++5S/wbjqcnTaE07KGRLukDuNEQvh/jDGLgO7GmJuAfcaY+MY3T/TETGNMDPAIMBsYDlxtjBl+yGqzgcGRx83A309g2z8Bd1lrxwK/jrwWERERkSj7+0db6JHxBm6Hg1un/S7a5XQoxw3hkWkn1wL1QH/gV8BaY8x6Y8xLLfisScA2a+0Oa20QeBG49JB1LgX+YcMWA+nGmB7H2dYCqZHnaUBhC2oSERERkTZQWOWjZMU/eS/NckH6CIZkHzr2emo7oRMzrbU7jDHnWmu3NC4zxiQDI1vwWb2Avc1e7wMmn8A6vY6z7Q+B940x9xH+peKMI3145CouNwN069aN/Pz8FpTedbjd7lP22NuC+tm61M/Wp562LvWzdamfrauj9fOZdT7yMt/BbwwTEi7pULWdqLbs6XFDuDHGREamtzRfbq11A4ubr3O8XR1h2aHbHG2dY237PeBH1tpXjTHfAJ4Ezj1sZWsfAx4DmDBhgp05c+Zxyu2a8vPzOVWPvS2on61L/Wx96mnrUj9bl/rZujpSP3eVeUhc8Cv+0w8uyhrLVedfG+2STkpb9vRE5oQvMMbcbozp03yhMSbOGHO2MeYZ4PoT2M8+oHez13kcPnXkaOsca9vrgX9Hnv+L8NQVEREREYmS+/6zntSs/1BnHHxv+t3RLqdDOpEQfiHh+eAvGGMKjTEbjDE7gK3A1cD/s9Y+fQL7WQYMNsb0N8bEAVcBbx6yzpvANyNXSZkCVFtri46zbSFwVuT52ZG6RERERCQKCnZXYrbN5T+pDi7vNok+6f2iXVKHdNzpKNZaP/A34G/GGCeQDfistVUt+SBrbcgYcxvwPhADzLXWrjfG3BJ5/1HgXeAiYBvgJXxToKNuG9n1d4AHjTGxgJ/IvG8RERERaV/WWu5/q4CM3HycxPH9s+6Ndkkd1gnfMTMSgv8ZGZk+KdbadwkH7ebLHm323BK+I+cJbRtZvhA4/WRrEhEREZHW8Z91B+hT9Sjv9Hby/f6Xkp2YE+2SOqwTmY7SqDuw3BjzcuTGOUc6WVJERERETkHBUAOPvfsZu3LWkYOTb079ebRL6tBOOIRba39J+CY6TwI3AFuNMb83xgxso9pEREREpJN46rOdjLCPsj7eye1jbiHRmRjtkjq0loyEN04XORB5hIAM4BVjjO5SKSIiInKKKqr28dZH77Msez+nxabw1dHfinZJHd4Jh3BjzA+MMQWEbwv/GTDKWvs9wvOxv9ZG9YmIiIhIB3fP2xsYmf4khc5Y7px2NzGOmGiX1OGd8ImZhK+KcoW1dnfzhdbaBmPMxa1bloiIiIh0Bp9vK8O75Z983j/ARalDmdj37GiX1Cm0ZDqK69AAboz5I4C1dmOrViUiIiIiHV4w1MA9byzD0f0/xGH47/MejnZJnUZLQvh5R1g2u7UKEREREZHO5cmFOxkR/BtLEmP53uBvkJvcI9oldRrHnY5ijPkecCswwBizptlbKYTnhouIiIjIKWZHqZvXPvwA03cbgxypXDP1zmiX1KmcyJzw54H3gD8Azbtba62taJOqRERERKTDamiw3PnqGkZlP848Zyxzp92D0+GMdlmdyonctr4aqAaubvtyRERERKSje37pHtJKnmF+nwBfSx/JxP7nRrukTue4c8KNMQsjX2uNMTXNHrXGmJq2L1FEREREOoqiah+PvfcZZT0+IZsYfnzBo9EuqVM6kZHwaZGvKW1fjoiIiIh0VNZa/r9/r2VS2l95Py6WB0//Kanx6dEuq1Nqyc16njHGpDd7nWGMmdsmVYmIiIhIh/PC0r0Ed7/MR5k1XJDUl7NHXhftkjqtllyicLS1tqrxhbW2EhjX6hWJiIiISIezs8zDA28vIdTzXRIw3HnB49EuqVNrSQh3GGMyGl8YYzJp2R03RURERKQTCtU38KOXVnFm1kNscMXwq1HfIztF1wT/MloSou8HPjfGvAJY4BvAPW1SlYiIiIh0GH9dsI3E8hf5uE8lFyX248LTb412SZ3eCYdwa+0/jDEFwCzAAFdYaze0WWUiIiIiEnUr91Ty/ILPyen/Idk4+cXF/4h2SV1Ci6aTWGvXA+vbqBYRERER6UCqvEFue24Fk3o8wsfOGJ6Y9GtSEzKjXVaXcCK3rV9orZ1mjKklPA2l6S3AWmtT26w6EREREYkKay0/+ddqhvM0H6f4+WbGWCYP+1q0y+oyTug64cYYA4yw1u5ph5pEREREJMqe+HQnO7Z9TKD/CoYSzw9mPxbtkrqUE7o6irXWAq+1cS0iIiIi0gEU7K7kof8UkN37GeqN4b7zH8PlTIh2WV1KSy5RuNgYM7HNKhERERGRqKvwBLn9+RXM6PEQG13w2+E30bfH6dEuq8tpyYmZs4DvGmN2Ax6+mBM+uk0qExEREZF2VVffwK3PFXAa/+ST1CrmJA/mwsk/jnZZXVJLQvjsNqtCRERERKLunnc2cmDfp9T1W8owG89Pv/rPaJfUZbVkOsqt1trdzR+ArtQuIiIi0gW8vGwvry5aRWqfZ8EY7j//MVzOxGiX1WW1JISfd4RlGh0XERER6eRW7KnkV6+vZkLfv7A71nLf2B/Su6fmgbelE7lO+PcIj3gPMMasaVwMJAOft2FtIiIiItLGiqp93PJsAbO6PcJnCX5+kjudqeO+He2yurwTmRP+PPAe8AfgzmbLa621FW1SlYiIiIi0uVp/HTc+tYzBsa/wWdo+Lnbm8s0LHol2WaeEE7lZTzVQbYy5EbgC6Ne4nTEGa+3v2rRCEREREWl14SuhrKCu+lO29/mUEQ1x/ObyVzGOlsxWlpPVkqujvA5UAwVAoE2qEREREZE2Z63ll6+tY+POFaQOeIFka3jwK88Qn5Ae7dJOGS0J4XnW2gvbrBIRERERaRePLNjGmyvWctqARyk2lmem3kO3brr1S3tqyd8bPjfGjGqzSkRERESkzf1r+V4emLee8f3/wu7Yeh4Y9h2GDL0s2mWdcloyEj4NuNEYs4PwdBTdMVNERESkE/nPugP87NXVzOr3IEvj/NzV/RzOmPLDaJd1StIdM0VEREROAQu3lvGDF1Yyq/cTLE0o4zuJg7jiggejXdYpqyXTUfYA04HrI3fLtEC3NqlKRERERFrNyj2V3PzscqZ2f46lSdv5miOT26/4V7TLOqW1JIT/DZgKXB15XQvoQpIiIiIiHdjmA7Xc8NQyxmX+mxWpqznfJvKrb7yLiWnJhAhpbS3p/mRr7XhjzEoAa22lMSaujeoSERERkS9pS3Et1zy+mEHJ77Iu/XPObIjj3jnvEeNKinZpp7yWjITXGWNiCE9DwRiTAzS0SVUiIiIi8qVsKa7l6scW0ythATuy5jOm3sEDV7yJMzEz2qUJLQvhDwGvAbnGmHuAhcDv26QqERERETlpmw+EA3iPhE/Zn/0mg+vhr1/9F4lpedEuTSJOeDqKtfY5Y0wBcA7hyxNeZq3d2GaViYiIiEiLbT4QnoLSI+ETCrNfY1A9PH7RP0nNGRLt0qSZFs3It9ZuAja1US0iIiIi8iWs21/NN+cupechATytx9holyaHOOHpKMaYZ4wx6c1eZxhj5rZJVSIiIiLSIpsr6rn6scXkJS78IoDPflYBvINqyZzw0dbaqsYX1tpKYFyrVyQiIiIiLfLRpmLuW+5ncNbH7M16lUEhGx4B76mo1lG1JIQ7jDEZjS+MMZm0cDqLiIiIiLSu11fu5+Z/FDAy5x22p77J6JDlia88rxHwDq4lIfp+4HNjzCuEL1P4DXR1FBEREZGoefqznfz2rQ2c2f911sQvZlqd4YFLXyEhZ2i0S5PjaMnVUf5hjFkOnE346ihXWGs3tFllIiIiInJEDQ2W37+7kScW7mD6oOdZ5VzLeX7441Xv4kzrHe3y5ASccAg3xjwD3GGt/WvkdYYxZq619qY2q05EREREDuIL1vOjl1bxn/X7OWvoU6ww25hTn8CM/r9QAO9EdGKmiIiISCdR5g5w1eOLeX/TTqYM+wsrzDZuIZ1fXPMhDXFp0S5PWqAlc8IdxpiMSPjWiZkiIiIi7WhbiZsbn15Kha+QEac9xGbr5X9dA7js6/+CWFe0y5MWaumJmYuMMf/iixMz72mTqkRERESkyYcbi/nhi6tISdhBTt//o7yhjr/nnMWUi/4KjpZMbJCOQidmioiIiHRQ1loeWbCN+z/Ywui8NRQlPU9yfYjHh9zE4Gk/jXZ58iW05MRMA4wHMq21vzPG9DHGTLLWLm278kREREROTZ5AiJ++spp31xYxfeh/WE0+w+rq+esZd5Mz4mvRLk++pJZMR/kb0EB4JPx3QC3wKjCxDeoSEREROWXtKfdy87PL2VJSxrSRz7CqfgdfCcBvL/oH8XmKXl1BS0L4ZGvteGPMSghfHcUYE9dGdYmIiIickuatP8BP/rUaR2wZI4Y+wtpQDT+xaXzz6n9jUrpFuzxpJS0J4XXGmBjCJ2VijMkhPDIuIiIiIl9SXX0Df3xvE08s3MmI3juoSnqC8lCIv6eP44yvzoVYjX12JS0J4Q8BrwG5xph7gCuBX7ZJVSIiIiKnkP1VPm57fgUr91Qwc+SHrAh9yMC6Oh4c/m36nPmTaJcnbaAlV0d5zhhTAJxD+Oool1lrN7ZZZSIiIiKngI82FfPjl1cTooYpo5+koG4fl/hD/PK8R0gceE60y5M20qKb7VhrNwGb2qgWERERkVOGv66ee9/bxNOf7+K0vL0EUx5jWyDAXWRx+TUvYVJ7RLtEaUO646WIiIhIO9tQWMMPX1rJluIaZo36hBV1/6FPsI7/y5vNkPP/DDGKaF2dvsMiIiIi7aShwfLEwh3c9/4WUlNqmTTmKZYH9zHbX8dvzvojScMujXaJ0k4UwkVERETaQVG1j/9+eTWfby9n0tBN7DXPssdfxz0NaVwy5wVMRt9olyjtSCFcREREpA1Za3l5+V7ufnsj9cbLWeNeYYV/DWP9Af7Q/2vknfM7iHFGu0xpZ472/DBjzIXGmM3GmG3GmDuP8L4xxjwUeX+NMWb8iWxrjLk98t56Y8yf2uNYRERERI5nX6WXb85dys9eXUu/XvvpMfBe1vhWc7sPnjp/Lnnn/0EB/BTVbiPhkRv9PAKcB+wDlhlj3rTWbmi22mxgcOQxGfg7MPlY2xpjZgGXAqOttQFjTG57HZOIiIjIkTQ0WJ5bspt739sExs85Ez5gqecT+gbqeDZtPCOv+jvEp0W7TImi9pyOMgnYZq3dAWCMeZFweG4ewi8F/mGttcBiY0y6MaYH0O8Y234PuNdaGwCw1pa00/GIiIiIHGZXmYefvbqGJTsrGDu4CLfrCZa53XzTE+S2M39Lwthrol2idADtGcJ7AXubvd5HeLT7eOv0Os62pwHTI3fx9AM/sdYua8W6RURERI4rEKrn0fwdPJK/DVdcgLMnvMcyz+f099Xxj4RBjP2vxyC1Z7TLlA6iPUO4OcIye4LrHGvbWCADmAJMBF42xgyIjKZ/sWNjbgZuBujWrRv5+fknXnkX4na7T9ljbwvqZ+tSP1ufetq61M/W1ZX6uaG8nn+sD3DAaxmetxFP8osUuP3cWOvnnOyrqOx+DvkrtgBb2qyGrtTPjqIte9qeIXwf0LvZ6zyg8ATXiTvGtvuAf0dC91JjTAOQDZQ237G19jHgMYAJEybYmTNnfplj6bTy8/M5VY+9LaifrUv9bH3qaetSP1tXV+hnSa2fu9/eyJurC8nL8TNjxL9Z6V7DoGCQvyUOY8QNj0JK93appSv0s6Npy562ZwhfBgw2xvQH9gNXAYdOinoTuC0y53syUG2tLTLGlB5j29eBs4F8Y8xphAN7WVsfjIiIiJy6QvUNPLdkD/fN20ygro5zJq5mjftfbKqp40eeEP81/bc4R18F5kh/zBdpxxBurQ0ZY24D3gdigLnW2vXGmFsi7z8KvAtcBGwDvMCNx9o2suu5wFxjzDogCFx/6FQUERERkdby6dZS/vftDWwpdjN2cCX1SU+x1F3CTI+Xn/c4m55z/giJmdEuUzq4dr1Zj7X2XcJBu/myR5s9t8D3T3TbyPIgcF3rVioiIiJysJ1lHu55ZwPzN5aQlx3i3Cn/YUn1Qnp4QjwYSubsix+HPlOiXaZ0ErpjpoiIiMgxVPvqePjDrTyzaBcuZwMXTFnDqpp/UVBVx421Pm4ZcwuJZ/5QN92RFlEIFxERETmCYKiBF5ft4S/zt1LpDTBj9AGK7D/4vLqSmR4vP0kfR9/r7oeMvtEuVTohhXARERGRZhoaLG+tKeT+eVvYU+FlzAAPQ7JeZoV7M4OCQf6vPpUzLnwYBsyMdqnSiSmEi4iIiADWWj7ZWsaf/rOJ9YU1DO4R4oIzP2JRxUekVdfzS3cdX5v8E2InfAtiFKHky9FPkIiIiJzyVu2t4o/vbWLRjnJ6ZjYwe9pSllW8xbLyENfVePjugMtIu+pXuuqJtBqFcBERETllrdtfzYMfbuWDDcVkJsNXzlzPyqqX+azMzyVuD7dmjKXXVb+H3GHRLlW6GIVwEREROeU0D98p8YbZU3azxfscn1TUMMPr4wfOXgy5+P+g35nRLlW6KIVwEREROWUcGr6/MnU/O/wvsrC6lNH+AH+sT2LCzAdg2CW626W0KYVwERER6fKOFL63+1/kk6pShgaCPOg3zJp6J2b8N3XSpbQL/ZSJiIhIl2St5fPt5Tz68XY+3VoWCd/72O5/6Yvw7bPMOv37mEnfhrikaJcspxCFcBEREelS6hss768/wN/zt7N2fzXZKTF8ddo+tnjC4XtYIMjPfA3MPP37mEnfUfiWqFAIFxERkS7BX1fPv1fs5/FPd7CzzEPfbAeXnbWNtVWvsKC8Khy+vfXMPP1WzKSbwZUc7ZLlFKYQLiIiIp1ahSfIC0v38PTnuyitDTA8z8GlM9ewrOw1PizxMdHn5y5fA2eOvwUz5RZwpUS7ZBGFcBEREemcNh+o5anPdvLayv0EQg1MGmyZfPoilpS8y74DIc72erkp5GL0pB/D+G9q2ol0KArhIiIi0mk0NFg+2lTCU5/v5LNt5cQ7Hcwa48ekLGBRcT5bD1gucbu5wZHJgKl3wqivQ2xctMsWOYxCuIiIiHR4tf46XinYx9Of72J3uZfuaU6umF7K/rp3+Kx6E8luy3W1tVwX34fuZ98DQy8GhyPaZYsclUK4iIiIdFjr9lfz3JI9vLFqP95gPWP7Opk4ZjMrK97gg7IK+taF+Hl1DZfmTibp4h9A/7N0kx3pFBTCRUREpEPxBet5a3Uhzy3dw+q9VeEpJ6PqcaZ/xucH/sP2ohBTfD5+4w4ybcgVOC69FXJOi3bZIi2iEC4iIiIdwpbiWp5fsodXV+yj1h9iUG4c3zirmH2hD1hYsQ6XFy6udXNtfTyDT78VTr8BEjOjXbbISVEIFxERkajxBkO8u/YALy/by9JdFcTFOJgxwpKWu4IlJe/yXombvqF6flJdw6XJA0mf+RsYfinEOKNdusiXohAuIiIi7cpay7JdlTy5NsD3P5qPJ1hPv2wXXz+rnJKGj1hSVkDsPjjb4+Eb3iCTBl6MueDb0Ot0zfeWLkMhXERERNpFYZWPf6/YxysF+9hV7iU+Bs4ZE0NS9mqWlrzHf0qq6FnfwA+qa7g8Npvs038MY6/RlBPpkhTCRUREpM346+qZt6GYfy3fy8JtZVgLEwfEM3XsftaWvsknvt3E7IXpHi9fd3s5s8/ZxJz17fBVTnSJQenCFMJFRESkVYXqG/h8ezmvr9rPvPXFuAMheqa7+NqZfoLxi1lU9BGbCoMMqAvx3zW1XOzIIHvsd8J3tUztGe3yRdqFQriIiIh8adZaVu6t4s1Vhby9ppAyd5CU+FhmjnCQkr2JFeX/4f3yAyRbuKS2lsu8QXJTxtP9sp9Avxka9ZZTjkK4iIiInLRtJbW8saqQN1YVsqfCS1ysgxlDXXTvuYftnk/4uGIdZi9M8gf4QU0t56SdRvyUO2Dk18hfspruA2ZG+xBEokIhXERERFpkR6mb99Yd4J01RWwoqsFhYMqgJGaML6Uo9BnLipfQsNcypK6eH9fWMDsUR/dRc2DstdB9ZLTLF+kQFMJFRETkuLaV1PLu2gO8u7aITQdqARjTO4lrZ9VQE7OURUWfsHZ/kF718K3aai7yhRg08HyY9g0YdC7ExkX5CEQ6FoVwEREROYy1li3Fbt5ZW8R7a4vYWuLGGDi9bwr/NcuLx7mCxQfyefOAm0xruKKmhos8Xsb0nIKZMQeGXQLxqdE+DJEOSyFcREREAGhosKwrrGbe+mLeXVfEjlIPxsCEfil882wP7thw8N58wEMKDma5a5nt9jA5/TScE26GkV+D1B7RPgyRTkEhXERE5BTmr6tn0Y5y5m8oZv7GYoprAjgMTB6YyvTRbqpjlrOk6FM2FYWD99m1tZzv8TAlsRdxw2+E0d+AnCHRPgyRTkchXERE5BRT4QmyYFMJ8zcW88mWUjzBehLjYjhzcDJ5vUqoZAWLij5hXaGPNBycFwnekxN74RxxEwy/DLqN0C3kRb4EhXAREZFTwM4yD/M3FPPBxmKW76qgwUJuiovzx8STnr2Fff7lLC9exuK9ITKsg9nuGs73eJmYmKfgLdIGFMJFRES6oMZpJh9vLiV/cwm7yr0ADOmezFXTHDhTNrGh6nPmV24CN/RtcHBdTRWzvD7GJPcmZsS3FLxF2pBCuIiISBexs8xD/uYS8jeXsnhHOYFQA65YB5MHpjJzbCWBuLUsL1nIW2UHMGUwJgQ/qqlkptfPgO7j4fQbYehXIHtwtA9FpMtTCBcREemkfMF6Fu8oDwfvLaXsjox2989O4tIJiaRmbqMouIplxUtYWegjHgdT/QFura1hRqCerP4zYexFMGQ2JOdG9VhETjUK4SIiIp1EQ4NlQ1ENC7eV8dm2MpbsrCAYaiDe6WDygDTOHVdDyLWR1eWLea96O1RDT2L5am010zxeJhNPwmkXwvSLYODZ4EqO9iGJnLIUwkVERDqwPeXeptD9+fYyKr11AAzOTebyCYmkZu2gKLiCZQeWsmK/l1gME+osl9dUMt3no3/6IMxp18Lg86HPVIjR//pFOgL9lygiItKBVHiCfL49HLoXbitjb4UPgG6pLmYMSaF79/34Yjaxpmw570ZGu3sQy8XuGqZ5PEwOGRL7nQVjzoPB50F6nygfkYgciUK4iIhIFFX76li+q4IlOyv4bFsZ6wtrAEhxxTJpQBoXnu6jIX4LW6pX8HHZWup31ROPg9NDlsurK5nm8zEgOQ9z2pXh0N13Gjjjo3xUInI8CuEiIiLtqMobZOnOcOhevKOcDUU1WAtxMQ7G9EnjxpkunCnb2etdTUHJcpbu8+HAMNLGclN1JVN9PsaEHMT1OxOmngOnXQBZA6N9WCLSQgrhIiIibajcHTgodG8ursVacMU6GNsnjetnJJGQuovSug0sL17KK8XlUAz9iONSdw1TPLVM9AdI7T4WRl4GA2ZB70kQ64r2oYnIl6AQLiIi0oqKqn0s31XJ0kjo3lriBiDe6eD0vhncNDSJ+NRdFAfXs6KkgA0lpVACWcbJZH+QKTUVTPH56ZGSBwMvCofu/jMgMTPKRyYirUkhXERE5CSF6hvYdKCWgt2VLN9dScGuCgqr/QAkxcUwvl8GZ420uJJ3ciASutcUl0Ex5Jg4JgSCTKipYKLfT7+YZEz/6TBuFgycBZkDonx0ItKWFMJFREROkC9k+WRLaThw765g1Z4qPMF6IHz1ktP7pHHJJEts4i6KAhtZUVzAqqJyAHJNHJMCdUysKWeCL0Df2GRMv2kwahr0mwa5I8DhiObhiUg7UggXERE5Amst+yp9kVHuCgp2V7GpyItlKQ4DQ7uncsm4TLKzigg6d7C9Zh3Ly9by6b7wJQW7GxdTA0EmVpczwR+gt0K3iDSjEC4iIkL4qiWr91Wzak8Vq/dVsXpvFeWeIADJrljG9knj/EG1DBkOVQ1bWF+xmncrt2JrLA4MQ0wCl3k9jKutYJw/QHdncvhygSMjobvbSIVuEWmiEC4iIqecQKiejUW1rNpTGQ7ee6vYWeYBwBgYlJPMzCGZ9MitgvidFPk3sap0JSXeEj7fCokmljHWybk1XsZ6ahgdCJKU0jN8R8rxU6D3ZOg2Ahwx0T1QEemwFMJFRKRLs9ayq9zLqr2VrN5bzcq9VWwsrCFY3wBAboqLMXlpXDDGSULyfmrZzubKDXxcsZHA7gAA3R3xnF5nGVdZybiAn8HBEDHdRsCAS8PBu/dkSO8dzcMUkU5GIVxERLqMhgbLznIP6/ZXs25/NWv3V7O+sIZafwiAxLgYRvVK49ozsklLP0Bd7G521W5kXfk6FhdWAOAyMQw38czx+BlVW85Yf4DuJg56nc7ujLH0nX4V5E2AhPQoHqmIdHYK4SIi0inVN1h2lrlZu7+atftqWFdYzYbCGtyBcOCOi3UwrEcqF4/OJSerHBO/h5LgVtaVreOVsl1QBgboH5PE9EAdo6urGOX3MShYhzO9D/SaCGMmQt5E6DEGYuPYmZ9P38Ezo3nYItJFKISLiEiHF6pvYHupp2l0e93+ajYU1eCNXB4w3hkO3JeOyyUnqwLj2k953Q42VWzgP1VbCdWGg3l2TDyj6mP5am2AUe4qRgSCpDiToNd4GP3VcODOmwDJudE8XBE5BSiEi4hIh1Ljr2NTUS0bi2qaHpuLa/HXhedwJzhjGNEzla+d3p3szHJM/H7KgtvZVLGRd6u2EqoOB+5Uh4thuPgvXwMjqkoZHQjSvb4BkzsM+lwAvSaEQ3fOEJ1AKSLtTiFcRESioqHBsrfSy8aiGjY0C937Kn1N62QkOhnWI5WrJvUgK7McE7ePkuB2NlZs5O3KbYSqvgjcw0083/RZhleVMjwYJC9Uj0nrHZ5KMvS/woG75ziIT43WIYuINFEIFxGRNucNhth0oPnodi2bimqa7jbpMNA/O4mxvdO59PRUkpKLqYvdT6F3O1sqt/B65U7qK8LrpkUC9/V+GF5VxvBggF6hekxKz3DIHnIt9BwLPcZCck70DlpE5BgUwkVEpNUEQw3sLPOwubiWrcW1bD5Qy5biWnZXeLE2vE6KK5ZhPVK5YnwPcrNqiI0/QE3DbnbUbGVdxWY+3l/etL/usUkMsU5mBRwMqypheMBPz1A9Jrkb9BwPp80JB+8eYyGlW3QOWkTkJCiEi4hIi4XqG9hd4Y0EbTdbSmrZcqCWnWUeQg3htB3jMPTLSmR4z1QuGpNBSmox9bGFlAR2srlyM+9WbiNYHb4jpdPEMCg2hWkhyxB3PUNqyxgSrCOtoQFSe4XvNjnoyi8Cd2qPKB69iMiXpxAuIiJH1dBg2V/lY0txLZuLw0F7S7GbbaVugqHwiZLGQJ/MRAbnpnD28HTSU6twuIqpbdjHzprtbKvcxqdFhVAU3mdmTAKnORK4pi6W06rLGeJ107+uDqcjFrKHQI/pMG4UdI88EjOj2AERkbahEC4iIgRDDeyp8LCtxM320savbraVuJsuAwjQMy2ewd1SOGNQOpnpNcTEH8Br97O7dgfbqraxrGQvDcXhcB5rYugfm8KYBsPXA7EMrSpiSCBAdn0DxpUaDthDZ34RtnOGQqwrSh0QEWlfCuEiIqeQGn8d2w8J2ttL3Oyu8FIfmUYC4bA9MDeZr5/ei9wsN874EvymkL3ucNh+tWIXobLwlUkcGPo40zjNxjK7zsWg6mIG+9z0rgvhBEjvC91GwOBvQPeR4cCd3jc8hC4icopSCBcR6WKstRyo8bO9xMMHu+v48PV1TYG7pDbQtJ4zxtAvK4kh3VO4YGQ2GWk1xLjK8FPEPs8utldt553qHQSqvtimlzOVwcbFWfWJDKouYbCnkn51dbgskJQLucNg2NnQbTjkDg+PbruSo9AFEZGOTSFcRKSTqvbWsaPMza5yDzvLvOws87CrzMOOUnfTpf8AUuL3Myg3mRmn5dAn20FiUgU4i6mp38/uml3sqN7B4pJ9hIpDTdt0d6YwwMQzpyGRQW4/g6uLGVBXR6K1EJccDtv9Z4eDdmPgTsqORhtERDolhXARkQ7MEwiFw3W5h52lHnaWe5rCdqW3rmk9h4G8jET6ZSdxZZ88umfVEZ9Qxpbdn5PS08mOmh2srN7J+3tKmraJNQ76ONMYRBzn2ST615QzoLaM/o1hO8YFWYOg+xQYPQxyR4TDd3ofTSUREfmS2jWEG2MuBB4EYoAnrLX3HvK+ibx/EeAFbrDWrjjBbX8C/BnIsdaWtfWxiIi0Fn9dPXsqvOwoDYftXWUedkSCdvPpIwA90uLpl5XE7FE96JlhSEqqBGcp7oYi9rn3sKdmDx/U7Ka2qrZpm6RtLgbEpjKlwUH/+hT6Vx9ggLeavMY523EpkHMa5M0Mf80ZCtmnQUY/3c5dRKSNtFsIN8bEAI8A5wH7gGXGmDettRuarTYbGBx5TAb+Dkw+3rbGmN6R9/a01/GIiLREtbeO3RUe9lR42V3uZW/k654KL4XVvqYb2QBkJ8fRLyuJs07LIS/LSWpyNQ5XOX57gELvXnZV7+Lzmt2U7/7ipjYGQ8+4NPo6EvgKifQPhhhQfYD+Pje59fUYgMTscMAePDl8KcDGwJ3SQyPbIiLtrD1HwicB26y1OwCMMS8ClwLNQ/ilwD+stRZYbIxJN8b0APodZ9v/B/wP8EZ7HIiIyKHqGyxF1T72VHjZU+5ld4W36fmeCi/VvrqD1s9OdtEnM4FJ/TPpnekiI9VDbHw5QVNCsW8vu2t2s7ZmN+/vLcTyRULPcqbQNzaZs6yLPmTQr6aMvjUl9A5FTo4ESOsdHsnufS6bK6DbmZdAzhBdb1tEpANpzxDeC9jb7PU+wqPdx1un17G2NcZ8FdhvrV1tNJIjIm3IGwyxt8LH7vLwiHbzUe19lT6C9Q1N68Y6DHkZCfTJSmJM7zR6pceSmFyNw1lOwJRS7N3PXvdettbuI//AfkKFX5wUmRwTT9+4NMZYJ5fG5NDHXUm/qiL6BHykNA6Zx6dD9mDImxGet934yBwAcYlN+yrKz2dI36nt1SIRETlB7RnCj5SQ7Qmuc8TlxphE4BfA+cf9cGNuBm4G6NatG/n5+cfbpEtyu92n7LG3BfWzdUW7n8F6S5nPUuZriHw9+HlN8OB/shJiITfRQU6CYUgfB+kJfuLiK7Cx5fhNOeX1ZZSFyljgLqOmuuagbROJo5tJJK/BwdS6JPoF3Az2ljPA7yaroQEDNBgnvoQeeBN74c0dy/6EnngTe+FL6EWdM+WLKSQNQClQWgYcfEpMtHva1aifrUv9bF3qZ+try562ZwjfB/Ru9joPKDzBdeKOsnwg0B9oHAXPA1YYYyZZaw8037G19jHgMYAJEybYmTNnfsnD6Zzy8/M5VY+9Laifraut++mvq2dfpY99ld7I14Ofl7kPPgnSGWPolZ5AXk4ikzIS6JkeR1qqD2dcBSFHKeWBIva697Kvdh+ravfirnOD74vtc+PSyItNZpQjnd7BBHp7a+hddYDe/lrSIkEb4whPH8nsD/3PDY9uZw2ErEE40nqT5Igh6Uscs35GW5f62brUz9alfra+tuxpe4bwZcBgY0x/YD9wFXDNIeu8CdwWmfM9Gai21hYZY0qPtK21dj2Q27ixMWYXMEFXRxE5Nfnr6tlf9UW43ltx/JDdMz2BvIwEzhmaS690F2mpfuJc1dTHlOOuL+GAp4j97v2sdRfyfmEhoX1fTBuJNbH0cqWTF5PImJhMetcn0ttdQe+qInoF/STYyLniMa7wlUYy+0OfcyCjf3jaSGb/cACPjWvHLomISEfQbiHcWhsyxtwGvE/4MoNzrbXrjTG3RN5/FHiX8OUJtxG+ROGNx9q2vWoXkeirb7CU1gbYX+WjsMpHUbWPwio/hVU+Cqt9FFX5KfcED9rm8JAdT3qqH6erChtbgTtUSpG3kEJ3Ievd+/mgqJC6/QefQJnpTCXPmcJw4+RcZ0/yAh5615TTu7qQ7qEQTRfwi0uBzH6QNQoGffWLkJ05AFJ6gsPRLn0SEZHOoV2vE26tfZdw0G6+7NFmzy3w/RPd9gjr9PvyVYpIe7PWUuMLsaemng83FlNY5WN/lT8StMNhu7jGT6jh4DnZSXEx9MpIoEdaAqN6pdMj1UVGagBnfDU2pgJPfSmFnv0UugvZ4N7PhweKCOw/eDQ805lKz7hUhpg4znbl0TPgpZenil7VB+jhd5PQ/NqBidnhYN1rKoxsNpqd0T98t0idHC4iIidId8wUkTbnC9ZzoMZPUZWP/VU+iqobR7AjX6t8eBtvs/75ciB8dZEe6fH0SAtfxq9HWjw5qQ4SEmqJiaumzlRQFSyh2FtMkbuIdZ4i5h84gH+//6DPTnem0DMujcEOF2e58ugV9NPLU0XPmgP09FSF7wzZyJUK6X0hYwj0Ox8y+kZe9w1PG3Elt1fLRESki1MIF5GTZq2l0ltHUbWP4ho/B6oDHKjxc6Dax4GaAMXV4dHsGn/osG2zk+PomZ7AoJxkpg/OpkdaHIX71zBiZDcaYirxNpRT7D1AkbuI/Z4iVlQXU1Fccfh+4tLp7kxhsCOO6fG96BUM0stXS8/qYnrWlpDcPGTHxodvuZ7eF3qd+UXATu8bXp6QodFsERFpFwrhInJEwVADJbX+o4frGh/FNQGCoYaDtjMGcpJddE+Lp09WIpP6Z9I9LZ70pHri4qshpoqgqaDMX0yRp4gDngMs8hygeFcxIRvi1RVf7CsxJoEernS6xyQywplNd5NOD7+XHt4qulcX081XTVzzG+XGxEFaXvgxcEwkYPf7Imgn5ypki4hIh6AQLnKKsdZS5a2jpDYQCdmBSLhuDNs+DlQHKPcEDrqVOoAr1kH3tHi6pcYzvk8G3VPjyU2JIyUpgNNVi42pJmArKPPvoMRbQrGnmLXeEuYXFYcv39dMrIkhNz6T7rEpjI1JoXtSMknlpQyJNXRzl9GjuoiU+tDBNwlIyAwH7PTToO/Z4SkiaXnhr+m9w3O2dQKkiIh0AgrhIl1EQ4Ol3BOkpNZPSU2g2dcvwnZpbfjR/M6OjdITnXRPjad7Wjwje6bRLTWe3NQYXPEeYpw1NDiqcNeXUOINP4q9xWz2llBSXkJdw8FXFHEYB9lx6eQ6U+gX42JSQi96OhvoHvDS3VNFj+oisr1VxLCz2Uax+OIySeh2GuQNgRF54WDdGLLT8iDuy1wxW0REpONQCBfp4EL1DZS5gxTX+JsCdWO4Lo2E65JaP2XuIPUNh96EFtISnHRLdZGbEs+A7CRyIs/TkkLEOmuxMVWEHFVU+Pc1Bewd3mIWVRRTUXj4HOz4GBe5rgy6OZMZG5tGbnI63UJBuvk95HqqyK0pJttbefg/LgmZkNoLUvtA3pnNwnUkYKd0Z8knn+pGEyIickpQCBeJAmstNf4QZe7wyHSZO0BZbYAyd7Bp1LoxZJd7godNCwHISoojJ8VFt9R4hnZPITfVRU6yi9TEEE6XG2JrCVFFdbCUEl8JZd4ySn2lbPWVUVJRgjfkPWyf6XFp5LrS6RabyPCEPLq5etEtGCDX7ybXU0m36gOk+qs5bFZ1Yjak9YK0AdB7GqT2jATuZl+dCW3SSxERkc5IIVyklTQG66ZQHQnWpe4AZbXBpmWltQHKPMHDTmgEcBjITnaRm+qiR1o8Y3unkZMST05yHEkJ4XnXxNQQpJrKQBGl3lJKfaUc8Jay1l1KWWkZgfrAYftNiIkn25VOjjOZ02ISmZY8kNz6eroFg+T6a+nmriSn5gDxwT2HbGnCJzOm9oSMwdD3rEiobgzYPSGlBzjj26irIiIiXZNCuMgxNN5EptQdODhER4L1pt1+Hli3sGkU+0hzrWMchsykOLKTXeSkuBiYm0xOsovsZBeZyTHEu/wYZw31pgZ/QyVl/j1No9a7fWUsry6l7EAZoYbDL/OX4kwOh+vYZMbEZZHjzCGnvp7soJ9cfy3Z7kpy3CUk+fYcPnod44KU7uEQnTsSBp578Oh1Wi9I7q5bqouIiLQBhXA55fiC9ZR7AlR4gpS7g5R7glR4As2eBw+aHnK0YJ2VFEc8ln4ZcQzOTSE7JY6cZBcZSbHEu/zEOGuxjlqCtpqKQDHlvnLK/GWU+8rZ5S+nvKicqkAVDfYIJ0m60shxZZATm0T/+B7kxPciJ1RPdp2fHF8tOZ5KsmuKSfAdOnINOJzhYJ3SHbKHQv9ZX4Tt5l/j03W5PhERkShRCJdOz19XT7knSLk7EPkaCdVNz4NN71d4gl/cmfEQcbEOspPiyEyOIyvJxeDcFHJSXGQnx5GZHEu8y4eJrYUYN4GGKioCB1i1ZRWJOYmU+crY4SunorCCSn8llsMnccfHxJPlSicrLoW8mHjGJPUjK5FwuA76yfG7yfFWkV1bitN7hHBtYiIhujtkDII+08LPkw8J2AkZukyfiIhIB6cQLh1OY6iucAcp8wSocAcpj4TqishodfPR66OG6hgHWclxZCbFkZXsYkB2UuR5HFlJcaQmGGLjPDhiPDQ4avGEqij3F1LuL6fcF37s8JdTvj88Yn2kYB1n4shtyCErLpk+scmMS8kgKwmyIlNCsgJesrxVZLnLSXTvxBxyKT8gMi2kGyR3g/QBkDfl8FHrlB6QmKVwLSIi0kUohEub8tfVU+kNUumpo9IbHpWu8gapiLz+YlkdFZ7w62OF6sYQnZkUR/+sRLKSXZH51pFQ7fTicLppMG589dVUBkqo9FdS4Q+PUG/yV1BRXkFlYSWeOs8RPychNoHs+Cyy4lLo60xjfFw2WdaEg3VdgKyAhyxvNVnuMlzV+4hp2Hb4TkwMJOWET2pM7ga5oyPPc79Yltwt/NyVqmkhIiIipxiFcDlhjYH60NDcPGA3BuvGZUcL1AAp8bFkJMaREQnRg3OTyUiKjFwnxZGaAM44H45YDw0xbvz11VT6S8OBOhAO1pv8FVRWVlJRVHHUUB3riCXTlUGmK42MmER6J/YgIzGPzAbIrA+RURcgM+Aly1dDlqeCRE8RBDYfueiEzC+CdN4A9iYF6T3s9C8CdXI3SMqFxExwxLRG20VERKQLUgg/BVlrcQdCVHnrqPbVUeU9+qh0Y6Cu8ATx1R07UGcmxZGRGD458bRuKWRGAnZGYhzpiTG44uqIcXrB4SFkanHX1VDlL6IyUElVoIoKfwWb/ZVNofrQ25w3ijWxZMRnkOlKJyM2iV6JPclK6kNGgyGzvp6MUJCsgJcMfy0ZnipSPGUY7y44wgmQGEf4GtdJOZCUDT37f/E8KeeQcJ1z2JVCtufn0/vMmV/iuyEiIiKnIoXwTqyhwVLrD1HlCwfnKl8dVd5gU7AOLwtS3ey90moP3nnvETrCnRUbpUYCdXpiHLkp8QzplkpGorMpUGckxhLvCuKI9WFivISopbaumqrAASr94UBdFahiq7+SquoqqkqqqA5UU2+PHOJjHbFkuDJIj0sj05nEiMSeZCb3JcM6yGywZIaCZAR9ZPrdZHirSPWUYzwboe7wm80AEJcSDtHJuZA5AHpPigTrZuG68ZGQoRFrERERaXcK4R1AqL6BGn+IKm+QKl9dJDQHm4J0dSRAV/kOfl3tq+MYWZpkVyxpCU7SE8OPod1T6R0fYPjAvpFlcaTGxxDnDOKM80GMlxBu3HU1VPqLqQpUhUep/VVsC1RSXVNNZUkl1YFqQvbwa1ZDeJQ6PT6ddFc66bFJDIzPIT2xF+k4yLCQHgqRXhcO1en+WjJ81SR6KjDeDRDyHflAHLEHB+is0w4O0geF62zdmVFEREQ6PIXwdrJqbxWvFuw7wmh1kBr/kQNto9T4WNIT40hPdJKW4KR3ZiLpTeE6jvQEJ6kJMbjigpiYyOi08eCpq6E6UE51sJrqQPhR5tzB6oZYqsurqS6spiZYc8TrVAPEmBjSXGmRUeoU+iXkMDYxj3QTQ7qFjPp60kNBMgI+0gMe0n01JHvLMUU7wVcFR7iaCBAZqc4KX+0juXv4RjGJmeEAnRhZntj4PDM8Wq0TF0VERKQLUQhvJ/sqvby9ppD0xDjSEpxkJsUxIDup6XVGJFCnJMTgcoaneliHh3o81DaG6UB1U6AuD1SzI1BNdWn49bHCNECyM5k0VxqOegd5ju70TO5NWspA0kwMaRbS6htIrw+RXhcgozFQeypxlJaAdyMc4VboQPjGMI3BOSkLuo8+PEQ3hevs8OtYVxt1WURERKRzUAhvJ+MGWO6+ro6qQCk1gZqmMF0cqGJroIZqdzVVgSpqg7XHDNMpzhRSXamkudJIi02mV0oGqWlO0giH6fQGS1p9iLRggLSgl7SAhxRfNc6aSvBtwXorMRx9/8SnfRGY03pDj7FfjFonZh8cuBOzdHk9ERERkZOgEN5OVhSv4P9b+P81vW4M0+mudNLiUugVn0mqI450YknDkFbfQFp9HWmhIGkBH2kBD6n+GmIrK8G3E3wVUB88+gc6k76YypGQAWl5kJDBnpJa+g4bF1me2WydTEhI10mKIiIiIu1AIbydzHCk8mb6VNL8blL9tcTWVIJ3L/jWwFGubw1ATNzBYTlr4MHB+aAQnfHF66NM+diZn0/fqTPb5iBFRERE5IQohLeTNE85aZvnfxGcU3tBt1GR4HxoiG723Jmo6R4iIiIiXYxCeHsZeUX4ISIiIiKnPEe0CxAREREROdUohIuIiIiItDOFcBERERGRdqYQLiIiIiLSzhTCRURERETamUK4iIiIiEg7UwgXEREREWlnCuEiIiIiIu1MIVxEREREpJ0phIuIiIiItDOFcBERERGRdqYQLiIiIiLSzhTCRURERETambHWRruGdmeMKQV2R7uOKMkGyqJdRBeifrYu9bP1qaetS/1sXepn61I/W19r9LSvtTbn0IWnZAg/lRljlltrJ0S7jq5C/Wxd6mfrU09bl/rZutTP1qV+tr627Kmmo4iIiIiItDOFcBERERGRdqYQfup5LNoFdDHqZ+tSP1ufetq61M/WpX62LvWz9bVZTzUnXERERESknWkkXERERESknSmEd3LGmK8bY9YbYxqMMYedvWuM6WOMcRtjftJs2enGmLXGmG3GmIeMMSay3GWMeSmyfIkxpl+zba43xmyNPK5vl4OLgqP10xhznjGmINK3AmPM2c3eUz+P4lg/n8aYn0d6s9kYc0Gz5ernCTLGjDXGLDbGrDLGLDfGTGr2Xqv191RijLk90rP1xpg/NVuufp4kY8xPjDHWGJPdbJn6eRKMMX82xmwyxqwxxrxmjElv9p562oqMMRdGernNGHNnm3yItVaPTvwAhgFDgHxgwhHefxX4F/CTZsuWAlMBA7wHzI4svxV4NPL8KuClyPNMYEfka0bkeUa0j709+wmMA3pGno8E9qufX6qfw4HVgAvoD2wHYtTPFvd3XrP+XATkt3Z/T6UHMAuYD7gir3PVzy/d097A+4TvzZGtfn7pfp4PxEae/xH4o3raJn2OifRwABAX6e3w1v4cjYR3ctbajdbazUd6zxhzGeFAsr7Zsh5AqrV2kQ3/pP0DuCzy9qXAM5HnrwDnRH5jvgD4wFpbYa2tBD4ALmyDw4m6o/XTWrvSWlsYebkeiI+MIqifx3CMn89LgRettQFr7U5gGzBJ/WwxC6RGnqcBjT+jrdnfU8n3gHuttQEAa21JZLn6efL+H/A/hH9WG6mfJ8laO89aG4q8XAzkRZ6rp61rErDNWrvDWhsEXiTcr1alEN5FGWOSgJ8Bdx3yVi9gX7PX+yLLGt/bCxD5j7wayGq+/AjbnIq+BqyM/I9a/Tw5R+uB+tkyPwT+bIzZC9wH/DyyvDX7eyo5DZge+dP8x8aYiZHl6udJMMZ8lfBfDVcf8pb62TpuIjyyDeppa2uX/6/EtvYOpfUZY+YD3Y/w1i+stW8cZbO7gP9nrXUf8ovtkX7Ltcd571jbdDon2c/GbUcQ/hPg+Y2LjrCa+hl2rH6eTG9OiX4e6lj9Bc4BfmStfdUY8w3gSeBcWre/Xcpx+hlLeErTFGAi8LIxZgDq51Edp5//H1/8W3nQZkdYpn5GnMi/qcaYXwAh4LnGzY6wvnp68tqlNwrhnYC19tyT2GwycGXkxKJ0oMEY4yc8Rzyv2Xp5fPEn7H2E5+/tM8bEEv7zdkVk+cxDtsk/iZo6hJPsJ8aYPOA14JvW2u2RxftQP0+mn429adTYt1O+n4c6Vn+NMf8A7oi8/BfwROR5a/a3SzlOP78H/DvyZ/ulxpgGIBv186iO1k9jzCjCc5NXRwaC8oAVJnzysPp5DMf7N9WETz6/GDgn8rMK6mlrO1o/W5Wmo3RR1trp1tp+1tp+wF+A31tr/2qtLQJqjTFTInO/vgk0jla+CTReWeJK4KPIf+DvA+cbYzKMMRmERzbeb8fDibrIGejvAD+31n7WuFz9PGlvAldF5tX3BwYDS9XPFisEzoo8PxvYGnnemv09lbxOuI8YY04jfEJWGepni1lr11prc5v9f2gfMN5aewD186QZYy4kPNX0q9Zab7O31NPWtQwYbIzpb4yJI3zi6put/ilf5qxOPaL/AC4n/I9bACgG3j/COr/l4KujTADWET7z9698cdOmeMKjadsIn009oNk2N0WWbwNujPZxt3c/gV8CHmBVs0fjlRPUz5P4+ST85+rtwGYiZ+urny3u7zSggPCZ+0uA09uiv6fKg3Do/mekPyuAs9XPVuvtLiJXR1E/v1QftxGeq9z4/6FH1dM26/VFwJZI337RFp+hO2aKiIiIiLQzTUcREREREWlnCuEiIiIiIu1MIVxEREREpJ0phIuIiIiItDOFcBERERGRdqYQLiIiIiLSzhTCRURERETamUK4iIhgjPmuMcYaY85qtuy2yLJzI6/vN8ZsMMY8HL1KRUS6hthoFyAiIh3CaGANMAz42BiTCHwLKAXWGmMGAGdaa4dHsUYRkS5DI+EiIgIwCngBGBp5/QPCt7RuANKBj4G+xpiVxpikqFQoItKFKISLiAiER8BfBoYaY9KAOcDnwDpr7WbgGeBX1tpx1lpPFOsUEekSNB1FROQUZ4zpDZRba3cYY3KB/wEeBk4jPEUFwiPlb0TWvwz4CpALPAIkNn9trZ3XrgcgItIJKYSLiMhoYG3keS1wITAJ+AuwIrJ8BLAewFr7OvC6MSYDuM9a+63mrwGFcBGR41AIFxGRUXwRwv9MeFS83hgzCnjaGJMC1FlrvYds90vCI+FHey0iIkehEC4iIqOAVwGstW83Wz4c2ACMBdY1LjTGGOBe4D1r7YpDX7dX0SIinZmx1ka7BhER6USMMT8ArgeWAauAuOavrbWPRq86EZHOQSFcRERERKSd6RKFIiIiIiLtTCFcRERERKSdKYSLiIiIiLQzhXARERERkXamEC4iIiIi0s4UwkVERERE2plCuIiIiIhIO1MIFxERERFpZwrhIiIiIiLt7P8HnVdc2hb5Xv0AAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots(figsize=(12, 6))\n", + "lstyles = [\"-\", \"--\", \":\"]\n", + "lwidths = [1, 2, 4]\n", + "for idx, averaging_method in enumerate(omega_averaging_methods):\n", + " fref_out, ecc, mean_ano, eccMethod = measure_eccentricity(\n", + " fref_in=fref_in,\n", + " dataDict=dataDict,\n", + " method=\"ResidualAmplitude\", \n", + " return_ecc_method=True,\n", + " extra_kwargs={\"debug\": False, \"omega22_averaging_method\": averaging_method})\n", + " ax.plot(eccMethod.t_for_omega22_average, eccMethod.omega22_average / (2 * np.pi), label=f\"{averaging_method}\")\n", + "ax.legend()\n", + "ax.set_xlabel(r\"$Mf_{22}$\")\n", + "ax.set_ylabel(\"eccentricity($f_{22}$)\")\n", + "ax.grid()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5c8dff5b-12b0-4a41-87d2-c2477cbc68af", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/test/test_interface.py b/test/test_interface.py index cd071cd4..4177427f 100644 --- a/test/test_interface.py +++ b/test/test_interface.py @@ -1,6 +1,7 @@ import measureEccentricity from measureEccentricity import load_data from measureEccentricity import measure_eccentricity +import numpy as np def test_interface(): @@ -37,3 +38,16 @@ def test_interface(): # Make diagnostic plots eccMethod.make_diagnostic_plots(usetex=False) + + # Try evaluating at single frequency + tref_out, ecc_ref, meanano_ref = measure_eccentricity( + fref_in=0.025 / (2 * np.pi), + method=method, + dataDict=dataDict) + + # Try evaluating at an array of frequencies + tref_out, ecc_ref, meanano_ref, eccMethod = measure_eccentricity( + fref_in=np.arange(0.025, 0.035, 0.001) / (2 * np.pi), + method=method, + dataDict=dataDict, + return_ecc_method=True)