diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..4192455 --- /dev/null +++ b/.gitignore @@ -0,0 +1,14 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# Environments +/env \ No newline at end of file diff --git a/bisection.py b/bisection.py new file mode 100644 index 0000000..26c2ae7 --- /dev/null +++ b/bisection.py @@ -0,0 +1,76 @@ +import numpy as np +from functions import * + + +class bisector: + xinf, xsup = 0, 0 + xmed, xmed_old = 0, 0 + err_rel, err_obj = float('inf'), 0 + iter_n, max_iter = 0, 500 + + def __init__(self, function, xinf, xsup, err_obj=1e-4, max_iter=500, bias=0, verbose=False, *args, **kwargs): + self.xinf , self.xsup = xinf, xsup + self.err_obj , self.max_iter = err_obj, max_iter + self.function, self.bias = function, bias + self.args , self.kwargs = args, kwargs + self.verbose = verbose + + def next_iter(self): + self.iter_n += 1 + self.xmed_old = self.xmed + self.xmed = (self.xinf + self.xsup) / 2 + self.calc_error() + + self.yinf = self.function(self.xinf, *self.args, **self.kwargs) - self.bias + self.ymed = self.function(self.xmed, *self.args, **self.kwargs) - self.bias + self.logs('iter') + + if self.ymed == 0: + self.err_rel = 0 + return + elif self.yinf * self.ymed < 0: + self.xsup = self.xmed + elif self.yinf * self.ymed > 0: + self.xinf = self.xmed + else: + raise Exception('Fail iter_xmed') + + def calc_error(self): + self.err_rel = 100 * (self.xmed - self.xmed_old) / self.xmed + self.err_rel = abs(self.err_rel) + + def run(self): + self.logs('header') + while((self.err_rel > self.err_obj) and (self.iter_n < self.max_iter)): + self.next_iter() + + if self.err_rel == 0: self.logs('exact') + else: self.logs('results') + return(self.xmed) + + def logs(self, item): + if self.verbose is True: + if item == 'header': + print("Iter | xinf{0:4} | xsup{0:4} | xmed{0:4} | eps{0:5} | y(xmed)".format('')) + print("{:=>60}".format('')) + if item == 'iter': + print("{:>4} | {:^8.4f} | {:^8.4f} | {:^8.4f} | {:_>8.4f} | {:^8.4f}". + format(self.iter_n, self.xinf, self.xsup, self.xmed, self.err_rel, self.ymed)) + if item == 'results': + if not (self.err_rel > self.err_obj): print("Finish by eps: {:.4e}".format(self.err_rel)) + if not (self.iter_n < self.max_iter): print("Finish by iter: {}".format(self.iter_n)) + print("Xroot: {0}".format(self.xmed)) + if item == 'exact': + print ("Exact root found: {}".format(self.xmed)) + + +def function(tInt_perc): + frec, ampl = 1, 1 + rms_unit = (np.sqrt(2) / 2) * ampl + rms_dimmer = fun_rms_simbolic(tInt_perc, frec=frec, amp=ampl) / rms_unit + return rms_dimmer + + +solver = bisector(function, xinf=0, xsup=100, bias=0.1, err_obj=1e-5, verbose=False) +sol = solver.run() +print("rms({:.4f}) = {:.4f}".format(sol, function(sol))) \ No newline at end of file diff --git a/functions.py b/functions.py new file mode 100644 index 0000000..f6c9bad --- /dev/null +++ b/functions.py @@ -0,0 +1,51 @@ +import numpy as np + + +def fun_Sine(t, f=1, a=1): + w = 2 * np.pi * f + return a * np.sin(w * t) + + +def fun_Dimmer(t, tInt_perc, f=1, a=1): + periode = 1 / f + tInt = 0.5 * periode * (tInt_perc / 100) + vt = fun_Sine(t, f, a) + vd = 0 + if (0 <= t and t < tInt): + vd = 0 + elif (tInt <= t and t < 0.5 * periode): + vd = vt + elif (0.5*periode <= t and t < 0.5*periode + tInt): + vd = 0 + elif (0.5*periode + tInt <= t and t < periode): + vd = vt + return vd + + +def fun_rms_numeric(data_y): + sum_y = 0 + for i in data_y: + sum_y = sum_y + i**2 + rms = np.sqrt(sum_y / len(data_y)) + return rms + + +def fun_rms_simbolic(tInt_perc, frec=1, amp=1): + periode = 1 / frec + tInt = 0.5 * periode * (tInt_perc / 100) + T1, T2 = 0, periode + + def sine_integrate(t): + w = 2 * np.pi * frec + int_a = amp ** 2 + int_b = t / 2 + int_c = (np.sin(2*t*w) / (4*w)) + sine_int = int_a * (int_b - int_c) + return sine_int + + rms_a = 1 / (T2 - T1) + rms_b1 = sine_integrate(T2 / 2) + rms_b2 = sine_integrate(tInt) + rms_b = (rms_b1 - rms_b2) * 2 + rms = np.sqrt(rms_a * rms_b) + return rms \ No newline at end of file diff --git a/main.py b/main.py new file mode 100644 index 0000000..4d382a7 --- /dev/null +++ b/main.py @@ -0,0 +1,117 @@ +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.widgets import Slider, Button +from scipy.optimize import bisect + +from functions import * + + +def main(): + for i in range(100+1): + print("duty([per/2] * {:>14.10f}%)= {:>6.2f}%".format(solve_tInt(i), i)) + + init_plot() + + +def gen_data(tInt_perc=50, frec=1, amp=1): + samples = 100 + periode = 1 / frec + t = np.linspace(0, periode, samples) + v_dimmer = [fun_Dimmer(ti, tInt_perc, frec, amp) for ti in t] + rms_dimmer = fun_rms_simbolic(tInt_perc, frec, amp) * np.ones(samples) + return [t, v_dimmer, rms_dimmer] + + +def init_plot(): + global fig, ax1, l1, l2 + t_series, vt_series, rms_series = gen_data() + plt.ion() + + fig, ax1 = plt.subplots(figsize=(6, 6)) + plt.subplots_adjust(bottom=0.35, top=0.95, right=0.89, wspace=0, hspace=0) + + l1, = plt.plot(t_series, vt_series) + l2, = plt.plot(t_series, rms_series) + + plt.grid(True) + update_plot() + setting_sliders() + ax1.set_ylabel(r'Amplitud [v]') + secax_y = ax1.secondary_yaxis( + 'right', functions=(rms_to_duty, duty_to_rms)) + secax_y.set_ylabel(r'Duty [%]') + plt.show(block=True) + + +def update_plot(tInt_perc=50, frec=1, ampl=1): + global fig, ax1, l1, l2 + t_series, vt_series, rms_series = gen_data(tInt_perc, frec, ampl) + # Plot data + l1.set_xdata(t_series), l2.set_xdata(t_series) + l1.set_ydata(vt_series), l2.set_ydata(rms_series) + # Update legend + sine_str = r'$v(t) = \mathcal{A} \mathrm{sin}(2 \omega t)$' + rms_str = "RMS: {:.5f}".format(rms_series[0]) + ax1.legend([sine_str, rms_str]) + # Update limits + ax1.set_xlim(np.min(0), np.max(1 / frec)) + ax1.set_ylim(-ampl, ampl) + + +def setting_sliders(): + def update_sliders(val): + update_plot(slide_tInt.val, slide_frec.val, slide_amp.val) + frec, ampl, tInt_perc = slide_frec.val, slide_amp.val, slide_tInt.val + rms_unit = (np.sqrt(2) / 2) * ampl + rms_value = fun_rms_simbolic(tInt_perc, frec=frec, amp=ampl) / rms_unit + slide_duty.set_val(rms_value * 100) + + def update_slider_duty(val): + tint = solve_tInt(slide_duty.val, slide_frec.val, slide_amp.val) + slide_tInt.set_val(tint) + + global ax_frec, ax_amp, ax_tInt, ax_duty, ax_solve + global slide_frec, slide_amp, slide_tInt, slide_duty, btn_solve + ax_frec = plt.axes([0.1, 0.10, 0.80, 0.03]) + ax_amp = plt.axes([0.1, 0.15, 0.80, 0.03]) + ax_tInt = plt.axes([0.1, 0.20, 0.80, 0.03]) + ax_duty = plt.axes([0.1, 0.25, 0.6, 0.03]) + ax_solve = plt.axes([0.8, 0.25, 0.1, 0.03]) + slide_frec = Slider(ax_frec, 'Frec', 1.0, 50, valinit=1, valstep=1) + slide_amp = Slider(ax_amp, 'Amp', 1.0, 311, valinit=1, valstep=1) + slide_tInt = Slider(ax_tInt, '%Int', 0.0, 100, valinit=50, valstep=1) + slide_duty = Slider(ax_duty, '%Duty', 0.0, 100, valinit=50, valstep=1) + btn_solve = Button(ax_solve, 'Solve') + slide_frec.on_changed(update_sliders) + slide_amp.on_changed(update_sliders) + slide_tInt.on_changed(update_sliders) + # slide_duty.on_changed(solve_tInt) + btn_solve.on_clicked(update_slider_duty) + + +def rms_to_duty(x): + vrms = x + # vrms = fun_rms_simbolic(slide_tInt.val, slide_frec.val, slide_amp.val) + vrms_max = np.sqrt(2) * 0.5 * slide_amp.val + vrms_duty = 100 * vrms / vrms_max + return vrms_duty + + +def duty_to_rms(x): + vrms_duty = x + vrms_max = np.sqrt(2) * 0.5 * slide_amp.val + vrms = vrms_duty * vrms_max / 100 + return vrms + + +def solve_tInt(duty, frec=1, amp=1): + def function(tInt_perc, bias=0): + rms_unit = (np.sqrt(2) / 2) * amp + rms_value = fun_rms_simbolic(tInt_perc, frec=frec, amp=amp) / rms_unit + return rms_value - bias + root = bisect(function, 0, 100, args=(duty / 100)) + return root + + +if __name__ == "__main__": + main() diff --git a/rms_vs_tinterrupt.ipynb b/rms_vs_tinterrupt.ipynb new file mode 100644 index 0000000..4ac8391 --- /dev/null +++ b/rms_vs_tinterrupt.ipynb @@ -0,0 +1,206 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 8, + "source": [ + "import numpy as np\r\n", + "import matplotlib.pyplot as plt\r\n", + "from functions import *\r\n", + "\r\n", + "frec, ampl = 1, 1\r\n", + "periode = 1 / frec\r\n", + "rms_unit = (np.sqrt(2) / 2) * ampl\r\n", + "\r\n", + "def function(tInt_perc, bias=0):\r\n", + " frec, ampl = 1, 1\r\n", + " rms_unit = (np.sqrt(2) / 2) * ampl\r\n", + " rms_dimmer = fun_rms_simbolic(tInt_perc, frec=frec, amp=ampl) / rms_unit\r\n", + " return rms_dimmer - bias\r\n", + "\r\n", + "def plot():\r\n", + " t_perc_serie = np.arange(0, 100, 1)\r\n", + " vt_serie = [(function(t) * 100) for t in range(100)]\r\n", + " plt.plot(t_perc_serie, vt_serie)\r\n", + " plt.xlabel(\"Porcentaje semi periodo [%]\")\r\n", + " plt.ylabel(\"Porcentaje RMS normalizado [%]\")\r\n", + " plt.grid()\r\n", + "\r\n", + "plot()" + ], + "outputs": [ + { + "output_type": "display_data", + "data": { + "text/plain": [ + "
" + ], + "image/svg+xml": "\r\n\r\n\r\n \r\n \r\n \r\n \r\n 2021-08-25T02:01:55.290193\r\n image/svg+xml\r\n \r\n \r\n Matplotlib v3.4.3, https://matplotlib.org/\r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n \r\n\r\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAykElEQVR4nO3dd3gVVfrA8e+bXkmooYWOAiIgoQmoYFuxoKtgpaO4FqzrimXV3XXdta5dQaSICCqioLI2BFGkdxBRitKb1IQSQt7fHzPhl43k5qbcO8m97+d55smdc2fuvIcJ983MmXOOqCrGGGMMQITXARhjjCk/LCkYY4w5wZKCMcaYEywpGGOMOcGSgjHGmBOivA6gNKpVq6YNGjQo0b5ZWVkkJiaWbUAVQDjWOxzrDOFZ73CsMxS/3osWLdqtqtVP9l6FTgoNGjRg4cKFJdp35syZdOvWrWwDqgDCsd7hWGcIz3qHY52h+PUWkV8Le89uHxljjDnBkoIxxpgTLCkYY4w5wZKCMcaYEywpGGOMOSFgSUFERonIThFZma+sioh8KSI/uz8ru+UiIi+KyFoRWS4ibQMVlzHGmMIF8kphDHBRgbJhwHRVbQpMd9cBegBN3WUI8FoA4zLGGFOIgPVTUNVZItKgQPHlQDf39VhgJnC/W/6WOuN4zxWRVBGpparbAhHbgl/2MPnnbBZnr/n/QpH/f+muCkKEQESEIAKRIkRGOEtUhBAVGUFUhBATFUFMZAQxURHERkUSFx1BXHQk8TGRJMZEkRDr/IyMkN8HY4wx5UiwO6+l5fui3w6kua/rAJvybbfZLftdUhCRIThXE6SlpTFz5sxiBzFtQzYfr8uGdWsBCNaMErGREB8lxEdBYrScWJJjIDlGqBQjpMQKKTFCaqxQKVaIkLJNJJmZmSX6N6vIwrHOEJ71Dsc6Q9nW27MezaqqIlLs72NVHQGMAGjXrp2WpPdit25wsY8egKqKqpMsVJVchVxVjucqx1U5flzJyVVycnPJOa5kH88lO8dZjubkcvjYcY64S9bR4xzKziHr6HEyjx7j4JEcDhw5xv7DzvJr1jH27Mrm8LGc38URFSGkVYqjVkocdSvHk14lgfTKCdSvmkCDaonUSI5Fipk0wrHHZzjWGcKz3uFYZyjbegc7KezIuy0kIrWAnW75FiA933Z13TJPiEi+u0nBueVzOPs4uzOPsivzKLsOHmXngSNs2+8sW/YdZsEve5m6bCu5+dJoQkwkjaon0qR6Ek3Tkjk1LZlTayZTt3J8sZOFMcZA8JPCVKA/8G/355R85beLyESgI7A/UO0J5VV8TKRzJVAlodBtjh3PZdu+I/y6J4tfdmexfncW63ZlMX/DHj5auvXEdsmxUTSvXYnT66Rwep0UWqen0qBqgiUKY0yRApYURGQCTqNyNRHZDDyKkwzeE5HBwK/A1e7m04CLgbXAIWBgoOKqyKIjI6hXNYF6VRM4q+n/DnCYeTSHn3Yc5MdtB1m97QArt+7n7bm/cjQnF4DUhGjapKdSJTeb+Hq/0To9lbjoSC+qYYwpxwL59NF1hbx13km2VeC2QMUSDpJio2hbrzJt61U+UZZzPJefd2aybNM+lmzcx+KNe5m58xiTf55LTGQEbeql0qVxNTo3qUqb9FSiI60vozHhrkIPnW18i4qMoHmtSjSvVYlrO9QD4JMvZhBbtwXzN/zGnPW/8fz0n/jPV05S6dy4KuecWp1up9agTmq8x9EbY7xgSSHMJMUI3VqkcUEL52ngfYeymbv+N2b9vJtv1uziix92ANCiViXOb5HGhS3SOK12JWuPMCZMWFIIc6kJMVzUshYXtayFqrJuVybTV+9k+uqdvPz1z7w4/WfqVo7n4tNr0aNlTdqkp1qCMCaEWVIwJ4gITWok06RGMjef05g9Wdl89cMOpq3cxujZGxgxaz31qiTQs3VtLm9Tm6ZpyV6HbIwpY5YUTKGqJMZwdft0rm6fzv7Dx/hi1XamLtvKqzPX8vKMtbSqm8JVbevSs3VtKifGeB2uMaYMWFIwfkmJj6Z3u3R6t0tn18GjfLxsK5MWbebRqav456er+UPLmlzXIZ0zG1W120vGVGCWFEyxVU+OZVDXhgzq2pBVW/fz/sLNTF68mY+XbaVhtURu6FiP3hnppCREex2qMaaYLCmYUjmtdgqn9UxhWI9mTFuxjfHzNvL4p6t59oufuOKMOgzq0sDaHoypQCwpmDIRFx3JlW3rcmXbuqzcsp9xc35l8uLNTJi/kXNOqc6NZzWka5NqdmvJmHLOurCaMteyTgpP9mrFnAfO494LTmHV1gP0fXM+l770HZ8s38rx3GANVm6MKS5LCiZgqiTGMPS8pswe1p2nrmrF4WPHuf2dJZz77EzeXbCRY8dzvQ7RGFOAJQUTcLFRkVzdPp0v7z6H1/u0pVJcNPd/sIJuT8/knXkbyc6x5GBMeWFJwQRNZIRwUctaTL29C6MHtKdaciwPfriCc5+dyfsLN5FjVw7GeM6Sggk6EaF7sxp8dGtnxgxsT+WEGO6btJwLn5/FtBXbcAbNNcZ4odCnj0TkRT/2P6CqD5dhPCaMiAjdTq3BOadU5/NVO3j2izXcOn4xbdJTeaBHMzo2qup1iMaEHV+PpF4OPFLE/sMASwqmVESEi1rW5IIWaXywaDPPffkT14yYy/nN03jokuY0rJbodYjGhA1fSeE/qjrW184iUtnX+8YUR2SEcHX7dC5rXZtRszfw6oy1XPDcN/Tv3IA7zm1qPaSNCYJC2xRU9fmidvZnG2OKKz4mktu6N2HGfd3olVGXUbM30P3ZmUycv5Fc6+NgTED53dAsIpeJyEwRmSsitwYyKGMAaiTH8e+rWvHJ0K40rp7IsMkr+OOrs1m2aZ/XoRkTsgpNCiLSpkBRX6A70Bm4JYAxGfM/Tqudwns3n8nz17Rh2/4jXPHqbP760Ur2Hz7mdWjGhBxfbQq3iEgE8FdV3Q5swmlUzgW2BiM4Y/KICFecUYfzmtfguS9/Yuz3v/Dfldv566XN6dm6to2pZEwZ8dWmcDPwMjBcRB7BeRJpDrAC6Bmc8Iz5X8lx0Tx62WlMvb0rdVLjuHPiUgaNWcCWfYe9Ds2YkOCzTUFVl6nq5cASYApQW1WnqurRoERnTCFa1klh8q1deOTSFszbsIcLnvuGMbM3WEO0MaXkq03hTyLyvYh8DyQCFwGpIvK5iJwdtAiNKURkhDCoa0O+uPts2jeowmMf/8C1I+byy+4sr0MzpsLydaVwq6p2xmlcvk9Vc1T1ReBa4IpgBGeMP+pWTmDMwPY807s1q7cf4KIXZjHqO7tqMKYkfDU0bxGRB4EE4Me8QlXdC9wT6MCMKQ4RoVdGXbo2qcYDk5fz909+4MsfdvDM1a29Ds2YCsXXlcLlOI3K3wH9ghOOMaVTMyWOUQPa8+RVp7N88z4u+s8svttyzAbZM8ZPvpJCFVX9WFU/U9XjJ9tARGoGKC5jSkxEuKZ9PT6762ya16rEyBXZ3P7OEvYfsn4NxhTFV1KY5sf+/mxjjCfSqyQwYUgnep8SzeertnPRC7OYs+43r8MyplzzlRRai8gBH8tBIC1YgRpTEpERwiWNYvjw1i7ER0dy/ci5PPvFGpvQx5hC+Oq8FqmqlXwsyapaJ5jBGlNSp9dN4ZM7utI7oy4vfb2Wa0fMtQ5vxpyEJzOvicjdIrJKRFaKyAQRiRORhiIyT0TWisi7IhLjRWwmdCXERPFUr9a8cG0bVm87wMUvfMtXP+zwOixjypWgJwURqQPcAbRT1ZZAJE7fhydx5nBoAuwFBgc7NhMeLm9Th0/vOIu6leO58a2F/Gvaao7Z7SRjAO/maI4C4kUkCqcfxDbgXGCS+/5YrIOcCaAG1RL54JbO3NCxHsNnref6N+ay48ARr8MyxnPiz/PbItIaOMtd/VZVl5XqoCJ3Av8EDgNfAHcCc92rBEQkHfiveyVRcN8hwBCAtLS0jIkTJ5YohszMTJKSkkpWgQosHOtdVJ3nbs1h9KqjxEYKt7aJpVmVyCBGFzh2rsNHcevdvXv3Rara7qRvqqrPBecLeyXwd3dZAQwtaj8fn1cZ+BqoDkQDHwF9gLX5tkkHVhb1WRkZGVpSM2bMKPG+FVk41tufOv+0/YB2f2aGNnrgUx3+zVrNzc0NfGABZuc6fBS33sBCLeR71Z/bR4OBjqr6iKo+AnQCbvI7Jf3e+cAGVd2lqseAyUAXnMH28obdqAtsKcUxjCmWpmnJTLmtCxc0T+OJaT9y+4QlHMrO8TosY4LOn6QgQP4ezcfdspLaCHQSkQRxZkY5D/gBmAH0crfpjzNUtzFBkxwXzWt92nL/Rc2YtmIbV776PRt/O+R1WMYElT9JYTQwT0QeE5HHgLnAmyU9oKrOw2lQXoxzKyoCGAHcD9wjImuBqqU5hjElJSLc0q0xYwZ2YOu+w/R85Ttmr93tdVjGBE2RSUFVnwMGAnvcZaCqPl+ag6rqo6raTFVbqmpfVT2qqutVtYOqNlHV3moT+RgPnXNKdT4e2pUaybH0GzWfMbM32KB6JiwUOnS2iFTJt/qLu5x4T1X3BC4sY7xXv2oik2/twl0Tl/LYxz+wZsdB/tazJTFRXj3JbUzg+ZpPYRGgOO0H9XA6lAmQitMu0DDQwRnjtaTYKEb0zeDZL9fwyox1rN+Vxet9MqicaB3uTWjyNfZRQ1VtBHwFXKaq1VS1KnApTt8CY8JCRIRw3x+a8fw1bViycR9/fHU263Zleh2WMQHhz3VwJ1U9MUS2qv4X6By4kIwpn644ow4ThnTk4JEc/vjKbL63BmgTgvxJCltF5GERaeAuDwFbAx2YMeVRRv0qfHRbF9IqxdFv1HwmLdrsdUjGlCl/ksJ1OL2PP3SXGm6ZMWEpvUoCk27pTKdGVfnz+8t47os19mSSCRm+GpoBcJ8yujMIsRhTYaTERzN6YHse+nAFL369lk17D/PkVa3sySRT4RWZFESkOvAX4DQgLq9cVc8NYFzGlHvRkRE8eVUr6lVJ4JkvfmLnwSO81ieDSnHRXodmTIn582fNeOBHnEdQ/4bTX2FBAGMypsIQEW4/tynP9m7NvPV7uPr1OWzfb0Nwm4rLn6RQVVXfBI6p6jeqOghn7gNjjOuqjLqMHtiezXsPc+Wrs1m786DXIRlTIv4khWPuz20icomInAFU8bWDMeHorKbVmTikE9nHlV6vz2HRr3u9DsmYYvMnKTwuIinAvcCfgZHA3QGNypgKqmWdFCbf0pnU+GhuGDmX6attDmhTsfiTFJap6n5VXamq3VU1A5gf6MCMqajqVXUeWT0lLZkh4xbxgfVlMBWIP0lhg4hMEJGEfGXTCt3aGEO1pFjeuakTnRpV4d73lzHy2/Veh2SMX/xJCiuAb4HvRKSxW1aaSXaMCQtJsVGMGtCeHi1r8vinq3n68x+tk5sp9/xJCqqqrwJDgY9F5DKc0VONMUWIjYrk5evbcl2HdF6ZsY5HpqwiN9f++5jyq8jOa7hXBao6W0TOA94DmgU0KmNCSGSE8MQfT6dSXDTDZ60n82gOT/VqRXSk9X425Y8/SeHivBequk1EumOjpBpTLCLCsB7NqBQfzdOfryHzaA4vXXcGcdGRXodmzP/wNfNaH1V9G7hO5KRNCLMCFpUxIUhEuK17E5LjonhkyipuemshI/q2Iz7GEoMpP3xdvya6P5MLWYwxJdDvzAY807s1s9fupv+o+Rw8cqzonYwJkkKvFFR1uPvzb8ELx5jw0CujLnHREdw1cSk3jJzHW4M6kJpgU3wa7/m6ffSirx1V9Y6yD8eY8HFpq9rERUVy6/jFXP/GPMYN7kDVpFivwzJhztfto0VFLMaYUjq/RRoj+7dj3a5MrntjLjsP2girxlu+bh+NDWYgxoSrs0+pzuiB7Rk8ZiHXDp/LOzd1omZKXNE7GhMART4oLSLVReQZEZkmIl/nLcEIzphw0blxNd4a3IEdB45w7Yg5bNt/2OuQTJjyd5Kd1dgkO8YEVPsGVXhrcEd2Z2ZzzfC5bNlnicEEn02yY0w5klG/MuMGd2DvoWyuGT6HzXsPeR2SCTM2yY4x5cwZ9Soz/saOHDh8jGtHzLXEYILKJtkxphxqVTeVty0xGA8UmRRU9ZOCk+yo6tRgBGdMOLPEYLzgz9NHDUXkORGZLCJT85bSHFREUkVkkoj8KCKrReRMEakiIl+KyM/uz8qlOYYxoSAvMew/fIzr3pjLVmt8NgHmz+2jj3CeOHoJeDbfUhovAJ+pajOgNc7TTcOA6araFJjurhsT9lrVTWXc4I7sy3ISw/b91sHNBI4/SeGIqr6oqjPcp4++UdVvSnpAt33ibOBNAFXNVtV9wOVAXoe5scAVJT2GMaGmTXoqYwd34LfMbK57Yy47DlhiMIHhT1J4QUQedW/xtM1bSnHMhsAuYLSILBGRkSKSCKSp6jZ3m+1AWimOYUzIaVuvMmMHtWfHgSNc/8Zcdh086nVIJgRJUXPGisi/gL7AOiDXLVZVLVFfBRFpB8wFuqjqPBF5ATgADFXV1Hzb7VXV37UriMgQYAhAWlpaxsSJE0sSBpmZmSQlJZVo34osHOsdanVes+c4zy48Qo0E4f4O8STHnHzK9FCrtz/Csc5Q/Hp37959kaq2O+mbqupzAdYCMUVt5+8C1AR+ybd+FvApsAao5ZbVAtYU9VkZGRlaUjNmzCjxvhVZONY7FOv83c+79JSHpmmP52fpvqzsk24TivUuSjjWWbX49QYWaiHfq/7cPloJpPqdgoqgqtuBTSJyqlt0HvADMBXo75b1B6aU1TGNCTVdmlRjRL92rN2ZSb/RNlGPKTv+JIVU4EcR+bysHkkFhgLjRWQ50AZ4Avg3cIGI/Ayc764bYwpxzinVefWGtqzasp9BYxZwKDvH65BMCCh06Ox8Hi3rg6rqUuBk97POK+tjGRPKzm+RxgvXnsHQCYu56a2FvNm/PXHRNuezKTmfSUFEIoHh6vQnMMaUQ5e0qsXRnNbc+/4ybh2/mNf7ZBAT5c9NAGN+z+dvjqoeB9aISL0gxWOMKYEr29bln1ecztc/7uSud5eQczy36J2MOQl/bh9VBlaJyHwgK69QVXsGLCpjTLFd37Eeh7JzePzT1cRFL+fS6r4fNzfmZPxJCn8NeBTGmDJx41mNOJR9nOe+/Il96VF076aInLwfgzEn488oqd8APwLJ7rJaSzHMhTEmsIae24Sbz2nE15tyeOrzNV6HYyoYf0ZJvRqYD/QGrgbmiUivQAdmjCkZEWHYRc3onh7FazPX8cqMtV6HZCoQf24fPQS0V9WdACJSHfgKmBTIwIwxJSci9G0RQ2q1NJ7+fA1JsVH079zA67BMBeBPUojISwiu3/Cv05sxxkMRIjzdqxVZR3N4dOoqkmKjuCqjrtdhmXLOny/3z9zezANEZADOOEXTAhuWMaYsREVG8OJ1Z9ClSVX+8sFyPl+13euQTDnnT0PzfcAIoJW7jFDV+wMdmDGmbMRFRzKibztOr5PC0HeW8N3Pu70OyZRjft0GUtUPVPUed/kw0EEZY8pWYmwUYwa2p1H1RIaMW8iSjXu9DsmUU/48fXSlO2/yfhE5ICIHReRAMIIzxpSd1IQY3hrUgerJsQwYvYA12w96HZIph/y5UngK6KmqKapaSVWTVbVSoAMzxpS9GpXieHtwR2KjIuj75jw27TnkdUimnPEnKexQ1dUBj8QYExTpVRJ4+8aOZB/P5YaR89hp8z2bfPxJCgtF5F0Ruc69lXSliFwZ8MiMMQFzSloyYwZ2YHfmUfqNms/+QzZJj3H4kxQqAYeAC4HL3OXSQAZljAm8NumpjOjbjvW7shg01ibpMY4iO6+p6sBgBGKMCb6uTavxwrVtuO2dxdzy9mLe6NfO5mIIc3b2jQlzPU6vxRN/PJ1vftrFn99fRm6uDbkdzvwZ5sIYE+Ku7VCPPYeyeeqzNVROiOaxnqfZkNthypKCMQaAW85pzN6sbN74dgOVE2O46/xTvA7JeKDQ20cicpmI1M+3/oiILBORqSLSMDjhGWOCRUR48OLm9Mqoy/Nf/cy4Ob94HZLxgK82hX8CuwBE5FKgDzAImAq8HvjQjDHBJiL8+8rTOb95DR6Zuoqpy7Z6HZIJMl9JQVU1r7vjlcCbqrpIVUcC1QMfmjHGC1GREbx8fVva16/Cve8tZdZPu7wOyQSRr6QgIpIkIhHAecD0fO/FBTYsY4yX4qIjeaN/O5rUSObmcYtsAL0w4ispPA8sBRbizMu8EEBEzgC2BTwyY4ynUuKjGTuoPdWTYxk0ZgFrd9oAeuGg0KSgqqOAc4DBwMX53toOWIc2Y8JAjeQ4xg3uQGREBP3enM/WfYe9DskEmK+nj9oCaYAAbUSkrVtWC6gWpPiMMR6rXzWRsYPac/BIDv1GzWdvVrbXIZkA8tVPYSGwEsibpil/TxYFzg1UUMaY8uW02im80b8d/UbNZ+CYBYy/sSOJsdbNKRT5alO4BzgAHAZGA5epand3sYRgTJjp1KgqL113Bss37+OW8YvJzsn1OiQTAL7aFJ5X1a7AUCAdmC4i74lIm2AFZ4wpX/5wWk3+deXpzLJxkkKWP6OkrheRKUA80Bc4BeepJGNMGLqmfT1+y3LGSaqSGMOjl7WwcZJCSKFJQUQaAdcClwObgInAE6paJo8fiEgkTrvFFlW91B06YyJQFVgE9FVVa9Eyphy65ZzG7D6YzajZG6ieHMtt3Zt4HZIpI76uFNYCy4EpOG0L9YBb8v4iUNXnSnnsO4HVOJP4ADwJ/EdVJ4rI6ziPwr5WymMYYwJARHj4kubsyTrK05+voWpiDNd2qOd1WKYM+Gpo/jvwIZALJAHJBZYSE5G6wCXASHddcJ5mmuRuMha4ojTHMMYEVkSE8HTv1pxzSnUe/HAFn63c7nVIpgyIavEbikQkUVWzSnxQkUnAv3CSy5+BAcBcVW3ivp8O/FdVW55k3yHAEIC0tLSMiRMnliiGzMxMkpKSSrRvRRaO9Q7HOkPw6n00R3lywRE2Hszl3ow4mleNDPgxC2Pn2j/du3dfpKrtTvqmqha6AHWAdkCMu14DeALY6mu/Ij7zUuBV93U34BOcznBr822TDqws6rMyMjK0pGbMmFHifSuycKx3ONZZNbj13pN5VM97dqae9shnumLzvqAdtyA71/4BFmoh36u+ejTfhfOU0UvAXBG5EacNIB7I8Dsl/V4XoKeI/ILTsHwu8AKQKiJ5bRx1gS2lOIYxJogqJ8bw1qAOJMdFMWD0fH7ZXeIbCcZjvtoUhgCnquqZOPf3XwYuVNW7VbXEA+Kp6gOqWldVG+A83fS1qt4AzAB6uZv1x2ngNsZUELVT4xk3uAM5uUq/UfPZeeCI1yGZEvCVFI6o6h4AVd0IrFHVRQGM5X7gHhFZi/NY6psBPJYxJgCa1Ehm9ID27Dp4lP6jF7D/8DGvQzLF5OuR1Loi8mK+9Vr511X1jtIeXFVnAjPd1+uBDqX9TGOMt86oV5nX+2Zw49gF3DR2IW8N7kBctHeNz6Z4fF0p3IfTiSxvKbhujDEndc4p1Xn26jYs+HUPQycsIee4jZNUURR6paCqY4MZiDEmtPRsXZu9Wdk8OnUVD364gievamXDYVQANvatMSZg+nduwG9Z2bw4/WeqJMYyrEczr0MyRbCkYIwJqLvPb8rerGxe/2YdVRKjGXJ2Y69DMj5YUjDGBJSI8FjP09h7KJsnpv1IakIMV7dL9zosUwhfDc0AiMgpIjJdRFa6661E5OHAh2aMCRWREcJzV7fhrKbVGPbBcj5fZeMklVdFJgXgDeAB4BiAqi7H6XRmjDF+i4mKYHjfDFrVTWXohCXMWfeb1yGZk/AnKSSo6vwCZTmBCMYYE9oSYqIYPaA99askcNNbC1mxeb/XIZkC/EkKu0WkMaAAItILKPEwF8aY8FY5MYZxgzuSEh9N/9HzWbcr0+uQTD7+JIXbgOFAMxHZAtwF3BLIoIwxoa1mShxv39iRCIG+I+exdV+ZTOhoykCRSUFV16vq+UB1oJmqdlXVXwIemTEmpDWslsjYQR04eCSHPm/O47fMo16HZPA9R3MfVX1bRO4pUA7OraQ9wFRV3RvYEI0xoeq02im8OaA9fd+cR//R85lwUyeS46K9Dius+bpSSHR/FpyGMxlnXuUM4L8Bjc4YE/I6NKzC630y+HHbQQaPXciRY8e9Dims+Rr7aLj782+FbSMifw9EUMaY8NK9WQ2eu6YNd05cwq3jFzO8bwbRkf40eZqyVmSPZhGJAwYDpwFxeeWqOkhVHwlgbMaYMNKzdW0OHjnGQx+u5J73lvH8NW2IjLAB9ILNn1Q8DqgJ/AH4BmeqzIOBDMoYE55u6FifYT2a8fGyrTz80cq8OdtNEPkz9lETVe0tIper6lgReQf4NtCBGWPC05/OaczBI8d4ZcY6KsVFMaxHMxtyO4j8SQp58+ntE5GWwHagRuBCMsaEuz9feCoHj+QwfNZ6kmKjGHpeU69DChv+JIURIlIZeBiYCiQBfw1oVMaYsCYiPHbZaWQezeHZL38iITaKwV0beh1WWPAnKUx3+yLMAhoBiIidHWNMQEVECE9d1YpDR4/zj09+ICk2kmva1/M6rJDnT0PzBycpm1TWgRhjTEFRkRG8cF0bzjmlOsMmr2DK0i1ehxTyfPVobobzGGqKiFyZ761K5Hs01RhjAik2KpLX+2QwYPR87nlvGbFRkVzUsqbXYYUsX1cKpwKXAqnAZfmWtsBNAY/MGGNc8TGRvDmgPa3qpjB0wmJmrtnpdUghy1eP5inAFBE5U1XnBDEmY4z5naTYKMYM7MD1b8zl5nGLGD2gPZ2bVPM6rJDjT5vCWhF5UERGiMiovCXgkRljTAEp8dG8NagD9asmMHjsQhb8ssfrkEKOP0lhCpACfAV8mm8xxpigq5oUy9s3dqRWShwDRy9gyUYbqLks+Tsd5/2q+p6qfpC3BDwyY4wpRI3kON65qRNVEmPoN2q+TetZhvxJCp+IyMUBj8QYY4qhZkocE4Z0IiU+mj5vzmPVVksMZcGfpHAnTmI4IiIHROSgiBwIdGDGGFOUOqnxTLipE4kxkfQZOY9NB3O9DqnC82c6zmRVjVDVOFWt5K5XCkZwxhhTlPQqCUwY0onYqEiemn+YNdttEOfSKDIpiKOPiPzVXU8XkQ6BD80YY/xTv2oiE4Z0IjJCuP6NuZYYSsGf20evAmcC17vrmcArJT2gm1RmiMgPIrJKRO50y6uIyJci8rP7s3JJj2GMCT8NqyUyrEOcJYZS8icpdFTV24AjAO7geDGlOGYOcK+qtgA6AbeJSAtgGM7ge02B6e66Mcb4rWZiBBPzXTH8uN2aP4vLn6RwTEQiAQUQkepAiVtzVHWbqi52Xx8EVgN1gMuBse5mY4ErSnoMY0z4alQ9iYlDOhEVKVz/xjx+2GqJoTikqOnuROQG4BqcMY/GAr2Ah1X1/VIfXKQBzpDcLYGNqprqlguwN2+9wD5DgCEAaWlpGRMnTizRsTMzM0lKSirRvhVZONY7HOsM4Vnv/HXekZXLkwuOcPS48pf2cdSvFOlxdIFT3HPdvXv3Rara7qRvqmqRC9AMuA24HWjuzz5+fGYSsAi40l3fV+D9vUV9RkZGhpbUjBkzSrxvRRaO9Q7HOquGZ70L1vnX3Vna+V/TtdVjn+vSjXs9iSkYinuugYVayPeqP08fdQK2qOorqvoysEVEOvqdkk7+mdE48zSMV9XJbvEOEanlvl8LsGEQjTGlUq9qAhOHdKJSfBR9Rs5joY2VVCR/2hRew3niKE+mW1Yi7q2hN4HVqvpcvremAv3d1/1xxlwyxphSSa+SwHs3n0m15Fj6jZrPnHW/eR1SueZPUhD3cgMAVc3Fv2k8C9MF6AucKyJL3eVi4N/ABSLyM3C+u26MMaVWKyWed4d0ok5qPANGz2eGzcdQKH+SwnoRuUNEot3lTmB9SQ+oqt+pqqhqK1Vt4y7TVPU3VT1PVZuq6vmqatd5xpgyU6NSHBOHdKJJjSSGvLWQaSu2eR1SueRPUvgT0BnYAmwGOuI+/WOMMRVJ1aRY3rmpE63rpnL7O4t5f+Emr0Mqd3zeBnL7J/xHVa8NUjzGGBNQKfHRvDW4AzePW8R9k5Zz8EgOg7o29DqscsPnlYKqHgfqi0hpejAbY0y5khATxcj+7bjotJr8/ZMfeO6LNeRrOg1r/jQYrwdmi8hUICuvsMCTQ8YYU6HERkXyyg1teXDyCl78ei17DmXzt54tiYwQr0PzlD9JYZ27RADJgQ3HGGOCJzJC+PdVp5OaGM3wb9azN+sYz13Tmtio0O39XJQik4Kq/g1ARJLc9UzfexhjTMUhIjzQoznVk2J5/NPV/JZ1lBH92lEpLtrr0DzhT4/mliKyBFgFrBKRRSJyWuBDM8aY4LnxrEY8f00bFv26l6tfn8OOA0e8DskT/jySOgK4R1Xrq2p94F7gjcCGZYwxwXfFGXUYNaA9m/Yc4o+vzA7LORn8SQqJqjojb0VVZwKJAYvIGGM8dFbT6rx785nk5Cq9Xvue2Wt3ex1SUPnbo/mvItLAXR6mFD2ajTGmvGtZJ4UPb+tC7dR4+o+aH1ad3PxJCoOA6sBknJFNq7llxhgTsuqkxvP+LWfSqVFV7pu0nKc++5Hc3NDvy1Do00ciEoczxEUTYAXOFJrHghWYMcZ4rVJcNKMHtueRKat4deY6NuzO4rmr2xAfE7qPrPq6UhgLtMNJCD2Ap4MSkTHGlCPRkRE88ceWPHxJcz5btZ3ew79n677DXocVML6SQgtV7aOqw3Gm4Dw7SDEZY0y5IiLceFYj3uzfjl93H6Lny9+xIEQn7PGVFE7cKlLVnCDEYowx5dq5zdL48LbOJMdFc/0bcxk/79eQGzPJV1JoLSIH3OUg0CrvtYgcCFaAxhhTnjSpkcxHt3ahc+NqPPThSu7/YDlHjh33OqwyU2hSUNVIVa3kLsmqGpXvdaVgBmmMMeVJSkI0owa0545zm/Dews30ev17Nu055HVYZcKfR1KNMcYUEBkh3HPhqYzs145ffzvEpS99x5c/7PA6rFKzpGCMMaVwfos0PhnalfQq8dz01kIe/+QHsnNyvQ6rxCwpGGNMKdWvmsgHt3Sm/5n1GfndBnoPn8Mvu7OK3rEcsqRgjDFlIDYqkr9d3pLXbmjLhl2ZXPzit7y3cFOFezrJkoIxxpShHqfX4rO7zqZV3RT+Mmk5t45fzG+ZR70Oy2+WFIwxpozVTo1n/I2dGNajGdNX7+TC/8zivyu2eR2WXywpGGNMAERGCH86pzEfD+1K7dR4bhm/mNveWczOg+V78h5LCsYYE0Cn1kxm8q2dufeCU/hy1Q7Of/Yb3pm3sdyOuGpJwRhjAiw6MoKh5zXlv3edRYvalXjwwxX0Hj6HFZv3ex3a71hSMMaYIGlcPYkJN3Xi6V6t+GV3Fj1f+Y77Jy1n18Hy0xBtScEYY4JIROjdLp0Z93Xjxq4N+WDxZro/M5Pnv/qJg0e8n7LGkoIxxnigUlw0D13Sgs/vPpuuTarx/Fc/c/ZTMxgxax1ZR70bmNqSgjHGeKhx9SRe75vBlNu60LJOCk9M+5EuT37Nc1+sYbcH/RvKVVIQkYtEZI2IrBWRYV7HY4wxwdI6PZVxgzvywS2d6dCgCi9+vZYu//6au99dyrz1vwWtZ3ShczQHm4hEAq8AFwCbgQUiMlVVf/A2MmOMCZ6M+pUZ0a8da3dmMub7DUxZspUPl2yhYbVEerSsyfkt0mhTN5WICAnI8cvTlUIHYK2qrlfVbGAicLnHMRljjCea1Eji8StOZ/5D5/Ns79bUrBTH8FnrufLV7+nwxFdMWbolIMeV8jJYk4j0Ai5S1Rvd9b5AR1W9vcB2Q4AhAGlpaRkTJ04s0fEyMzNJSkoqXdAVUDjWOxzrDOFZ71Cvc9YxZfmu4yzdmcO59aI5tUokUPx6d+/efZGqtjvZe+Xm9pG/VHUEMAKgXbt22q1btxJ9zsyZMynpvhVZONY7HOsM4VnvcKjzJScpK8t6l6fbR1uA9Hzrdd0yY4wxQVKeksICoKmINBSRGOBaYKrHMRljTFgpN7ePVDVHRG4HPgcigVGqusrjsIwxJqyUm6QAoKrTgGlex2GMMeGqPN0+MsYY4zFLCsYYY06wpGCMMeYESwrGGGNOKDc9mktCRHYBv5Zw92rA7jIMp6IIx3qHY50hPOsdjnWG4te7vqpWP9kbFToplIaILCysm3coC8d6h2OdITzrHY51hrKtt90+MsYYc4IlBWOMMSeEc1IY4XUAHgnHeodjnSE86x2OdYYyrHfYtikYY4z5vXC+UjDGGFOAJQVjjDEnhGVSEJGLRGSNiKwVkWFexxMIIpIuIjNE5AcRWSUid7rlVUTkSxH52f1Z2etYy5qIRIrIEhH5xF1vKCLz3PP9rjs0e0gRkVQRmSQiP4rIahE5M0zO9d3u7/dKEZkgInGhdr5FZJSI7BSRlfnKTnpuxfGiW/flItK2uMcLu6QgIpHAK0APoAVwnYi08DaqgMgB7lXVFkAn4Da3nsOA6araFJjuroeaO4HV+dafBP6jqk2AvcBgT6IKrBeAz1S1GdAap/4hfa5FpA5wB9BOVVviDLl/LaF3vscAFxUoK+zc9gCaussQ4LXiHizskgLQAVirqutVNRuYCFzucUxlTlW3qepi9/VBnC+JOjh1HetuNha4wpMAA0RE6uLMWDjSXRfgXGCSu0ko1jkFOBt4E0BVs1V1HyF+rl1RQLyIRAEJwDZC7Hyr6ixgT4Hiws7t5cBb6pgLpIpIreIcLxyTQh1gU771zW5ZyBKRBsAZwDwgTVW3uW9tB9K8iitAngf+AuS661WBfaqa466H4vluCOwCRru3zUaKSCIhfq5VdQvwDLARJxnsBxYR+ucbCj+3pf5+C8ekEFZEJAn4ALhLVQ/kf0+d55FD5plkEbkU2Kmqi7yOJciigLbAa6p6BpBFgVtFoXauAdz76JfjJMXaQCK/v80S8sr63IZjUtgCpOdbr+uWhRwRicZJCONVdbJbvCPvctL9udOr+AKgC9BTRH7BuS14Ls699lT39gKE5vneDGxW1Xnu+iScJBHK5xrgfGCDqu5S1WPAZJzfgVA/31D4uS3191s4JoUFQFP3CYUYnIapqR7HVObce+lvAqtV9bl8b00F+ruv+wNTgh1boKjqA6paV1Ub4JzXr1X1BmAG0MvdLKTqDKCq24FNInKqW3Qe8AMhfK5dG4FOIpLg/r7n1Tukz7ersHM7FejnPoXUCdif7zaTX8KyR7OIXIxz7zkSGKWq//Q2orInIl2Bb4EV/P/99Qdx2hXeA+rhDDt+taoWbMSq8ESkG/BnVb1URBrhXDlUAZYAfVT1qIfhlTkRaYPTuB4DrAcG4vzRF9LnWkT+BlyD87TdEuBGnHvoIXO+RWQC0A1neOwdwKPAR5zk3LrJ8WWc22iHgIGqurBYxwvHpGCMMebkwvH2kTHGmEJYUjDGGHOCJQVjjDEnWFIwxhhzgiUFY4wxJ1hSMKUmIsdFZKk7UuX7IpLgQQxX+DOwoYj8SUT6BSMmHzH0DNTovCIyTURSi7H9YyLy52JsP0ZENojIn9z1oe55n5Y3GqmIdBWR/+Tbp7H7+5FZjKoYj9gjqabURCRTVZPc1+OBRQU6zBW2X1S+MWpKG8MY4BNVnVTUtqHIfT5dVDW3yI3/d7/HgExVfcbP7ceQ799ZROYCnXH6wCwDPgE+A64r2Cci/++JKb/sSsGUtW+BJu547x+5Y7rPFZFWcOIv03EiMhsYJyJpIvKhiCxzl87udn1EZL77F+Zwd8hzRCRTRP7pbjvX3b8z0BN42t2+sYjcJCIL3O0+yLt6yf+XsbvdZyKySES+FZFmBSsjIue4n7nUHWwu2S2/z/385W4HKkSkgTjzGYwRkZ9EZLyInC8is8UZ976Du90AEXn5JMfK+7eZ425/U773CjveGhF5C1gJpIvILyJSzX3/Hvev+JUicle+z3rIje874NR85W3cf9Pl7jnxZ/4FAaJxRig9BvQB/htqneTCiqraYkupFpy/NMEZmG0KcAvwEvCoW34usNR9/RjOSJbx7vq7OIP1gdPDPAVoDnwMRLvlrwL93NcKXOa+fgp42H09BuiVL6aq+V4/DgzNd/w/u6+nA03d1x1xhsUoWLePgS7u6yS3jhfiTJQuOH9YfYIzdHUDnJ61p7vli4BR7naXAx+5nzMAePkkx3oM56/teJzeq5twBnrzdbxcoFO+z/jF3TcDpzd7ohv3KpyRcvPKE4BKwNp8/x7LgXPc138Hnj9JjAX/nfvi9Bp+G0gGvs47b4X9nthSvpe8QaOMKY14EVnqvv4WZ8ylecBVAKr6tYhUFZFK7jZTVfWw+/pcoJ+73XFgv4j0xfnyWuDcFSGe/x/wKxvnSxGcL90LComppYg8DqTifCl+nv9NcUaP7Qy87x4DIPYknzMbeM69LTZZVTeLyIU4X9RL3G2ScCY12YgzQNsK9xircCZCURFZgfMlXpQp7r/NYRGZgTP/R1cfx/tVnXHzC+oKfKiqWW4sk4GzcJLKh6p6yC2f6v5MAVJV9Rt3/7HA+0UFq6rjgHHuZzwCvAj0EKfdZhPORE/FuqVlvGVJwZSFw6raJn9Bvi/ak8kq4vMEGKuqD5zkvWPq/tkJHKfw3+ExwBWqukxEBuCMHZNfBM64+218BaKq/xaRT4GLgdki8gc3vn+p6vD/CdqZtyL/GDu5+dZzfcT6P4c8ybqv4xX1bxkUIlIb6KCqfxeRb3CS/cM4g9R96WlwplisTcEEyrfADXBicLrdWmA+B9d0nNtNeXMrp7hlvUSkhlteRUTqF3G8gzi3L/IkA9vEGT78hoIbu7FsEJHe7jFERFoX3E5EGqvqClV9EmeE3WY4Vx2D3KsNRKROXqxl4HJx5hmuipPIFpTweN8CV4gzgmgi8Ee3bJZbHu+2j1wGoKr7gb0icpa7f1/gm5N8bmH+ATzivo7HSWa5OLepTAViVwomUB4DRonIcpzRGvsXst2dwAgRGYzzl/8tqjpHRB4GvhCRCJwGzNtwRoMszETgDRG5A2fY5L/i3MLa5f7MnzDy/hq/AXjNPVa0+xnLCnzuXSLSHecLbhVOI+pREWkOzHGviDJxGliP+4jPX8txhn6uBvxDVbcCW4t7PFVdLM6TQvPdopGqugRARN5167kTJ+nk6Q+87jbK5420WiQROSPvmG7ROzjtFptw2n1MBWKPpJqwIiIvAYtVdbTXsRQkxXw81AtSikd/xR5JrRDs9pEJGyLyD5ynjEJuUqUg2g/8Q9zOa/5wH/1dijMXgCnn7ErBGGPMCXalYIwx5gRLCsYYY06wpGCMMeYESwrGGGNOsKRgjDHmhP8DAlLEPxzD6eYAAAAASUVORK5CYII=" + }, + "metadata": { + "needs_background": "light" + } + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 19, + "source": [ + "from scipy.optimize import newton, bisect\r\n", + "root = bisect(function, 0, 100, args=(0.10))\r\n", + "\r\n", + "for i in range(100):\r\n", + " root = bisect(function, 0, 100.0, args=(i / 100))\r\n", + " print(\"rms({:.4f}%)= {:.2f}%\".format(root, function(root) * 100))\r\n", + "\r\n", + "\r\n", + "print(\"rms(half_periode * {:.2f}%)= {:.4f}%\".\r\n", + " format(root, function(root) * 100))" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "rms(1/2 per * 100.0000%)= 0.00%\n", + "rms(1/2 per * 97.5220%)= 1.00%\n", + "rms(1/2 per * 96.0640%)= 2.00%\n", + "rms(1/2 per * 94.8385%)= 3.00%\n", + "rms(1/2 per * 93.7422%)= 4.00%\n", + "rms(1/2 per * 92.7319%)= 5.00%\n", + "rms(1/2 per * 91.7846%)= 6.00%\n", + "rms(1/2 per * 90.8861%)= 7.00%\n", + "rms(1/2 per * 90.0268%)= 8.00%\n", + "rms(1/2 per * 89.1999%)= 9.00%\n", + "rms(1/2 per * 88.4004%)= 10.00%\n", + "rms(1/2 per * 87.6242%)= 11.00%\n", + "rms(1/2 per * 86.8684%)= 12.00%\n", + "rms(1/2 per * 86.1305%)= 13.00%\n", + "rms(1/2 per * 85.4083%)= 14.00%\n", + "rms(1/2 per * 84.7001%)= 15.00%\n", + "rms(1/2 per * 84.0045%)= 16.00%\n", + "rms(1/2 per * 83.3202%)= 17.00%\n", + "rms(1/2 per * 82.6461%)= 18.00%\n", + "rms(1/2 per * 81.9813%)= 19.00%\n", + "rms(1/2 per * 81.3248%)= 20.00%\n", + "rms(1/2 per * 80.6759%)= 21.00%\n", + "rms(1/2 per * 80.0340%)= 22.00%\n", + "rms(1/2 per * 79.3984%)= 23.00%\n", + "rms(1/2 per * 78.7686%)= 24.00%\n", + "rms(1/2 per * 78.1440%)= 25.00%\n", + "rms(1/2 per * 77.5243%)= 26.00%\n", + "rms(1/2 per * 76.9088%)= 27.00%\n", + "rms(1/2 per * 76.2973%)= 28.00%\n", + "rms(1/2 per * 75.6894%)= 29.00%\n", + "rms(1/2 per * 75.0847%)= 30.00%\n", + "rms(1/2 per * 74.4829%)= 31.00%\n", + "rms(1/2 per * 73.8836%)= 32.00%\n", + "rms(1/2 per * 73.2866%)= 33.00%\n", + "rms(1/2 per * 72.6916%)= 34.00%\n", + "rms(1/2 per * 72.0983%)= 35.00%\n", + "rms(1/2 per * 71.5064%)= 36.00%\n", + "rms(1/2 per * 70.9157%)= 37.00%\n", + "rms(1/2 per * 70.3259%)= 38.00%\n", + "rms(1/2 per * 69.7368%)= 39.00%\n", + "rms(1/2 per * 69.1482%)= 40.00%\n", + "rms(1/2 per * 68.5598%)= 41.00%\n", + "rms(1/2 per * 67.9714%)= 42.00%\n", + "rms(1/2 per * 67.3828%)= 43.00%\n", + "rms(1/2 per * 66.7937%)= 44.00%\n", + "rms(1/2 per * 66.2039%)= 45.00%\n", + "rms(1/2 per * 65.6133%)= 46.00%\n", + "rms(1/2 per * 65.0215%)= 47.00%\n", + "rms(1/2 per * 64.4283%)= 48.00%\n", + "rms(1/2 per * 63.8336%)= 49.00%\n", + "rms(1/2 per * 63.2371%)= 50.00%\n", + "rms(1/2 per * 62.6385%)= 51.00%\n", + "rms(1/2 per * 62.0377%)= 52.00%\n", + "rms(1/2 per * 61.4343%)= 53.00%\n", + "rms(1/2 per * 60.8281%)= 54.00%\n", + "rms(1/2 per * 60.2189%)= 55.00%\n", + "rms(1/2 per * 59.6064%)= 56.00%\n", + "rms(1/2 per * 58.9903%)= 57.00%\n", + "rms(1/2 per * 58.3703%)= 58.00%\n", + "rms(1/2 per * 57.7461%)= 59.00%\n", + "rms(1/2 per * 57.1174%)= 60.00%\n", + "rms(1/2 per * 56.4839%)= 61.00%\n", + "rms(1/2 per * 55.8453%)= 62.00%\n", + "rms(1/2 per * 55.2010%)= 63.00%\n", + "rms(1/2 per * 54.5509%)= 64.00%\n", + "rms(1/2 per * 53.8944%)= 65.00%\n", + "rms(1/2 per * 53.2311%)= 66.00%\n", + "rms(1/2 per * 52.5605%)= 67.00%\n", + "rms(1/2 per * 51.8822%)= 68.00%\n", + "rms(1/2 per * 51.1956%)= 69.00%\n", + "rms(1/2 per * 50.5000%)= 70.00%\n", + "rms(1/2 per * 49.7950%)= 71.00%\n", + "rms(1/2 per * 49.0797%)= 72.00%\n", + "rms(1/2 per * 48.3535%)= 73.00%\n", + "rms(1/2 per * 47.6155%)= 74.00%\n", + "rms(1/2 per * 46.8649%)= 75.00%\n", + "rms(1/2 per * 46.1006%)= 76.00%\n", + "rms(1/2 per * 45.3215%)= 77.00%\n", + "rms(1/2 per * 44.5264%)= 78.00%\n", + "rms(1/2 per * 43.7139%)= 79.00%\n", + "rms(1/2 per * 42.8826%)= 80.00%\n", + "rms(1/2 per * 42.0306%)= 81.00%\n", + "rms(1/2 per * 41.1559%)= 82.00%\n", + "rms(1/2 per * 40.2563%)= 83.00%\n", + "rms(1/2 per * 39.3292%)= 84.00%\n", + "rms(1/2 per * 38.3713%)= 85.00%\n", + "rms(1/2 per * 37.3791%)= 86.00%\n", + "rms(1/2 per * 36.3482%)= 87.00%\n", + "rms(1/2 per * 35.2733%)= 88.00%\n", + "rms(1/2 per * 34.1480%)= 89.00%\n", + "rms(1/2 per * 32.9641%)= 90.00%\n", + "rms(1/2 per * 31.7113%)= 91.00%\n", + "rms(1/2 per * 30.3761%)= 92.00%\n", + "rms(1/2 per * 28.9402%)= 93.00%\n", + "rms(1/2 per * 27.3781%)= 94.00%\n", + "rms(1/2 per * 25.6521%)= 95.00%\n", + "rms(1/2 per * 23.7027%)= 96.00%\n", + "rms(1/2 per * 21.4258%)= 97.00%\n", + "rms(1/2 per * 18.6098%)= 98.00%\n", + "rms(1/2 per * 14.6679%)= 99.00%\n", + "rms(half_periode * 14.67%)= 99.0000%\n" + ] + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": null, + "source": [], + "outputs": [], + "metadata": {} + } + ], + "metadata": { + "orig_nbformat": 4, + "language_info": { + "name": "python", + "version": "3.9.1", + "mimetype": "text/x-python", + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "pygments_lexer": "ipython3", + "nbconvert_exporter": "python", + "file_extension": ".py" + }, + "kernelspec": { + "name": "python3", + "display_name": "Python 3.9.1 64-bit ('env': venv)" + }, + "interpreter": { + "hash": "a4046a09f199e48e61396275ea0e14da7d6ee49ffd95343aa82a9cc54eb636e7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} \ No newline at end of file