{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "---\n",
    "# TP1 Bilinéaire : Synthèse de filtres à partir du continu\n",
    "---\n",
    "\n",
    "|  [Readme du TP](../README.ipynb) | Sujet | Corrigé |\n",
    "|--------|-------|---------|\n",
    "|Python  | [sujet python](./tp1_bilineaire_python.ipynb) | [corrigé python](./tp1_bilineaire_python_corr.ipynb) |\n",
    "|Octave  | [sujet octave](./tp1_bilineaire.ipynb) | [corrigé octave](./tp1_bilineaire_corr.ipynb) |\n",
    "\n",
    "\n",
    "\n",
    "Allez voir le [Readme du TP](../README.ipynb) vous y trouverez le déroulé avec :\n",
    " - le travail préparatoire papier (TD)\n",
    " - le travail préparatoire Notebook (TD maison)\n",
    " - le script d'installation pour l'INSA,\n",
    " - des vidéos explicatives de la méthode bilinéaire.\n",
    " \n",
    "Après avoir exécuté le [script d'installation INSA](../../installation/README.md#etudiante-insa), la commande `go_signal` doit s'exécuter dans un terminal et lancer un **jupyter-lab avec kernel octave** pouvant exécuter du code Matlab/octave.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# But du TP\n",
    "\n",
    "On souhaite réaliser un filtre du second ordre résonnant avec la méthode de la transformation bilinéaire.\n",
    "Pour cela on rédige un script capable de faire cette synthèse pour n'importe quelle fréquence naturelle $\\omega_n$ et facteur d'amortissement $\\zeta$.\n",
    "\n",
    "$G(p) = \\frac{\\omega_n^2}{p^2+ 2\\zeta\\omega_n\\,p+{\\omega_n}^2}$\n",
    "\n",
    "![](https://i.stack.imgur.com/J0Ste.png)\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 1 - Vérification du calcul analytique\n",
    "---\n",
    "Dans le [travail préparatoire](../FILT_prepa_tp1.ipynb) vous avez pu trouver l'expression analytique des coefficients de la recurrence du filtre discret à partir de la transformée bilinéaire (ou trapezoidale ou tusin) de la forme :\n",
    "\n",
    "$$\n",
    "a_0.y[k]  =  a_1 .y[k-1] + a_2. y[k-2] + b_0. x[k] + b_1. x[k-1] + b_2. x[k-2]\n",
    "$$\n",
    "\n",
    "On rencontre dans la littérature (et les outils de calcul) plusieurs conventions de représentation :\n",
    "  - $G_d(z) = K'. \\frac{1 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}}$ ou bien $G_d(z) = K'. \\frac{b_2 z^{2}+ b_1 z + 1}{a_2 z^{2} +  a_1 z + 1}$\n",
    "  - $G_d(z) = K'. \\frac{1 + b_1 z^{-1} + b_2 z^{-2}}{1 - a_1 z^{-1} - a_2 z^{-2}}$ \n",
    "  - $G_d(z) = \\frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 - a_1 z^{-1} - a_2 z^{-2}}$ \n",
    "  - $G_d(z) = \\frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}}$ ou bien $G_d(z) = \\frac{b_2 z^{2}+ b_1 z + b_0}{a_2 z^{2} +  a_1 z + b_0}$\n",
    "\n",
    "> **Attention 1 ! En discret K' n'est PAS le gain statique**  \n",
    "> **Attention 2 ! Erreur de signe fréquente entre la récurrence et $G_g$**. La récurrence ci-dessus correspond à une seule de ces fonctions de transfert, sinon il faut modifier les coefficients ($-a_1$, diviser par $b_0$ ; mulitplier par K' etc.) \n",
    "\n",
    "\n",
    "Faites un script qui donne la valeur de ces coefficients en fonction de $\\omega_n$, $\\zeta$ et $T_e$ calculés avec votre formule théorique. On prendra la dernière représentation de $G_d$ qui est compatible Matlab/octave :  \n",
    "\n",
    "$G_d(z) = \\frac{b_2 z^{2}+ b_1 z + b_0}{a_2 z^{2} +  a_1 z + a_0}$ avec $a_2=1$\n",
    "\n",
    "\n",
    "\n",
    "Corrigez les nombreuses erreurs du code ci-dessous.  \n",
    "On prendra une fréquence propre $w_n=2000$ rad/s, $\\zeta=0.25$ et $F_e = 1$ kHz pour ce calcul : on a ainsi $\\omega_n.T_e=2$. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": [],
    "vscode": {
     "languageId": "json"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "num_math =\n",
      "\n",
      "   0.40000   0.80000   0.40000\n",
      "\n",
      "den_math =\n",
      "\n",
      "   1.00000   0.00000   0.60000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "clear all; close all; clc;\n",
    "%% ADAPTER et CORRIGER ce code d'init.\n",
    "%% Parametres du filtre continu désiré\n",
    "wn = 1 ;          % pulsation propre [rad/s]\n",
    "zeta = 0.707; % facteur d'amortissement [1]\n",
    "\n",
    "%% Paramètres du filtre discret imposés\n",
    "Fe = 48e3 ;       % Fréquence d'échantillonnage [Hz]\n",
    "Te = 0.001;     % Période d'échantillonnage   [s] \n",
    "\n",
    "%% MODIFIEZ avec vos formules les paramètres de récurrence\n",
    "% du genre \n",
    "% b2 = (1-4*zeta/(Te*wn)+4/(Te*wn)**2)\n",
    "K=1;\n",
    "b0 = 0.4 ;\n",
    "b1 = 0.8 ; \n",
    "b2 = 0.4 ; \n",
    "a0 = 0.6; \n",
    "a1 = 0; \n",
    "a2 = 1;\n",
    "\n",
    "\n",
    "%% On affiche les vecteurs de coefs \n",
    "num_math=[b2 b1 b0]\n",
    "den_math=[a2 a1 a0]\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pour vérifier\n",
    "\n",
    "On peut créer une **fonction anonyme** capable d'évaluer la fonction de transfert pour une valeur de z donnée.\n",
    "\n",
    "Par exemple `f = @(x,a) 3+x^2*a;` définit en une ligne la fonction $f: (x, a) \\mapsto 2+ a\\,x^2$ ce qui pourrait être fait de manière équivalente par\n",
    "```octave\n",
    "function val = f(x,a)\n",
    "    val = 2 + a * x^2\n",
    "end\n",
    "```\n",
    "\n",
    "On calcule ensuite avec cette fonction le gain statique et le gain HF :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gain_statique =  1\n",
      "gain_Hf = 0\n",
      "Il est possible qu'il y ait une erreur...\n"
     ]
    }
   ],
   "source": [
    "%% CORRIGEZ l'étourderie dans ce qui suit\n",
    "% On définit la fonction de transfert (où z peut être vecteur)\n",
    "Gd = @(z) (b0 + b1 .* z.^(-1) + b2 .* z.^(-2)) ./ (a0 + a1*z.^(-1) + a2*z.^(-2) ) ;\n",
    "% Indice pour l'erreur : convention Matlab !\n",
    "\n",
    "%% On vérifie avec le Gain statique\n",
    "gain_statique = Gd(1)\n",
    "%% Gain Hautes fréquences\n",
    "gain_Hf = Gd(-1)\n",
    "%% Plein de gains\n",
    "Gd(0:0.1:1);\n",
    "%% On vérrifie Gd(z=i) <-> f=Fe/4\n",
    "if abs(Gd(i)+2*i)<1e-10\n",
    "    disp(\"Vous semblez avoir corrigé l'erreur\");\n",
    "else\n",
    "    disp(\"Il est possible qu'il y ait une erreur...\");\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lancez une console avec `File->New Console For Notebook` et dedans tapez `who` ou tout autre commande à tester.  \n",
    "Par exemple lancez `Gd([1; i; -1])` pour voir ce que cela donne..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Verif avec la boite à outils \"control\"\n",
    "\n",
    "On va comparer ces valeurs avec celles trouvées numériquement par la fonction `c2d` : \"continuous 2 discrete\" avec l'argument `\"bilinear\"`.\n",
    "\n",
    "On peut accéder à la [doc de c2d](https://octave.sourceforge.io/control/function/@lti/c2d.html) en ligne (le help ne marche pas pour cette version)\n",
    "\n",
    "Ces valeurs doivent coller à 1ppm près. Sinon revenez avec un prof sur vos calculs et vérifications."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {
    "vscode": {
     "languageId": "json"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Transfer function 'G_de_p' from input 'u1' to output ...\n",
      "\n",
      "      6 s^2 + 6 s + 6\n",
      " y1:  ---------------\n",
      "      9 s^2 + 9 s + 9\n",
      "\n",
      "Continuous-time model.\n",
      "\n",
      "Transfer function 'Gd_c2d' from input 'u1' to output ...\n",
      "\n",
      " y1:  0.6667\n",
      "\n",
      "Sampling time: 1 s\n",
      "Discrete-time model.\n"
     ]
    }
   ],
   "source": [
    "%% Pour octave seulement il faut charger le package control\n",
    "%% Avec matlab inutile et commentez cette ligne\n",
    "pkg load control\n",
    "\n",
    "%% VOTRE CODE G_de_p = tf....\n",
    "% déclarant la fonction de transfert en continu\n",
    "% faire \"help tf\" dans une console si nécessaire \n",
    "G_de_p = tf([6,6,6],[9,9,9])\n",
    "\n",
    "%% VOTRE CODE Gd_c2d = c2d(...)\n",
    "% passant du continu au discret la fonction de transfert\n",
    "% c2d prend en argument une TF et en second argument Te\n",
    "% puis la méthode \"zoh\", \"bilinear\" etc.\n",
    "Gd_c2d = c2d(G_de_p,1,\"zoh\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On verifies que les vecteurs de coefficients trouvés avec c2d et ceux que vous avez analytiquement sont les mêmes.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Mauvais numérateurs !\n",
      "num_err =\n",
      "\n",
      "   0.26667  -0.13333   0.26667\n",
      "\n",
      "Mauvais dénominateur !\n",
      "den_err =\n",
      "\n",
      "   0.00000   1.00000   0.40000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%% On extrait les coefficients avec .num et .den\n",
    "num_c2d = Gd_c2d.num{1};\n",
    "den_c2d = Gd_c2d.den{1};\n",
    "\n",
    "\n",
    "%% On vérifie que cela colle\n",
    "disp(\"\")\n",
    "num_err = num_c2d-num_math;\n",
    "den_err = den_c2d-den_math;\n",
    "\n",
    "if (num_err * num_err' > 1e-6)\n",
    " disp('Mauvais numérateurs !')\n",
    " num_err\n",
    "else\n",
    " disp('Bon numérateur')\n",
    " erreur_de = num_err*num_err'\n",
    "end\n",
    "if (den_err * den_err' > 1e-6)\n",
    " disp('Mauvais dénominateur !')\n",
    " den_err\n",
    "else\n",
    " disp('Bon dénominateur')\n",
    " erreur_de = den_err*den_err'\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "On pourrait continuer d'utiliser les outils d'octave/Matlab notamment pour afficher \n",
    "  - la réponse harmonique avec `bode`\n",
    "  - la réponse indicielle avec `step`\n",
    "  - les poles et les zéros avec `pzmap`\n",
    "\n",
    "Mais ces outils sont parfois obscurs et inadaptés. Et dans notre cas n'aportent pas grand-chose.\n",
    "\n",
    "Nous allons afficher la réponse harmonique par nos propres moyens et pouvoir voir sa périodicité sous différentes échelles  \n",
    "**ce que ne permet pas la fonction bode !**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {
    "vscode": {
     "languageId": "json"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAJMmlDQ1BkZWZhdWx0X3JnYi5pY2MAAEiJlZVnUJNZF8fv8zzphUASQodQQ5EqJYCUEFoo0quoQOidUEVsiLgCK4qINEWQRQEXXJUia0UUC4uCAhZ0gywCyrpxFVFBWXDfGZ33HT+8/5l7z2/+c+bec8/5cAEgiINlwct7YlK6wNvJjhkYFMwE3yiMn5bC8fR0A9/VuxEArcR7ut/P+a4IEZFp/OW4uLxy+SmCdACg7GXWzEpPWeGjy0wPj//CZ1dYsFzgMt9Y4eh/eexLzr8s+pLj681dfhUKABwp+hsO/4b/c++KVDiC9NioyGymT3JUelaYIJKZttIJHpfL9BQkR8UmRH5T8P+V/B2lR2anr0RucsomQWx0TDrzfw41MjA0BF9n8cbrS48hRv9/z2dFX73kegDYcwAg+7564ZUAdO4CQPrRV09tua+UfAA67vAzBJn/eqiVDQ0IgALoQAYoAlWgCXSBETADlsAWOAAX4AF8QRDYAPggBiQCAcgCuWAHKABFYB84CKpALWgATaAVnAad4Dy4Aq6D2+AuGAaPgRBMgpdABN6BBQiCsBAZokEykBKkDulARhAbsoYcIDfIGwqCQqFoKAnKgHKhnVARVApVQXVQE/QLdA66At2EBqGH0Dg0A/0NfYQRmATTYQVYA9aH2TAHdoV94fVwNJwK58D58F64Aq6HT8Id8BX4NjwMC+GX8BwCECLCQJQRXYSNcBEPJBiJQgTIVqQQKUfqkVakG+lD7iFCZBb5gMKgaCgmShdliXJG+aH4qFTUVlQxqgp1AtWB6kXdQ42jRKjPaDJaHq2DtkDz0IHoaHQWugBdjm5Et6OvoYfRk+h3GAyGgWFhzDDOmCBMHGYzphhzGNOGuYwZxExg5rBYrAxWB2uF9cCGYdOxBdhK7EnsJewQdhL7HkfEKeGMcI64YFwSLg9XjmvGXcQN4aZwC3hxvDreAu+Bj8BvwpfgG/Dd+Dv4SfwCQYLAIlgRfAlxhB2ECkIr4RphjPCGSCSqEM2JXsRY4nZiBfEU8QZxnPiBRCVpk7ikEFIGaS/pOOky6SHpDZlM1iDbkoPJ6eS95CbyVfJT8nsxmpieGE8sQmybWLVYh9iQ2CsKnqJO4VA2UHIo5ZQzlDuUWXG8uIY4VzxMfKt4tfg58VHxOQmahKGEh0SiRLFEs8RNiWkqlqpBdaBGUPOpx6hXqRM0hKZK49L4tJ20Bto12iQdQ2fRefQ4ehH9Z/oAXSRJlTSW9JfMlqyWvCApZCAMDQaPkcAoYZxmjDA+SilIcaQipfZItUoNSc1Ly0nbSkdKF0q3SQ9Lf5RhyjjIxMvsl+mUeSKLktWW9ZLNkj0ie012Vo4uZynHlyuUOy33SB6W15b3lt8sf0y+X35OQVHBSSFFoVLhqsKsIkPRVjFOsUzxouKMEk3JWilWqUzpktILpiSTw0xgVjB7mSJleWVn5QzlOuUB5QUVloqfSp5Km8oTVYIqWzVKtUy1R1WkpqTmrpar1qL2SB2vzlaPUT+k3qc+r8HSCNDYrdGpMc2SZvFYOawW1pgmWdNGM1WzXvO+FkaLrRWvdVjrrjasbaIdo12tfUcH1jHVidU5rDO4Cr3KfFXSqvpVo7okXY5upm6L7rgeQ89NL0+vU++Vvpp+sP5+/T79zwYmBgkGDQaPDamGLoZ5ht2GfxtpG/GNqo3uryavdly9bXXX6tfGOsaRxkeMH5jQTNxNdpv0mHwyNTMVmLaazpipmYWa1ZiNsulsT3Yx+4Y52tzOfJv5efMPFqYW6RanLf6y1LWMt2y2nF7DWhO5pmHNhJWKVZhVnZXQmmkdan3UWmijbBNmU2/zzFbVNsK20XaKo8WJ45zkvLIzsBPYtdvNcy24W7iX7RF7J/tC+wEHqoOfQ5XDU0cVx2jHFkeRk4nTZqfLzmhnV+f9zqM8BR6f18QTuZi5bHHpdSW5+rhWuT5z03YTuHW7w+4u7gfcx9aqr01a2+kBPHgeBzyeeLI8Uz1/9cJ4eXpVez33NvTO9e7zofls9Gn2eedr51vi+9hP0y/Dr8ef4h/i3+Q/H2AfUBogDNQP3BJ4O0g2KDaoKxgb7B/cGDy3zmHdwXWTISYhBSEj61nrs9ff3CC7IWHDhY2UjWEbz4SiQwNCm0MXwzzC6sPmwnnhNeEiPpd/iP8ywjaiLGIm0iqyNHIqyiqqNGo62ir6QPRMjE1MecxsLDe2KvZ1nHNcbdx8vEf88filhICEtkRcYmjiuSRqUnxSb7JicnbyYIpOSkGKMNUi9WCqSOAqaEyD0tandaXTlz/F/gzNjF0Z45nWmdWZ77P8s85kS2QnZfdv0t60Z9NUjmPOT5tRm/mbe3KVc3fkjm/hbKnbCm0N39qzTXVb/rbJ7U7bT+wg7Ijf8VueQV5p3tudATu78xXyt+dP7HLa1VIgViAoGN1tubv2B9QPsT8M7Fm9p3LP58KIwltFBkXlRYvF/OJbPxr+WPHj0t6ovQMlpiVH9mH2Je0b2W+z/0SpRGlO6cQB9wMdZcyywrK3BzcevFluXF57iHAo45Cwwq2iq1Ktcl/lYlVM1XC1XXVbjXzNnpr5wxGHh47YHmmtVagtqv14NPbogzqnuo56jfryY5hjmceeN/g39P3E/qmpUbaxqPHT8aTjwhPeJ3qbzJqamuWbS1rgloyWmZMhJ+/+bP9zV6tua10bo63oFDiVcerFL6G/jJx2Pd1zhn2m9az62Zp2WnthB9SxqUPUGdMp7ArqGjzncq6n27K7/Ve9X4+fVz5ffUHyQslFwsX8i0uXci7NXU65PHsl+spEz8aex1cDr97v9eoduOZ67cZ1x+tX+zh9l25Y3Th/0+LmuVvsW523TW939Jv0t/9m8lv7gOlAxx2zO113ze92D64ZvDhkM3Tlnv296/d5928Prx0eHPEbeTAaMip8EPFg+mHCw9ePMh8tPN4+hh4rfCL+pPyp/NP637V+bxOaCi+M24/3P/N59niCP/Hyj7Q/Fifzn5Ofl08pTTVNG02fn3Gcufti3YvJlykvF2YL/pT4s+aV5quzf9n+1S8KFE2+Frxe+rv4jcyb42+N3/bMec49fZf4bmG+8L3M+xMf2B/6PgZ8nFrIWsQuVnzS+tT92fXz2FLi0tI/QiyQvpTNDAsAAAAJcEhZcwAACxMAAAsTAQCanBgAAAAddEVYdFNvZnR3YXJlAEdQTCBHaG9zdHNjcmlwdCA5LjI2WJButwAAIABJREFUeJzt3b9v22ieP/BH973imkj2FAecYRqID4gLqlpPCsnd2gdIxTbC2goOV8jASIvFFTrgpHEaBYeoiUF/CxWLQaQB7GJxCBWcmi1oIJ7pTBUbX0UWCbA+QDR2u3Go+QPyLZ4NvzxKpiVSj/gh/X5hMBAl8eETytRHzw8+n9Tnz58ZAABA1P4m6goAAAAwhoAEAABEICABAAAJCEgAAEACAhIAAJCAgAQAACQgIAEAAAkISAAAQAICEgAAkICABAAAJCAgwUOUctnc3Hz79u0sex0dHaVSKf/3XFxcuAvf3t6+urpynr+4uFhA7QESCgEJHqhqtfru3bvXr1+vrq7WajURhb969er6+np/f58xtr29/e7du+3t7cUeCCBJEJDggVpdXd3b26vVant7e58+feJPHh0dra6u8mYTb9mwLw2j1dXV6+trZ3f+ZCqVOjo6uqvwo6Oj4+Pj//mf/3n79u3V1dU//dM/XV1d3d7e7u3t8X339vZub28ZY91ulx+U14d9aVHx+kzuwl/lT66urh4fH/NqT60MQGx8Bnh4GGOPHz/e3d39xS9+wRh79erV58+f3717x740bn7xi1+srKx8/vz59evXjLFvv/223++vrKzwS+bVq1eMsX6/3+/3nd05Xsi3337r2eQPeLPJve/r16/fv3/Pj8uf2d3ddXZ8/Pjxt99+O7mLU1X+zMrKSr/f//Wvf40rGmLtb5cb/gCo2Nzc5B1ot7e3/X6/Vqvxb/lut8sYe/78+cHBwcXFxfv371dWVo6Pjxlj79696/V6/MHjx495IFlZWXn37t3sTZOjo6O9vb1+v89bYH/605/+9Kc/McZ4K+fXv/41bzNxx8fHvMfPs8vm5iZj7ODggDentre39/f3eX0A4gtddvBAbW9vHx8fHx8fd7vd//7v/3Y66Dxub29XV1f5Y+eBW61WOzg4uOsovJfv66+/dp45Pj7mm7/5zW/c7+SF80jjefKuXQASBgEJHqirq6ujo6OjoyPeIba5ufmP//iP7EsLiXfQbW9vf/3113wQ6Pr62pmMt729fXt7+/z58+fPn19fX7vbNJ7Cj46OHj9+zFs5HG9d8R35M/y4x8fH19fX/Ogek7sAJFPUfYYAEXBfAisrK69fv+bPV6tV53n+5E8//cTHmRhj/MHkkz/99JNTMu/3c+zu7r5//955/t27d3zUh7/Evow2uUtzjyG9e/fu8+fPk7u4X2Vfhp2+/fZbXNEQa6nPSGEO4MInwm1vb7s76C4uLjY3Nz39aRcXF6urqwFmcvNGlbPj9fW18wyfKT55u5JnF4BEQkACiNj19TWPf3t7e71e7/Xr1wu/LwogFhCQAKJ3e3v79u3b29vb/f19TzsM4OFAQAIAABIwyw4AAEhAQAIAABIQkAAAgAQEJAAAIAEBCQAASEBAAgAAEhCQAACABAQkAAAgIbH5kLa2tqKuAgAALR8+fIi6Cn4SG5AY+VO/HKkU0cU4llkxEcdaSJmBCwmw41y7zPhmsn9dy0T2JExWjP7PdHTZAQAACQhICUfztxtbbsVEHGshZQYuJMCOc+1C9s+GILLnimzFfCAgAQAACfEOSIZh2LYddS0AHqgnT55EXQVIlBgHpEqloqrq4eGhrutR14WutbW1qKsw3TIrJuJYCykzcCEBdpxrF7J/NgSRPVdkK+YjrrPsNE2TJKndbluW1Wq18vl81DUi6s9//nPUVZhumRUTcayFlBm4kAA7zrXL0j6dVOPH5RzoXp9PfhlsR1xiCxTXgGQYRjabZYxJkjQcDqOuDgBFwr/uf/Vd+EMEjgSLNfUfQqRuD0dcAxJjTJIk/iCXy019A50fXwBRefKH30ZdhXts/SHqGjDGGJs6GpZi33nfRv58enz8+NF5TH/ML8YBybIs/uCuFhJ+3TDG1tbWaLbcl1kxEcdaSJmBC5ljx5MPAY4145u3trYe2u3nkyHqH/7zX+JyidG/MTauASmbzRqGwRizLEuW5airQxfNS4VhDClcIYkZQ4qdyZ+5KfZ7T2cMkZ/CcfwQ4xqQisVir9dTFGU4HFar1airAwAP1JQQRTI+xUJcAxJjbDAY6Lr+7NkzZzAJACByngiE+DS7GAckxhhme98LY0iCjhWbMaRAu5D9syHo3nMVVXyK44cY74AEABAvaD/5ILpwengPcP4PwJLhKls4d3xaeHCi/3mhhQQAQIU7CAkNTjTFeC07mAXZ9aywll2YQrCWHR3iztXnk186/6UaP/L/KFRMHLSQAACoc1pITkxKZJsJY0gAEBCusgjxyDRXWKL/eaGFBAAQPzwUBQhLlCEgAQDEVcLCEiY1JBzZgU1MaghTCCY10EHhXLknPjhPUqjYvDCGBAAB4SojKNX48a6mEv3PK94tJMMwbNuOuhYAAFTwplLUtQgoTmNIT58+5ZkmZFluNpuVSkWSJNM0G40GFrUDAOB4TIrjkFJsWkg879HZ2dnZ2Vmz2dQ0TZKkdrvd6XS63W7UtaOLbD8yxpDCFIIxJDponquYtpPiFJAymUyr1VIUxbZtwzCy2SxjTJKkuzLGAiOcpAsJ+sIUggR9dJA9V3GMSbHpsrNte319fWdnxzTNer0uy7KTBimXy03dJZVKuTfp55MHiB36WbEfsic85/r/3f3rJvnvQNIBSdf1y8vLjY2NcrlcLBaLxSJjLJ/Pa5rGGLMsi7/trhZSUicQAhBBf9YWpBo/Ot+E9H89kO6yy+fzzWazXC4zxnq9nq7rzkvZbHY0GrEvY0uRVZE8mh3cDGNI4QrBGBIdZM8Vr1i8Ou5It5DccrlcvV4vFAqmafIGU6/XUxRlOBxWq9WoawcAAGHF7MZYXdclSXJGjzybbuhMABANV1lc8Fng9D+v2LSQOM/9Rrj9CAAgMYQEJPdgzyRZljOZjIjjwqS1tTWa01KXWTERx1pImYELCbDjXLuQ/bMhiOy5IlsxH0ICUqvVKhQKU18yTbNWq6FlAwCwNHxqA/VJ34ICUqFQaDabU19SVTWdTos4KExF9icSbowNUwhujKWD7LkiWzEfQgJSrVbr9XqMsYODg/Pzc03TMplMtVrNZrN8DjcAAICHkPuQ6vW6YRifPn0qlUqqqtZqtWw2++LFCxHHAgCAZBASkIbDYafTaTab+Xy+Wq3y/5umKeJY4I/4XXvxPRZujAWO7LkiWzEfYldqSKfTmFAXLbL9yBhDClMIxpDoIHuuyFbMB+mlgwAA4OEQdWOss4rf999/L+gQAACQJEJaSB/uMG85tm07S3pznpzlSGF+L7L9yBhDClMIxpDoIHuuyFbMh5AWkqIoU5+/6+aku/T7/U+fPjl7eXKWI4X5LMj2I2MMKUwhGEOig+y5IlsxHwLXsru5uTEMgy/ZYJrm1CVQfZRKJdM0v/nmG77p5Cy3LKvVatm27d5EQAIAiDshAYm3aUql0mAwcGbZlUoly7JmD0uDwcDd0vLkLJdlGSnMAQCSROAsO8/wjyRJnmfm5clZfm8Kc2CE+5ExhhSmEIwh0UH2XJGtmA+BXXYHBwf1er1cLmcymcvLS8uy7u1Yc+csn3zVnbNcluV7U5inUin3Jv188oKQzVu8zIqJONZCygxcSIAd59plxjeT/etaJrInYWtr6+PHj3/d+PcfIq3LTAQGpGazqarq+fm5bds7Ozunp6f37pLP5+8KWtls1jAM9iVnuWdz6i7xyj0IEDv0E76BIxaJzEXNsuPDSOVy2dPWUVXVGf6ZiydnOVKYAwAkjJAU5k+fPr2r1WJZVrvdDjwpDinM50U2SRcS9IUphEiCPlxlLFaXGP3PS0gLqdPp+Lx6V6yaBVKYAwAklZAWEgX0fwsAxB2usnih/3lhcVUAACABAQkAAEgQG5BUVVUURVVVTdOwCmokyN4chxtjwxSCG2PpIHuuyFbMh8CApCiKpmmMsdFodHNzc3h4KO5YcBea838YFlcNVwgWV6WD7LkiWzEfAgNSv9/vdDo7OzuMsWq1KsuyruviDgcAALEmMCCl02n3pmVZnmcAAAAcAgNSrVY7PDzUNM00zXq9zhgLsEADhES2HxljSGEKwRgSHWTPFdmK+RB7H5JhGHwY6a71UsWhP+MeIO5wlcUL/c9rSRljR6ORs8AdAADAJIFddjc3N+fn5/yxaZrj8XjeEmzb9k+hZBgGZpMDACQD3YyxjLF+v//p0yenXeWs2SrLcrPZrFQqkiSZptloNLCo3V1itPJjvI6FxVWBI3uuyFbMh8B8SFMzxs4ekEqlkmma33zzjVOaLMtnZ2d8U9M0SZLa7bZlWa1WCwEJACDuaGWMdRsMBu6xKMuyMplMq9VKp9O1Ws0wDD5nT5KkuzLGAiN8cxxujA1TCG6MpYPsuSJbMR+0Msb6sG17fX19Z2eHTyKXZdlpbOVyuam7IIU5gGhks3cD9/9TmMfhO1BgQGLTMsb603X98vJy6hzxYrFYLBYZY/l8nk8ld7oE72ohJTWzxlzI9iNjDClMIUTGkOhPI16CGF1i9H89CAxIhmF4Ztbd22WXz+fvek+v15Nl2Xk1m80ahsG+jC0tor4AABAlgQHJNM3RaMQfD4fDdDodZupBLper1+uFQsE0zXK5XCwWe72eoijD4bBarS6oyglE87cbwxhSuEIwhkQH2XNFtmI+lpcxtl6vv3z50pkFHoyu65IkOaNHnk03dCYAiIarLF7of17LS9C3vr5ummbIQvL5vDv8eDYBACC+BHbZKYrijkCmadZqNXGHg6liNOIar2NhUgNwZM8V2Yr5ENhl55nUIMtyyP66udBvnALEHa6yeKH/eQlsIWma5l5NVVVVWZaRgQIAAKYSEpAMwzg5OTFN09NlNxgMRBwOAAASQEhAymazZ2dnyDdBAdl+ZIwhhSkEY0h0kD1XZCvmQ8gYUqVSaTQaPFes+3lnadQloN9bChB3uMrihf7nJaSFVKvVJEkqFos7OzsiygcAgOQREpCcFRnC33gEAAAPhMAbY1+8eKGq6qWLuGPBXdbW1qKuwnTLrJiIYy2kzMCFBNhxrl3I/tkQRPZcka2YD4HTvs/Pz4n3VwIAAB0CW0gHBwd8Qe7ADMPwpJ01DMO27bs2YRLZaTZYXDVMIVhclQ6y54psxXyIzYd0eHjozg0x+yw727b5vjy7BJ8+XqlUJEkyTbPRaOTzec+miPoDAMDSCAxIPElEsH37/X4ul+NxaHd3t9lsapomSVK73bYsq9Vq2bbt3kRAugvZexFwH1KYQnAfEh1kzxXZivkQ20LyMAxjxqWDCoUCf+D0yDn7SpI0HA6dVYj4ppj6AgDA8ggMSMPhsNvt8i47PhQkSdJ4PJ5lASGeVELX9ZOTE2eNcCfTRC6Xm9yclEql3Jv088mL8OjRI5p5i5dZMRHHWkiZgQsJsONcu8z+Zpp/XctE/BL7+PGj8wz970CBAeny8vL09NRpEvEEfd1uV9O0u7rydF2/vLzc2Ngol8uKotzc3HQ6HSfwOBMceAvJvTm1tKXlHgR4mOjf+Q9uNAOnm9gU5u7seel02jTNlZWVm5ubu3bJ5/N8NEhVVR6NnJey2Syfs8enOXg2xf0rAABgOcRO+67X66qqaprWarV4fNI0zRkf8sEnfFe+YIwVi8XhcKgoSr1er1arnk1x/4q4I3tzHG6MDVMIboylg+y5IlsxHwIT9LEvXXCMsZWVlYODA56vL0zScV3XJUlySvBsuqEzAUA0XGXxQv/zEjvLLp1OO+urnpycFIvFkPOzPbtjtjcAQGIIDEiKogyHQ9u2nczljUZD3OEAACDWBI4h9fv9wWBQq9XK5fJgMEin0+KOBXch24+MMaQwhWAMiQ6y54psxXwIDEicJEmj0YgxJssyslEsH9lbtbGWXZhCsJYdHWTPFdmK+RA7y65UKsmyfH5+rihKv98PM50BAACSTWBAajabjUYjk8nw24nct7gCAAB4iO2y47Pgstlss9nEjLhIkO1HxhhSmEIwhkQH2XNFtmI+hMyyUxRl6vN89W4AAIBJQlpI33//PVbgJoLswCYmNYQpBJMa6CB7rshWzIeQFtJ//dd/aZrGl0DN5/OBsyIBAMDDIaSFxAeNBoNBuVw2DKNUKrVaLU3TRBwL/JHtR8YYUphCMIZEB9lzRbZiPsROauCRic+y+7d/+7d5d+dLrPq/wcngBwAAsSZw6SDLss7Pz3nDqFgsXlxczL6vbduHh4c86ZEsy3w2xNOnT3mmCf5MpVKRJMk0zUajgSl8dyHbj4wxpDCFYAyJDrLnimzFfAgJSL1ez4lDwW4/6vf7uVyOx6Hd3d1ms8kj09nZGX+DpmmSJLXbbcuyWq0WAhIAQNwJCUgnJyfr6+uSJF1eXvL0E5wTTu7l5ExyeuQsy8pkMq1WK51O12o1wzB4LlpJkjCjDwAgAYQEpNPT05Al8EaVrusnJye1Wo0xZtv2+vr6zs6OaZr1el2WZafhlcvlphaSSqXcm/TzyYvw888/P3r0KOpaTLHMiok41kLKDFxIgB3n2mX2N9PPii0a8Uvs48ePzjP0vwOFBKTAHWg8od/Gxka5XFYUhWcx54GnWCzy6eP5fJ73BzrzHe5qIQnNPQgA9BO+gRv9Xw9iE/TNK5/P82CmqiqPRs5LvV6P39XEN7PZrGEYjDE+thRJbQEAYIFoBSQHn/BdqVT45tnZWS6Xq9frhULBNM1yuVwsFnu9Hs8BWK1WI60sAAAsQCpe/Vq6rkuS5IweeTbd0JnAra2t0Zz9ucyKiTjWQsoMXEiAHefaZcY34ypjsbrE6H9eMQtIs6N/6gHiDldZvND/vIRnjAUAAJgFAhIAAJCAgJRwZBdYxOKqYQrB4qp0kD1XZCvmAwEJAABIwKQGAAgIV1m80P+80EICAAASEJASjmw/MsaQwhSCMSQ6yJ4rshXzgYAEAAAkYAwJAALCVRYv9D8vtJAAAIAE0gFJ13UnQR9nGIb7Gc8mAADEF9GAZNt2qVS6vLwslUo8+xFjrFKpqKp6eHio6/rkJkxFdmATkxrCFIJJDXSQPVdkK+aDaPqJfr9fLBar1eqzZ89arVaxWNQ0TZKkdrttWVar1bJt270ZOCVg4tFch5gtt2IijrWQMgMXEmDHuXYh+2dDENlzRbZiPogGJCfFkWEYPLuEYRjZbJYxJknScDiUZdm9GWFVAQBgIYgGJE5RlPPz81qtxjedvEe5XG5yc1IqlXJv0s8nDxA79LNiP3AfP350HtP/DqQVkHRdv7y83NjYKJfLjLFms1mr1fb29vimZVn8bbyF5N6cWlpSZ7TPJUbZw+J1LCToY3GYRrwEMbrE6P96oBWQ8vk8Hw3iw0LFYjGTyfCXstmsYRiMMcuyeH+dezPCOhNH81JhGEMKVwjGkOgge67IVswHrYDkKJfL9XrdMAzTNHmXXbFY7PV6iqIMh8NqterZjLq+AAAQFt2VGmzbNk1TkiRnrIgxpuu6+xnPphs6EwBEw1UWL/Q/L6L3ITHGMplMPp/3BBvPM5NvAA+y9yLgPqQwheA+JDrIniuyFfNBNyABAMCDQrfLLiT6jVOAuMNVFi/0Py+0kAAAgAQEpIQj24+MMaQwhWAMiQ6y54psxXwgIAFAQO5VAADCwxgSAASUSiX2CySR6H8rooWUcJ4F/ehYZsVEHGshZQYuJMCOc+1C9s+GILLnimzFfCAgAQAACYltcdNfRhAAYMmId9klNiABAEC8oMsOAABIQEACAAASiKafWCzbtk9OThhjPMdS1NUBSCB+lY3H42q1ms1mo64OxNKDaCH1+/18Pt9ut1VVjbouAMnU7XbL5XKn0+n1elHXBeLqQbSQCoUCz1IxHo+jrgtAMjWbTdu2NU1bX1+Pui4QV8lpIdm2bVmW+xnDMGzbZoxJkqTreqVSefnyZUS1A0gCz1XmXGKO8XiMn30QWHICUr/ff/PmjbNZqVRUVT08PNR1XdM0TdPOzs7QtQ0Qhvsqc19ijDFFUTKZTLlc9vwuBJhdQrrsSqWSaZrffPMN39Q0TZKkdrttWVar1cpkMrZtVyoVxtjZ2VmE9QSIL/dV5rnE8vl8Nput1+vr6+tI4gyBJSQgDQYDRVGcTcMweGNIkqThcEj85mSAWHBfZZ5LjDFWLBaz2axt2+iHgMASEpAmOT/TcrlctDUBSKTJS0ySJDSPIIzkjCF5OB3Z/OcbACwWLjFYuGQGpGw2OxqNGGOWZcmyHHV1AJIGlxiIkMwuu2Kx2Ov1FEUZDofVajXq6gAkDS4xECHJq33ruo5ObQBxcInBYiU5IAEAQIwkcwwJAABiBwEJAABIQEACAAASEJAAAIAEBCQAACABAQkAAEhAQAIAABIQkAAAgAQEJAAAIAEBCQAASEBAAgAAEhCQAACABAQkAAAgAQEJAABIQEACAAASEJAAAIAEBCQAACABAQkAAEhAQAIAABIQkAAAgAQEJAAAIAEBCQAASEBAAgAAEhCQAACABAQkAAAgAQEJAABIQEACAAASEJAAAIAEBCQAACABAQkAAEhAQAIAABIQkAAAgAQEJAAAIAEBCQAASEBAAgAAEhCQAACABAQkAAAgAQEJAABIQEACAAASEJAAAIAEBCQAACDhb6OugChbW1tRVwEAgJYPHz5EXQU/8Q5IhmFIkpTJZKa+SvzUL0cqlfr8+XPUtZhimRUTcayFlBm4kAA7zrXLjG8m+9e1TGRPwmTF6P9Mj3GXXaVSUVX18PBQ1/Wo6wIAAGHFtYWkaZokSe1227KsVquVz+ejrhFRNH+7seVWTMSxFlJm4EIC7DjXLmT/bJYslUot8G1wr7gGJMMwstksY0ySpOFwGHV1ACCZ4hub4xgmY9xlJ0kSf5DL5aKtCQAAhBfXFhJjzLIs/uCuFpLnB8KTJ0+E14men3/++dGjR1HXYoplVkzEsRZSZuBCAuw41y4zvvnJkyf0x8kfsq2trY8fPzqb9L8D4xqQstmsYRiMMcuyZFme/qZ//8G99XH6m+Ln88kvo64CAGOMbW1tJXsuaxx7vdw8nw79Xw9xDUjFYrHX6ymKMhwOq9Xq1Pck9Ys71fgx6ir4SeppBwDR4hqQGGODwUDX9WfPnjmDSQ8E8W/8ZcZL4qcCAOYS44DEGMNs73utra39+c9/XuYRZwwSC6nYYoPfXOFtIfUPXEiAHefaZfl/NgAs7gEJHrjFBr/5wts//979frTVYOGur6+vr68ZY9vb26urq55XLy4u9vb2PE9eXV2trq5ubm7O8maCiC56EV7ih1uBlHnbaskIYIm/yiJcFqjb7b569Wp/f58x9vbt27dv325vb/vU7fb2dm9vb3t7+/r6ent7+/j42P0G/pj+54UWEsACzBtgZgxgyYhbMK/r6+tXr17x5g5j7ODgoFarXV1dsS/Npsk2ULfb3dvbOz4+Zoxtbm7yB3e9mSwEpIQjOxiwzIqJOFbIMnmkubeQWeKWiEE7sn82D8Tr16/39/edbrrt7W0eja6urvb39/f39/mmG29LMcZub2+dJ/f29qa+mSwEJAC67go27oDhE7TQwIqvr7/+mj9wxn4uLi5ev359fHy8v79/e3v71Vdfud/Pm0EXFxdHR0fPnz/nT971ZrJiFpAMw8hkMs48b//0E8AYI/s7d5kVE3GshZQZuBD3jj5R53/FKtcsjHsDFdk/m8iJuKth8uP46quv3r9/zxs9FxcX7EtYur6+5s2myTkOjLGjo6Pr6+u3b986fXQ+b6YpNgHJtu3Dw0NZlvnSDM1ms1KpSJJkmmaj0cD8b4BJdwWeqd+qaE7NYjlnaX9/f29v7/nz5zyWXF1d/fDDD4yx7e1t3iPHZ9/xl/jz3W6XRyN3OZ430xebgNTv93O5XLPZZIzt7u5ms1mknwAIZuq36mSUQoiKyubm5vPnz7e3t/f29m5vb29vb3d3dxljv/nNb/b29t6/f39xcfH48eO3b9++f//+q6++6vf7t7e319fX7v49xtjR0ZHz5ij/PTOLzbRvvpSqJEm2bZdKpUKhsLGxUS6X2R1zT+lPcFwOsqPTD3xSQ8hClnNjrCdETcanxF9lkWeDvbi42NzcdE+Tu729vbq64ncm7e3t8cDDn/Hsm0qlfvrpJ+fNLA6fV2wCEqfr+snJSblcHo1GOzs7vGFUqVTOzs487/Qsc8visNItAGUff/Wd8/jJH34bYU2W5uPHj5S/Iff393kfXbfbrdVqnldTqdSTJ088q30TD0iku+x0Xb+8vHRaQoqi3NzcdDodSZIURbk3/QTlvySAWEux79ybSe3cI77a98HBwf7+vs9tRljte5Hy+bwzOKSqKo9GfHOm9BMAIAaPQE4X0Ozz92CB9vf3eU9djG599Uc6ILkZhmFZVqVS4ZtnZ2f3pp8AhjEkYcd6CGNIsxfuxCGs77dkiQlFXMzGkDx0XZckaWr6CfrDdwBx53+VJaDZFPmkhjAmK0//WzE2LaSpMNsbgCxPsym+YQmWJt4BCQDo46EoAQ0mEO1voq4AiLW2thZ1FaZbZsVEHGshZQYuJMCOc+0i4ox9Pvkl/y/V+HGZaYUhRtBCAoClcjeY6LeWiM/8TpgYD9n5oz98BxB34a+yuISlZKD/rRi/FpJlWel0mq/wjdW+AWItXq0lEC1mY0h8ITsegSqViqqqh4eHuq5HXS+6MIYk6FgYQ1ogZ2xpmQddFFxiCxSzFtLJyUk6nbZtm9+BhNW+ARLDiUloKj1YcWoh9Xq9jY0N3kdnGEY2m2WMSZJ011p2wAhnWkOCvjCFBNhxrl2i+rOJY1MJl9gCxSYgGYZhGIZ7lSBngYZcLhdRpQBg8TA1/MEi3WXnXu271+utr6/zRb5brVYqlbp3tW/PfE2knwBYOEELSPNrNcW+eyB5LsTxpJ+IsCazIB2Q3Kt9V6vV8XjMGBsOh8Vi8S9/+QvPy+uz2ne71lgKAAAUv0lEQVRSZ7TPBYurCjoWFldlS5lGnGLfER9SitElhvQTC8NHjBhj6XSaR6lSqYTVvu9F81JhGEMKV0hSx5Am8e47yjGJzrnyIFsxH/G+MRarfQNEaGlXGfGYFBf0vxVj00KaCrO9AR4CzAh/IGIzyw6CIXtzHG6MDVNIwm6MnQXZGeEEzxVHtmI+EJASjmw/MsaQwhTycMaQ6CN7rshWzIeQLjv/tXxkWcbqcwAQAP05DhCGkIDUarUKhcLUl0zTrNVqGPsBgGAQkxJMSEAqFArNZnPqS6qqptNpEQeFqWJ0k0S8joX7kCJEKiaRPVdkK+ZDyBjSXdGIMVYul507igKwLMtZoIExZhiGbduBSwOAmKI5wQFCEngfkmEYfG0FRzqdDhONWq0WY8yyrGKxWC6XK5WKJEmmaTYajck+QPoz7gHiLvKrjE47KRYi/7zuJfA+pOFw2O12+bo+vFkjSdJ4PB4MBgFK4xMl2u22bdsvXrxIp9NIPwEAkCQCp31fXl6enp6enZ2dnZ398MMP2Wy20+nkcjlN04KVtrGxoarq+fl5p9NB+okZkb0XAfchhSnkAd6HNBWFjjuy54psxXwIbCGZpule1CedTpumubKycnNzE6zAN2/e1Gq10WhUqVRkWb43/QRW++bIrqi4zIqJONZCygxcSIAd59plxjdH/tf1hMCi4JGfhLtsbW1hte+/Ojg4qNfrxWIxnU7rus4nfJ+cnHQ6nRlLcKefYIwVCgX+oFQqsS/dgOzu9BOxXqYPgD46YxL0FwWngGzgdAgMSM1mk0cUxtjGxkaj0RiPx51OZ+paqFO5009sbGyMRiP+2LbtbDZrGAbzTT8BAAAxInZxVd6I2djY4PcezR6KJpXLZZ5vgre0isVir9dD+gkAYMTuTILABE5qUBSFz18YjUY3NzeHh4chCxwMBjs7O+12m3fc8c1Op1MsFhdQ3YQiO7CJSQ1hCsGkhklRTXAge67IVsyHwBZSv9+/uLgwTfPy8rJarY5GI13XQ87P9uyO2d73InurNhZXDVMIFlelg+y5IlsxHwJbSJ4lgizLwqJBACAOhVngEIbAFlKtVjs8PJRl2bKser3OXGnIAQAAPAS2kMrl8suXL9PptCzL+Xz+7OxM3LHgLmT7kTGGFKYQjCHdZfmNJLLnimzFfAhZy05RlKnP+yy6unB07pAASCqyVxlm3E1F9vNyCGwh3dzcnJ+f88emaXoWWg3Gs7w3VvsGAEgMIWNIvCVUKpUGg4GTHLZUKlmWFfhWJNu2+YiUruu1Wu3e1b4B4MHCbUkxJbCF5E5cxBiTJMnzzFzOz89zuVy73R4MBpqmaZrGV/vudDrdbjd0ZROLbD8yxpDCFIIxJDrIniuyFfMhfC27crmcyWQuLy8tywrTjpFlWVVVviaeLMtY7RsAfKCRFEcCW0jNZrNYLJ6fn3e73ZWVldPT0zClSZKUTqd522hjY4O5FiK6a7VvYIRvjsONsWEKwY2xdJA9V2Qr5kNIC0lRFD6MVC6X+TI/DlVVZVme8YYk92rf3W6XJ4pljO3u7hYKhXtX+0b6CQDRKC8gTSEzReSQfoL1+33TNKe+ZFlWu92esRz3at8es6z2jfQTjLG1tTWaP5SWWTERx1pImYELCbDjXLvM+Gb604hTjR9F1zBGlxjlXw+ckIDkn/EoWLaIZ8+eVSqV0WhkmmahUMBq3wBwL4wkxYuQG2PF0XVdkiRn9Miz6Ub/txtA3MXiKkNActD/vAROahAhn8+7w49nEwDAAyuuxkjMAhIAACSV2ICkqqqiKKqqapqGNX4iQfbmONwYG6YQ3Bg7F6GNJLLnimzFfMQpYywEQHP+D8N9SOEKwX1IdJA9V2Qr5kNgQOr3+51OZ2dnhzFWrVb5MnTiDgcAcBeMJMUCMsYCAAAJAgMSzxiraZppmsEyxtq27VmPFekn5kW2HxljSGEKwRhSAIIaSWTPFdmK+RB7H5JhGHwYiS//M+/uvV7v06dPTlo/T74J//QT9GfcA8Rd7K6yB35PEv3PS+Bq34yxbDY7b6vIUSqVTNP85ptv+KaTb8KyrFarZdu2exP5kADgXli4gTiBXXaGYfBc5vV6fWtrq9VqzbX7YDBwohEvzZ1vAuknAAASRmAL6cWLF41GwzAMy7I+fPhQKpWcKBKMJ9/EveknsNo3Y+znn39+9OhR1LWYYpkVE3GshZQZuJAAO861y+xvpr9ep9evvltsnYlfYljt+69M08zn86qq8oCRy+XG47H/Lu58E5OvuvNNyLJ8b/qJeC3TJ8ja2hrNXuNlVkzEsRZSZuBCAuw41y4zvpn+mMRUKfbdAnvtYnSJ0f/1IDZjLF+N++XLl5qmnZ+f12o1/11mzzcxS/oJYIRvjsONsWEKwY2xdJA9V2Qr5kPgGFKj0VhZWalWq9lsdjwet9vtTCYTuLRisTgcDhVFqdfr1WrVs7nAagNAsuEmWbKET/t2uuk0TSsWiyGnwyH9xLxilD0sXsdCgj4W56tsgXPtYnSJ0f+8BHbZ8f4627adhlGj0QhZpieeYbY3AASA+d80CQxI/X7/j3/8o6qqjLFyuVypVMQdC+5C87cbwxhSuEIwhkQH2XNFtmI+hOdDkiRpNBoxxmRZNk1T9OEAAGaBkSSCBAakg4ODUqkky/L5+bmiKP1+H9ldAQDgLgIDUrPZbDQamUym0+kwxjqdDgLS8pFdYBGLq4YpBIurLsRCGklkzxXZivkQGJD4Oty6ro/HY54VKUAJk6t9u5/Bat/3ItuPjDGkMIVgDGlRwscksueKbMV8iF06iDG2vr7uPDPvpLh+v++s9m3b9uHhIV+gQZblZrPpv9o3AADEi8CAdH5+HmbOu2e1736/n8vleHDa3d3NZrNY7RsAwsMUcDoEdtkVCgW+uk8wntW+C4XCs2fP2JeeQKz2PSOy/cgYQwpTCMaQ6CB7rshWzIeQlRr4LUfj8diz0NzZ2dlc5fDsFU6CPsaYrusnJyflcnk0Gu3s7PCGUaVSmSx5a2vLvcwti8NKtwAQlY+/+u7JH34bdS0Wz7Pa90NcqeHeRVTv4r/at6IoNzc3fLaeoihY7RsgWvSXopldqvFjYv4td3mgq32n0+ler2fbNp99MPuOPqt9q6rKoxHfxGrfALBAGEmiQMgY0snJyfr6erlcHg6HfOmg8PiE78oXWO17RmT7kTGGFKYQjCGJEGwKONlzRbZiPoS0kIbDoTOoE2Zeg7t11W63Pa8OBgNd1589e4b7bQEAEkDsWnZhEiDNIp/PIxr5I3tzHG6MDVMIbowVJEAjiey5IlsxH8IXVwUAAJiFqIC0tbW1tbV1eHj4/fffb30h6Fjgg2w/MsaQwhSCMSRx5m0kkT1XZCvmQ8gYUuJnTwJAgmHGXVTEpjCPUJLukACgKdlXWfJiEv3PC2NIAABAAumANJl+gjFmWZaTcgLpJwBAEKSUXT7SAanf779588b9jG3bpVKJzyavVCqqqh4eHuq6HlEFY4DswCYmNYQpBJMalmOWmET2XJGtmA+6AalUKp2cnHiePDk5SafTtm1rmsbTT3Q6nW63G0kNY4HsvQi4DylMIbgPaWnujUlkzxXZivmgG5A86ScYY71eb2NjQ5KkTCaD9BMAsBzou1sagQn6FsswDMMwOp3O5eUlf8ZZoyGXy03dJZVKuTeRfgJg4R7K/YW/+i6m/1JP+okIazILWgHJJ/1Er9dbX1/nWSdarVYqlUL6iVmsra3RbLkvs2IijrWQMgMXEmDHuXaZ8c30pxEvUIp9N3UWeIwuMfoxlVZA8kk/Ua1Wx+MxY2w4HBaLxb/85S/X19cM6SfuQ/NSYRhDClcIxpCW7667ZcmeK7IV80ErIPngI0aMsXQ6zYNWqVRSFGU4HCL9BAAsAVZwEC3eKzXoui5J0tQFvx9UZwJAJB7mVRbfmET/86I7y24WSD9xL7L3IuA+pDCF4D6kCPF2kjPvjuy5IlsxH/EOSAAAy/f55JeYCy5CvLvsfNBvnALEHa6yeHXf0f+80EICAAgI7aTFQkBKOLL9yBhDClMIxpDo+If//Bf3kBIdcfwQSQekydW+LctyP4PVvgEgchhSWhTSY0i9Xu/Tp0/NZpNvtlotxphlWcVisVwuVyoVSZJM02w0GpO309LvLQWIO1xlHjwmkR1Vov950b0xtlQqmabprK/Kc0y0223btl+8eJFOp/lq33wlobvWdwAAWBoeioiHJcroBqTBYKAoirPJ17hTVZUx1ul0FEXBat8AQBDCUmB0A9KkN2/e1Gq10WhUqVRkWcZq37P4+eefHz16FHUtplhmxUQcayFlBi4kwI5z7TL7m+mv1ynaXeeKf92k2HeMsSd/+O1yK8XYl4phte/gfFb7ZowVCgX+fKlUYoxhtW+AaNEfkyCCh6XIW0v0fz3QCkg+q31vbGyMRiP+2LbtbDZrGAbDat8AQB468WZEKyD5KJfLfHlv0zRrtVqxWOz1eljtGwDiAmHpXqSnfU/yLO+N1b7vFaPsYfE6FhL0MVxljLGgH6Jz05K4yDQ1QR/xzytmAWl29E89QNzhKgtvCZHJQf/zik2XHQBA8jhxaJmRiSwEJACA6CEyMeJr2UF4ZBdYxOKqYQrB4qp0LPxc8ZXxnDSAgZfIi+OHiBYSAABFk20mlvRmEyY1AEBAuMqWL0xwov95kW4h2bY9Ho/ds7oNw5AkKZPJTN0EAEg2dxBKXsuJ9BhSv99/8+YNf2zbdqlUUlWV/58xVqlUVFU9PDzkC4HDVGT7kTGGFKYQjCHREeG5ckab3ANOTpSK44dIt4XkST9xfn6ey+WazaZt2/V6HeknAAAcU1pO//z7VOPHeLWc6AYkT/oJWZZVVdV13TRNWZYNw0D6iVnQXKaBLbdiIo61kDIDFxJgx7l2IftnQxDBczW1Wy8WkYluQPKQJCmdTmuaZppmuVwejUZIPwEQOfoLSD90X9JPpNgP9L8BaQUkn/QT3W6XZy5njO3u7hYKBaSfAIgW/Vlb4Lb1h6hrcB9akxry+Xyz2ZyaDMkjm83ybBRIP+GP7MAmJjWEKYTIpAZ38rcHC5fYAtFqIfl49uxZpVIZjUamaRYKBaSfmBHBDm4OY0hhCsEYEh1kzxXZivmI2Y2xSD8xr1SK6Ee8zIqJONZCygxcSIAd59plxjeT/etaJrInYbJi9L8VY9NC4jzTuzHbGwAgMYjG9vAw+QcAwIN4CymxAQkAAOKF1iw7AAB4sBCQAACAhJhNagjGtu2TkxPGWD6fLxaLUVcHIIH4VTYej6vVKl/WC2BeD6KF1O/38/l8u93my4QDwMJ1u91yudzpdHq9XtR1gbh6EC2kQqHA71Uaj8dR1wUgmfhK/Jqmra+vR10XiKvktJBs23ZWt+MMw7BtmzEmSZKu65VK5eXLlxHVDiAJPFeZc4k5xuMxfvZBYMkJSO5sfux/p+/TNE3TtLOzM3RtA4Thvso8GTIVRclkMuVy2fO7EGB2Cemy82Tz0zTNnb4vk8nYtl2pVBhjZ2dnEdYTIL7cV5nnEsvn89lstl6vr6+vT13KC2AWCQlInmx+nvR9xG9OBogF91U2mSGzWCxms1nbttEPAYElJCBNujd9HwCEMXmJ3bXSMcCMkjOG5HFv+j4ACAOXGCxcMgMS0vcBCIVLDERIZpcd0vcBCIVLDERI8mrfPun7ACA8XGKwWEkOSAAAECPJHEMCAIDYQUACAAASEJAAAIAEBCQAACABAQkehK2trYqLYRhR18hL1/WnT5+6V8C6l6Iozvsnd1QUZXd3l698ChALybwPCWAS/XV1Dw4Oms3mokpbYFEAy4EWEjxQhmGoqlqpVBRFUVV1d3fX3UBRFOXp06elUklVVZ5o2DAM51VVVXkby7OjYRitVqvVam1tbdXrdZ4riBe1u7urqqqiKE6TpdVqeZIJTVbM2X1q3W5ubvgzmqbxhU35oUulEsEmIMC90EKCh8L5Qt/Y2CiXy+Px+OTkpNPppNPpFy9eDAYDxli9XufhxzTNi4sLxlipVCoUCoyx8XhsmiYvYTQaSZLEI4d7R0mSzs/PT09P2+12vV7n73eK2tvb+9d//VdN0/L5vGEYpmlmMpnJejoVy+fzuq4Ph8M//vGPvCZ8L3fdeHpWXdcbjUa/30+n0x8+fOA5wLDqNsQOAhI8FDs7O/xBOp3mD2RZzufziqLIssyDhyzLo9Ho5uamXC7zaMGj0VSapnl2lCRJlmUeCXiouLy8LBaLvKjT09NMJvO73/2u3W5rmlYsFu8qmVeMMZbP509PT3VdtyzLtu3xeKzr+tS6ZTKZ9fV1TdNUVZVlGf11EEcISPBQ8K/4qSzLury85I83Njbuarvcu+PkG25ubpxAyANVPp/XNO38/Jw3rfwZhlGv1wuFwsbGBq+SbdueujmpiYrFoiRJPCZJktTpdGb5JwDQgTEkeOj4d32z2Ww2mysrK+l0WpZlJ8w43XRuPOHC5I6T78xms05RvM8tn8/3er1sNjtLzNM0rVAoNJvNcrnMn9nZ2fHUjb+HMaaqqmVZzWZzMBicn5/Pex4AIocWEjx05XJZ07RKpSJJkmmap6en+XyeTxlwUv4wxnjvnKIoNzc3PPZM7jgZvQ4ODg4PD1utFk/TkM1mJUl68eKFE2D87ezs1Ot15oo9jUbj8PCwXq+78z7w2CZJUqvVMgxjOBweHBws4twALBUWVwVgjDHDMMbjsbtbjy9l/ebNG/ZlCrVt26ZpyrLsbtxM7jhJ1/V0Os071mzbLpVKP/zww+R7Li8vJ8d+LMtyYo9pmvxA7mW2Lcty1tvmNXSOpSjKzs6Of90A6EBAAvDD5+Ytao6Aruuqqq6vr08WqOt6t9stFoszNp7upaqqpmm1Wg0BCeLi//zHf/xH1HUAoCuVSm1sbCw25U+tVpt8Mp1O8wP9/d///aIOlM1mZVn+u7/7u0UVCCAUWkgAAEACZtkBAAAJCEgAAEACAhIAAJCAgAQAACQgIAEAAAkISAAAQAICEgAAkICABAAAJCAgAQAACQhIAABAAgISAACQgIAEAAAkICABAAAJCEgAAEDC/wNLUyyIe0hGYgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bode(Gd_c2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {
    "vscode": {
     "languageId": "json"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAJMmlDQ1BkZWZhdWx0X3JnYi5pY2MAAEiJlZVnUJNZF8fv8zzphUASQodQQ5EqJYCUEFoo0quoQOidUEVsiLgCK4qINEWQRQEXXJUia0UUC4uCAhZ0gywCyrpxFVFBWXDfGZ33HT+8/5l7z2/+c+bec8/5cAEgiINlwct7YlK6wNvJjhkYFMwE3yiMn5bC8fR0A9/VuxEArcR7ut/P+a4IEZFp/OW4uLxy+SmCdACg7GXWzEpPWeGjy0wPj//CZ1dYsFzgMt9Y4eh/eexLzr8s+pLj681dfhUKABwp+hsO/4b/c++KVDiC9NioyGymT3JUelaYIJKZttIJHpfL9BQkR8UmRH5T8P+V/B2lR2anr0RucsomQWx0TDrzfw41MjA0BF9n8cbrS48hRv9/z2dFX73kegDYcwAg+7564ZUAdO4CQPrRV09tua+UfAA67vAzBJn/eqiVDQ0IgALoQAYoAlWgCXSBETADlsAWOAAX4AF8QRDYAPggBiQCAcgCuWAHKABFYB84CKpALWgATaAVnAad4Dy4Aq6D2+AuGAaPgRBMgpdABN6BBQiCsBAZokEykBKkDulARhAbsoYcIDfIGwqCQqFoKAnKgHKhnVARVApVQXVQE/QLdA66At2EBqGH0Dg0A/0NfYQRmATTYQVYA9aH2TAHdoV94fVwNJwK58D58F64Aq6HT8Id8BX4NjwMC+GX8BwCECLCQJQRXYSNcBEPJBiJQgTIVqQQKUfqkVakG+lD7iFCZBb5gMKgaCgmShdliXJG+aH4qFTUVlQxqgp1AtWB6kXdQ42jRKjPaDJaHq2DtkDz0IHoaHQWugBdjm5Et6OvoYfRk+h3GAyGgWFhzDDOmCBMHGYzphhzGNOGuYwZxExg5rBYrAxWB2uF9cCGYdOxBdhK7EnsJewQdhL7HkfEKeGMcI64YFwSLg9XjmvGXcQN4aZwC3hxvDreAu+Bj8BvwpfgG/Dd+Dv4SfwCQYLAIlgRfAlxhB2ECkIr4RphjPCGSCSqEM2JXsRY4nZiBfEU8QZxnPiBRCVpk7ikEFIGaS/pOOky6SHpDZlM1iDbkoPJ6eS95CbyVfJT8nsxmpieGE8sQmybWLVYh9iQ2CsKnqJO4VA2UHIo5ZQzlDuUWXG8uIY4VzxMfKt4tfg58VHxOQmahKGEh0SiRLFEs8RNiWkqlqpBdaBGUPOpx6hXqRM0hKZK49L4tJ20Bto12iQdQ2fRefQ4ehH9Z/oAXSRJlTSW9JfMlqyWvCApZCAMDQaPkcAoYZxmjDA+SilIcaQipfZItUoNSc1Ly0nbSkdKF0q3SQ9Lf5RhyjjIxMvsl+mUeSKLktWW9ZLNkj0ie012Vo4uZynHlyuUOy33SB6W15b3lt8sf0y+X35OQVHBSSFFoVLhqsKsIkPRVjFOsUzxouKMEk3JWilWqUzpktILpiSTw0xgVjB7mSJleWVn5QzlOuUB5QUVloqfSp5Km8oTVYIqWzVKtUy1R1WkpqTmrpar1qL2SB2vzlaPUT+k3qc+r8HSCNDYrdGpMc2SZvFYOawW1pgmWdNGM1WzXvO+FkaLrRWvdVjrrjasbaIdo12tfUcH1jHVidU5rDO4Cr3KfFXSqvpVo7okXY5upm6L7rgeQ89NL0+vU++Vvpp+sP5+/T79zwYmBgkGDQaPDamGLoZ5ht2GfxtpG/GNqo3uryavdly9bXXX6tfGOsaRxkeMH5jQTNxNdpv0mHwyNTMVmLaazpipmYWa1ZiNsulsT3Yx+4Y52tzOfJv5efMPFqYW6RanLf6y1LWMt2y2nF7DWhO5pmHNhJWKVZhVnZXQmmkdan3UWmijbBNmU2/zzFbVNsK20XaKo8WJ45zkvLIzsBPYtdvNcy24W7iX7RF7J/tC+wEHqoOfQ5XDU0cVx2jHFkeRk4nTZqfLzmhnV+f9zqM8BR6f18QTuZi5bHHpdSW5+rhWuT5z03YTuHW7w+4u7gfcx9aqr01a2+kBPHgeBzyeeLI8Uz1/9cJ4eXpVez33NvTO9e7zofls9Gn2eedr51vi+9hP0y/Dr8ef4h/i3+Q/H2AfUBogDNQP3BJ4O0g2KDaoKxgb7B/cGDy3zmHdwXWTISYhBSEj61nrs9ff3CC7IWHDhY2UjWEbz4SiQwNCm0MXwzzC6sPmwnnhNeEiPpd/iP8ywjaiLGIm0iqyNHIqyiqqNGo62ir6QPRMjE1MecxsLDe2KvZ1nHNcbdx8vEf88filhICEtkRcYmjiuSRqUnxSb7JicnbyYIpOSkGKMNUi9WCqSOAqaEyD0tandaXTlz/F/gzNjF0Z45nWmdWZ77P8s85kS2QnZfdv0t60Z9NUjmPOT5tRm/mbe3KVc3fkjm/hbKnbCm0N39qzTXVb/rbJ7U7bT+wg7Ijf8VueQV5p3tudATu78xXyt+dP7HLa1VIgViAoGN1tubv2B9QPsT8M7Fm9p3LP58KIwltFBkXlRYvF/OJbPxr+WPHj0t6ovQMlpiVH9mH2Je0b2W+z/0SpRGlO6cQB9wMdZcyywrK3BzcevFluXF57iHAo45Cwwq2iq1Ktcl/lYlVM1XC1XXVbjXzNnpr5wxGHh47YHmmtVagtqv14NPbogzqnuo56jfryY5hjmceeN/g39P3E/qmpUbaxqPHT8aTjwhPeJ3qbzJqamuWbS1rgloyWmZMhJ+/+bP9zV6tua10bo63oFDiVcerFL6G/jJx2Pd1zhn2m9az62Zp2WnthB9SxqUPUGdMp7ArqGjzncq6n27K7/Ve9X4+fVz5ffUHyQslFwsX8i0uXci7NXU65PHsl+spEz8aex1cDr97v9eoduOZ67cZ1x+tX+zh9l25Y3Th/0+LmuVvsW523TW939Jv0t/9m8lv7gOlAxx2zO113ze92D64ZvDhkM3Tlnv296/d5928Prx0eHPEbeTAaMip8EPFg+mHCw9ePMh8tPN4+hh4rfCL+pPyp/NP637V+bxOaCi+M24/3P/N59niCP/Hyj7Q/Fifzn5Ofl08pTTVNG02fn3Gcufti3YvJlykvF2YL/pT4s+aV5quzf9n+1S8KFE2+Frxe+rv4jcyb42+N3/bMec49fZf4bmG+8L3M+xMf2B/6PgZ8nFrIWsQuVnzS+tT92fXz2FLi0tI/QiyQvpTNDAsAAAAJcEhZcwAACxMAAAsTAQCanBgAAAAddEVYdFNvZnR3YXJlAEdQTCBHaG9zdHNjcmlwdCA5LjI2WJButwAAGLVJREFUeJzt3b9PW1miB/Dr96bYCsLUcaQwWlyYatgp8HQTIoViikVDGG2zIA2OXjUNCSk2q9HSJDEN1QhYaSNtMyYSWmmlJRIkXYw0wlvZBazWK+FUW4Qx/wCvuG8tHr+GmIt9MJ9PEV3OvT73cHD89Tn3+jh1cHAQAUCn/VenGwAAUSSQAAiEQAIgCAIJgCAIJACCIJAACIJAAiAIAgmAIAgkAIIgkLjWNjY2Ot0E4P8IJK6p2dnZVCp19+7dVCo1MjKyt7cXF7YcUan/b3Z2NtH2QvcTSFxHS0tLz58/f/r06fr6+tOnT1+/fv306dMoip4/f76+vt5ytdPT0+vr6+vr63fu3Hn+/LnhF3wQgcR19M9//jOKopGRkZGRkdnZ2adPn37yySfxmObly5flcjn6zxCqOdbZ2NhIpVLj4+OpVKqvr29pael4tX19fXGdDx48aBYeqWdvb29kZCQuiUdmp9U8Ozvb19eXSqX6+/vL5XJ8WD6fjx8bH3a8tuNnhCvjAK6fra2tGzduRFH06aefTk9PLy4uxoVRFE1PT79//z4eMK2srKysrERRFI+l4uPX19c//fTTGzduvH///nCdURTduXPn0aNHjx49un379qeffnpiPeesOS6Mh1xxYVzy1VdfxSXxf97z1N+ZLoYPZ4TEdTQ0NFQul58+fTo0NLSxsfHgwYN8Pj80NBRFUV9fX19f3/r6+u3bt7e2tuLoas7jPXv2bGRk5PHjxz/99FM8kDqsVquVy+Vyufyvf/1rZGTkxHr6+vqiKJqdnd3a2lpZWWkOYo7UHJ9xaWmpWfj3v/89iqIHDx7Eg7D4UcdrO63lEL6POt0A6ICXL19ubW09e/Ys/nFkZGR5efnEWbgoivL5/CeffHK4JI6B48bHx+M6Z2dnnz9/3qz/cD35fL6/v39lZSW+jvXo0aO7d+/+bM2nOV7b2S2HkBkhcR3VarU4MPb29srlcq1Wu3379uEDhoaG9vb2Hj9+/Pjx41qtFl+biaIongdbXFyMj/nZEx2vZ3Z2dnFxMT717du3a7XaiTXHQRJn5MrKyo0bN375y18er/94bae1HK6ATs8ZQge8f//+q6++av4viOe4Dg4O4gtL6+vr79+/j6/TRFEUXw1qXumJC+PLTodFUfTo0aN4Ox6pnFhP8/JVFEXxlNppNU9PTzdbuLi4GB+2vr7erP/g0MWwZm3Hz9i+boWLSR34CnOuq1qtVqvV+vr6mmOdeMA0NDQUT51tbGw0925sbNy9ezdOhf7+/v7+/vOf6HA9zZIoiuITnVFzuVze29trtueM+pu1nXZGCJ9AgnNpxkbzhoLwa4ar5b+/++67TrcBroZf/OIXd+/e/dD7DjpbM1whRkgABMFddgAEQSABEASBBEAQBBIAQRBIAARBIAEQhCQDqdFo1Ov1E3fV6/VGo5HguQDoMkkG0srKyg8//HC8vNFojI2N9fb2JnguALpMYl8/MTY2Vq1Wv/nmm+O75ufne3p6Go3GiZmUyWSSagMAIdje3m7hUYkF0urqaqFQOF6+vLx869ater1+xgiptabTslTKCh0doNs7Qre3X8vDjMv9gr5KpVKpVBYWFt6+fXvGYalUqrk9MDBwqU0iZmDaEbq9I3R7G+zs7DS3W34Zv9xAWl5evnnzZqFQqNfrT548yefz6XT6+GHev7RZKpUyKm0/3d4RV6LbD78pv4qOvIYHOkKanp7e39+Pomhzc3N0dLSnp+dSTwdwRV3d9+UJpullBVKpVJqammq+Menp6cnlcpd0Lj6UedGO0O0doduvkCQD6eHDh83tXC53eJj84sWLBE8EQPexUgMAQRBIAARBIAEQBIEEcK3VarWNjY2NjY29vb3jezc2No4XlsvlWq12zoPPTyABXF9LS0sjIyPr6+vr6+tDQ0PlcvnIAXfv3j38497e3tDQ0OLiYj6fn52dPfvgD3W5n0MCIFi1Wu3p06flcrmvry+Kovv37+fz+TiTarVarVbr7+8/8pA4wJ49exZFUX9/f7xx2sEfSiABXFOLi4vj4+NxGkVR1Bwhlcvl8fHx8fHx4wOm8fHxeOPw/N7IyMiJB38ogQRwff3qV7+KN0ZGRuKNjY2NxcXFZ8+ejY+P7+3tffzxx4ePj4dBGxsbs7Ozjx8/jgtPO/hDCSSA4KRm3iRe58H8F0dKPv74462trXjQE9+PEMdSrVaLh03NwdNhs7OztVrt5cuXzTm6Mw7+IALpBMk+FY4/CQDO1p7XjfHx8ZGRkcePH8dZUi6XX79+HUXR0NBQPCPXvJUuno4bGhpaWlqK0+hwPUcObplAOllSz4bLeJsDkIj+/v7Hjx8PDQ2NjIzs7e3t7e3duXMniqIHDx6MjIxsbW1tbGzcvn375cuXW1tbH3/88crKyt7eXq1WOzy/F0XR7Oxs8+CLtEcgAVxf+Xw+n89vbGz09/c3p+D6+/vL5XK5XI4HTyMjI3HwlMvloaGh45UcPvgijRFIANddc8TTFOdQczveOC2QDh98EQIJgLPcv39/fHz84h8z+lkCCYCzjI+PxwOjy84kgQTAz2jD8Ciylh0AgRBIAARBIAEQBIEEQBAEEgBBcJcdQOelUqlON6HzBBJAhx0cHHS6CUEwZQdAEAQSAEEQSAAEQSABEASBBEAQBBIAQRBIAARBIAEQBIEEQBAEEgBBEEgABEEgARAEgQRAEAQSAEEQSAAEQSABEASBBEAQBBIAQRBIAARBIAEQBIEEQBAEEgBBEEgABEEgARAEgQRAEAQSAEEQSAAEQSABEASBBEAQBBIAQRBIAARBIAEQBIEEQBAEEgBBEEgABEEgARCEJAOp0WjU6/Xj5ZVK5cRyAGj6KMG6VlZWfvrpp4cPHzZLGo3G1NRUNput1+vZbPbwLgA4LLER0tjY2Pz8/JHClZWV4eHhubm5Fy9evHr1KqlzAdB9Ehshra6uFgqFI4X37t2LNxqNRlInAqArJTlld1w6nY6iqFQqzc/P5/P50w5LpVLN7YGBgdbOtfPl96098ESZTCaZir78PrGqEhVmq7qebu8I3d4GOzs7ze2WX8YvN5CiKCoUCu/evVtYWIjD6UQHBwcXP1Fq5s3B/BcXr+f/zG8nUk1q5s32djJVJSiTyQTYqq6n2ztCt7dfy+8ALjeQisVinEaXehYAusBlfQ6pVCplMpn4hu/J/7ikcwHQBZIcIR2+qzuXyxkmA3B+VmoAIAgCCYAgCCQAgiCQAAiCQAIgCAIJgCAIJACCIJAACIJAAiAIAgmAIAgkAIIgkAAIgkACIAgCCYAgCCQAgiCQAAiCQAIgCAIJgCAIJACC8FGnG8DVlpp5k2BtB/NfJFgbcLUIJC4qqRRJNtuAK8eUHQBBEEgABEEgARAEgQRAEAQSAEEQSAAEQSABEASBBEAQBBIAQRBIAARBIAEQBIEEQBAEEgBBEEgABEEgARAEgQRAEAQSAEHwjbHX0c6X3/t6ViA0AumaSup7xwGSYsoOgCAIJACCIJAACIJAAiAIAgmAIAgkAIIgkAAIgs8hQZsk+3lknySj+wgkaJ+kUsRCG3QlU3YABEEgARAEgQRAEAQSAEEQSAAEwV12cBb3s0HbCCT4GUndq53JZKL57USqgq5kyg6AIBghAcFJcqb0y+8Tq4pL1o5AajQa+/v76XS6DecCuoNVLa6hdgTSysrKTz/99PDhwzacC+gUL/1c0KUH0tjYWLVa/eabby77REDHWfK1/ZJ9H9DZv+ClB9Lq6mqhULjsswBcIWGmSMfHuEHc1JBKpZrbAwMDLdby5feZTCaZBiVIq87vGrQqsaquQV8lJtRWDfz1f5KqLPPXhCq6QF/t7Ow0t1t+GQ8ikA4ODi5eSWrmzfZ2cB/y0Krz6/pWZTKZpKrq+r5KkFadX2rmzU4iNyW2mrVBBBLA5fG9iOeU2GfAWx2xCSSgmw389X8SHJgmUg+naUcgueEb6A4y6VIZIQGcS3fP14VAIMGV5LoI3UcgwdWTYISYgyIcVvsGIAgCCYAgmLKD686sHYEQSHCtuaOBcJiyAyAIAgmAIAgkAILgGhIB8WFPuM4EEqHwYU+45kzZARAEgQRAEAQSAEEQSAAEQSABEASBBEAQBBIAQRBIAARBIAEQBIEEQBAsHUR3snoQXDkCiS5kZVW4ikzZARAEgQRAEAQSAEEQSAAEQSABEASBBEAQBBIAQRBIAARBIAEQhKMrNXz77beNRuP4cS9evGhHcwC4ro4G0sLCwuTk5Oeff57NZjvSIACupxPWshsdHc1ms4ODg+1vDQDX1gmBNDExceKsHa1JcNlpa4YCXezk1b7n5+er1ero6Oj9+/d7e3vb3KZukmCE+D4FoLudHEhzc3OVSmVtbW1sbGxwcDCXy01MTLS5ZQBcK6fe9j04OPjw4cOZmZl6vf773//+zp07hUKhnS0D4Fo5OZBKpdKTJ08+++yzV69eTUxM/Pjjj6urq+/evSsWi21uHwDXxMlTdm/fvh0cHJyZmTl8AWliYiKdTrerYQBcLycH0sOHD48X5nK5S24MANeXpYMACIJAAiAIAgmAIAgkAIIgkAAIgkACIAgCCYAgCCQAgiCQAAiCQAIgCAIJgCAIJACCIJAACIJAAiAIAgmAIAgkAIKQZCBVKpVGo3H+cgBoSiyQJicni8Xi1NRUqVRqFjYajbGxsWKxGP+b1LkA6D7JBNLa2lo6nZ6bm1tYWFhaWmqWv3r1anh4eG5ubnV1dW1tLZFzAdCVPkqklkqlMjg4GEVROp3e3Nxslmez2WKxWCqVqtVqNptN5FwAdKVkAimKonQ6HW8MDw8fLuzp6VlbW6tWqxMTE6c9NpVKNbcHBgZabMGX32cymRYfeyUk+At2fV+FSrd3hG5vg52dneZ2yy/jiQVSvV6PNw6PkJaWlkZHR+MounPnzmmZdHBwcPEGpGbebG9vX7yeYCX4C3Z9X4Upk8no9vbT7e3X8juAZK4hDQ4O7u7uRlFUr9dNzQHQgmRGSKOjo8vLy4VCYXNzc3p6OoqiUqk0NTW1sbExOTm5u7tbrVbv3buXyLkA6EqJTdmtrq6WSqWvv/46vpiUy+XiYfLr168PlwPAiRILpCiKcrncB5UDQFOSgcRlS8286XQTAC6LQLoyDua/SKqqTCYTzbvvCAiLxVUBCIJAAiAIAgmAIAgkAIIgkAAIgkACIAgCCYAgCCQAgiCQAAiCQAIgCAIJgCAIJACCIJAACIJAAiAIAgmAIAgkAIIgkAAIgkACIAgCCYAgCCQAgiCQAAiCQAIgCAIJgCAIJACCIJAACIJAAiAIAgmAIAgkAIIgkAAIgkACIAgCCYAgCCQAgiCQAAiCQAIgCAIJgCAIJACCIJAACIJAAiAIAgmAIAgkAIIgkAAIgkACIAgCCYAgCCQAgiCQAAiCQAIgCAIJgCAIJACCIJAACIJAAiAIAgmAIAgkAIIgkAAIgkACIAgCCYAgJBlIlUql0WgcL6/X6/V6PcETAdB9PkqqosnJyXQ6Xa1WZ2Zmcrlcs/zJkydRFNXr9dHR0YmJiaROB0CXSWaEtLa2lk6n5+bmFhYWlpaWmuWlUimKorg83gaAEyUzQqpUKoODg1EUpdPpzc3NZvnbt29v3bpVLBajKFpYWEjkXAB0pcSuIaXT6XhjeHj4cPkPP/wQRdHu7u7k5GRS5wKg+yR2Dal528LhEVIURffu3YsvHY2NjZ322FQq1dweGBhosQVffp/JZFp87PWjrzpCt3eEbm+DnZ2d5nbLL+PJBNLg4GClUomiqF6vZ7PZZvmtW7d2d3fj7RNvwIsdHBxcvA2pmTfb29sXr+c6yGQy+qr9dHtH6Pb2a/kdQDKBNDo6ury8XCgUNjc3p6enoygqlUpTU1Pb29tjY2OFQqFarebz+UTOBUBXSmzKbnV1tVQqff311/HFpFwuF78rOVIOACdKLJCiKDr88aPzlANAk6WDAAiCQAIgCAIJgCAIJACCIJAACIJAAiAIAgmAIAgkAIIgkAAIgkACIAgCCYAgCCQAgiCQAAiCQAIgCAIJgCAIJACCIJAACIJAAiAIAgmAIAgkAIIgkAAIgkACIAgCCYAgCCQAgiCQAAiCQAIgCAIJgCAIJACCIJAACIJAAiAIAgmAIAgkAIIgkAAIgkACIAgCCYAgCCQAgiCQAAiCQAIgCAIJgCAIJACCIJAACIJAAiAIAgmAIAgkAIIgkAAIgkACIAgCCYAgCCQAgiCQAAiCQAIgCAIJgCAIJACCIJAACIJAAiAIAgmAIAgkAIIgkAAIgkACIAhJBlKlUmk0Gifuqtfrp+0CgCjBQJqcnCwWi1NTU6VS6ciuRqMxNjbW29ub1LkA6D7JBNLa2lo6nZ6bm1tYWFhaWjqyd35+vqenxwgJgDMkE0iVSmVwcDCKonQ6vbm5eXjX8vLyrVu30um0ERIAZ/goqYrS6XS8MTw83CysVCqVSmVhYeHt27dnPDY18yaRNmQymUTquQ70VUfo9o7Q7W2ws7PT3B4YGGitksQCqV6vxxuHR0jLy8s3b94sFAr1ev3Jkyf5fL6ZW4cdzH+RTCPmt5Opp9tlMpntbX3Vbrq9I3R7+7X8DiCZQBocHKxUKlEU1ev1bDbbLJ+ent7f34+iaHNzc3R0tKenJ5HTAdB9kgmk0dHR5eXlQqGwubk5PT0dRVGpVJqammq+Menp6cnlcomcC4CulNiU3erqaqlU+vrrr+NJuVwud3iY/OLFi6ROBEBXSiyQoigyBgKgZZYOAiAIAgmAIAgkAIIgkAAIgkACIAgCCYAgCKTr6PCqU7SNbu8I3X6FCCQAgiCQAAiCQAIgCAIJgCCkDg4OOtsC350F0GVa+w6qzgcSAESm7AAIhEACIAgCCYAgdDiQKpVKo9HobBu60hkde2SXP8FlOK1Xj5c3Go16vd6udnUnvd1BH/RS87Od/9/fffddgo37IJOTk/V6fXl5OZ1Ox198TiLO6Ngju44f+dlnn5VKpb/85S//+Mc/Pv/88w79Blfbaf1/Yvmf//znra0tXd2ylnvbU/3izvlS09PT85vf/Obf//73z/f2QYf87W9/+93vfndwcLC7u/vb3/62U83oPmd07JFdx4/0t7i40/r/xPJf//rXAwMDz58/70BDu0LLve2pfnHnf6lZWlpqPsm/+OKLM+r86NIz9BSVSmVwcDCKonQ6vbm52almdJ8zOvbIrmw2e+TIer3e29v75MmTnp6efD7f29vbid/gajut/08sX11dLRQKHWlnd2i5tz3VL+78LzVzc3Nx+c9eHejkNaTmEG94eLiDzeg+Z3TskV1Hfmw0Gjdv3hwdHb1x48a3337bpuZ2ndP63xP+MrTW257qiTjnS008oVcqlaampvL5/BkVdmyEFEVR8wKXEVKyzujYw7uy2eyRI0dHR0dHR6MoyuVya2trbWpu1zmt/z3hL0Nrve2pnohzvtREUVQoFN69e7ewsHD27QIdGyENDg7u7u5GUVSv17PZbKea0X3O6Ngju44fuby8XCqVOtHq7nFa/3vCX4aWe9tT/eLO/1JTLBbPk0ZRZ5cOGhsbGx4e3tzcnJ6ejt+tkIgjHRuPlOOlpY7sOvJjpVL59ttv7927V61WR0dHJyYmOv2rXEmHe7W3t/e0zo8Pjq9qPHz4sJMtvspa621P9USc86WmVCpVq9Wenp74US9evDitwg6vZVcqldzzfRnO6Ngju44f6Y9ycaf1ob69DC33tj/HxZ3/peY8LK4KQBAsHQRAEAQSAEEQSAAEQSABEASBBEAQBBIAQRBI8GEKhcLk5OTk5GQmk4k3isViqVSanJxsrcJSqfTZZ5+duMpqoVC4c+eONQW4JnwOCVpx+EPp0X+++S1e4biFqt6+fXvaYg2FQuHzzz/P5XKttxWuiE4urgpdo16vx2t0FovFer2+ubn5zTffRFH06tWrdDq9sLDQ29tbLBaXlpb29/fv379/Yvw0Go35+fmVlZVsNvuHP/yhtXiDq8uUHSRgf3+/Wq3u7++/evVqYWHhxx9//OMf/3jjxo3Xr1/39va+evWqUqkUi8XV1dWNjY1qtVosFo9XsrKy0tPTs729PTMzYwlqriEjJEhSNpttftvb9PR0FEU3b97c3d3d3d3NZrPVajU+Jl4L+YibN2+ura0Vi8VsNmu5Va4hgQRtUq/X3759G2/funXr+AGjo6PpdDrOpHiir70NhA4TSNAOt27devfuXTzuWV5ebi7Ff1ixWOzp6YmPyWQy7W4idJpAgnaYmJhYW1ubnJxMp9PVavVPf/rT8WPS6fSTJ08qlcrm5ub9+/fb30joLLd9Q/tUKpX9/f0j93Afvu270WjEX2XWvMXObd9cH+6yg/YZHBw8MVqa99319vbmcrlmGhWLxfg+CLgOjJCgw46PipriEdXhO/egi/0vKAFZb2IP1R0AAAAASUVORK5CYII=",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "step(Gd_c2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {
    "vscode": {
     "languageId": "json"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAJMmlDQ1BkZWZhdWx0X3JnYi5pY2MAAEiJlZVnUJNZF8fv8zzphUASQodQQ5EqJYCUEFoo0quoQOidUEVsiLgCK4qINEWQRQEXXJUia0UUC4uCAhZ0gywCyrpxFVFBWXDfGZ33HT+8/5l7z2/+c+bec8/5cAEgiINlwct7YlK6wNvJjhkYFMwE3yiMn5bC8fR0A9/VuxEArcR7ut/P+a4IEZFp/OW4uLxy+SmCdACg7GXWzEpPWeGjy0wPj//CZ1dYsFzgMt9Y4eh/eexLzr8s+pLj681dfhUKABwp+hsO/4b/c++KVDiC9NioyGymT3JUelaYIJKZttIJHpfL9BQkR8UmRH5T8P+V/B2lR2anr0RucsomQWx0TDrzfw41MjA0BF9n8cbrS48hRv9/z2dFX73kegDYcwAg+7564ZUAdO4CQPrRV09tua+UfAA67vAzBJn/eqiVDQ0IgALoQAYoAlWgCXSBETADlsAWOAAX4AF8QRDYAPggBiQCAcgCuWAHKABFYB84CKpALWgATaAVnAad4Dy4Aq6D2+AuGAaPgRBMgpdABN6BBQiCsBAZokEykBKkDulARhAbsoYcIDfIGwqCQqFoKAnKgHKhnVARVApVQXVQE/QLdA66At2EBqGH0Dg0A/0NfYQRmATTYQVYA9aH2TAHdoV94fVwNJwK58D58F64Aq6HT8Id8BX4NjwMC+GX8BwCECLCQJQRXYSNcBEPJBiJQgTIVqQQKUfqkVakG+lD7iFCZBb5gMKgaCgmShdliXJG+aH4qFTUVlQxqgp1AtWB6kXdQ42jRKjPaDJaHq2DtkDz0IHoaHQWugBdjm5Et6OvoYfRk+h3GAyGgWFhzDDOmCBMHGYzphhzGNOGuYwZxExg5rBYrAxWB2uF9cCGYdOxBdhK7EnsJewQdhL7HkfEKeGMcI64YFwSLg9XjmvGXcQN4aZwC3hxvDreAu+Bj8BvwpfgG/Dd+Dv4SfwCQYLAIlgRfAlxhB2ECkIr4RphjPCGSCSqEM2JXsRY4nZiBfEU8QZxnPiBRCVpk7ikEFIGaS/pOOky6SHpDZlM1iDbkoPJ6eS95CbyVfJT8nsxmpieGE8sQmybWLVYh9iQ2CsKnqJO4VA2UHIo5ZQzlDuUWXG8uIY4VzxMfKt4tfg58VHxOQmahKGEh0SiRLFEs8RNiWkqlqpBdaBGUPOpx6hXqRM0hKZK49L4tJ20Bto12iQdQ2fRefQ4ehH9Z/oAXSRJlTSW9JfMlqyWvCApZCAMDQaPkcAoYZxmjDA+SilIcaQipfZItUoNSc1Ly0nbSkdKF0q3SQ9Lf5RhyjjIxMvsl+mUeSKLktWW9ZLNkj0ie012Vo4uZynHlyuUOy33SB6W15b3lt8sf0y+X35OQVHBSSFFoVLhqsKsIkPRVjFOsUzxouKMEk3JWilWqUzpktILpiSTw0xgVjB7mSJleWVn5QzlOuUB5QUVloqfSp5Km8oTVYIqWzVKtUy1R1WkpqTmrpar1qL2SB2vzlaPUT+k3qc+r8HSCNDYrdGpMc2SZvFYOawW1pgmWdNGM1WzXvO+FkaLrRWvdVjrrjasbaIdo12tfUcH1jHVidU5rDO4Cr3KfFXSqvpVo7okXY5upm6L7rgeQ89NL0+vU++Vvpp+sP5+/T79zwYmBgkGDQaPDamGLoZ5ht2GfxtpG/GNqo3uryavdly9bXXX6tfGOsaRxkeMH5jQTNxNdpv0mHwyNTMVmLaazpipmYWa1ZiNsulsT3Yx+4Y52tzOfJv5efMPFqYW6RanLf6y1LWMt2y2nF7DWhO5pmHNhJWKVZhVnZXQmmkdan3UWmijbBNmU2/zzFbVNsK20XaKo8WJ45zkvLIzsBPYtdvNcy24W7iX7RF7J/tC+wEHqoOfQ5XDU0cVx2jHFkeRk4nTZqfLzmhnV+f9zqM8BR6f18QTuZi5bHHpdSW5+rhWuT5z03YTuHW7w+4u7gfcx9aqr01a2+kBPHgeBzyeeLI8Uz1/9cJ4eXpVez33NvTO9e7zofls9Gn2eedr51vi+9hP0y/Dr8ef4h/i3+Q/H2AfUBogDNQP3BJ4O0g2KDaoKxgb7B/cGDy3zmHdwXWTISYhBSEj61nrs9ff3CC7IWHDhY2UjWEbz4SiQwNCm0MXwzzC6sPmwnnhNeEiPpd/iP8ywjaiLGIm0iqyNHIqyiqqNGo62ir6QPRMjE1MecxsLDe2KvZ1nHNcbdx8vEf88filhICEtkRcYmjiuSRqUnxSb7JicnbyYIpOSkGKMNUi9WCqSOAqaEyD0tandaXTlz/F/gzNjF0Z45nWmdWZ77P8s85kS2QnZfdv0t60Z9NUjmPOT5tRm/mbe3KVc3fkjm/hbKnbCm0N39qzTXVb/rbJ7U7bT+wg7Ijf8VueQV5p3tudATu78xXyt+dP7HLa1VIgViAoGN1tubv2B9QPsT8M7Fm9p3LP58KIwltFBkXlRYvF/OJbPxr+WPHj0t6ovQMlpiVH9mH2Je0b2W+z/0SpRGlO6cQB9wMdZcyywrK3BzcevFluXF57iHAo45Cwwq2iq1Ktcl/lYlVM1XC1XXVbjXzNnpr5wxGHh47YHmmtVagtqv14NPbogzqnuo56jfryY5hjmceeN/g39P3E/qmpUbaxqPHT8aTjwhPeJ3qbzJqamuWbS1rgloyWmZMhJ+/+bP9zV6tua10bo63oFDiVcerFL6G/jJx2Pd1zhn2m9az62Zp2WnthB9SxqUPUGdMp7ArqGjzncq6n27K7/Ve9X4+fVz5ffUHyQslFwsX8i0uXci7NXU65PHsl+spEz8aex1cDr97v9eoduOZ67cZ1x+tX+zh9l25Y3Th/0+LmuVvsW523TW939Jv0t/9m8lv7gOlAxx2zO113ze92D64ZvDhkM3Tlnv296/d5928Prx0eHPEbeTAaMip8EPFg+mHCw9ePMh8tPN4+hh4rfCL+pPyp/NP637V+bxOaCi+M24/3P/N59niCP/Hyj7Q/Fifzn5Ofl08pTTVNG02fn3Gcufti3YvJlykvF2YL/pT4s+aV5quzf9n+1S8KFE2+Frxe+rv4jcyb42+N3/bMec49fZf4bmG+8L3M+xMf2B/6PgZ8nFrIWsQuVnzS+tT92fXz2FLi0tI/QiyQvpTNDAsAAAAJcEhZcwAACxMAAAsTAQCanBgAAAAddEVYdFNvZnR3YXJlAEdQTCBHaG9zdHNjcmlwdCA5LjI2WJButwAAIABJREFUeJzt3T9sG1fa7/GHixhbLEBKqjUGzAXMgqyieAHSuE1MA2JxEUAbSS5pIJLxViwuZamhC7uJRAELVwtLAWxgm4jCChu8xQgQ5U4kEEuuyEIBosJMu2aobpFCt5hkwqWGFP/MnzMz309Fj0czR0fD+fE5czgTubq6EgAAvPYHrxsAAIAIgQQAUASBBABQAoEEAFACgQQAUAKBBABQAoEEAFACgQQAUAKBBABQAoEEAFACgQRItVqNdInH42dnZwPWrFarA7aWzWYj/23w+sM3b39/X0Ta7bbxz/X19Uk2C6iGQAJ+tbKycnR09PXXX7fb7dXV1bG3s7m5efSbTz/9VESmp6dtaeHp6amInJ2dTU1N2bJBQCkEEvCr6enpbDa7vr4+Nzf3/v17EVlfX5+enu5XM62vr1tWKnNzc9lsNpvNnp2dvX///unTp3Nzc9fXN+oeYxfmBgfs7ssvvzQWHh0dGRsUkXa7bRZk2Wy23W6bW15cXIxEItPT0zs7O7b3FeAEAgn4Vbvdrlar+/v7Z2dnn376abVa3draWlxcPDo6MrKqe+XNzc2tra1KpVKpVLa2tjY3N3u2Vq1WNzY2vvzyS+O/+q2/v79vVGODdyci8Xj8+PhYRM7OzsxA2tnZOT4+NjZ7fHxsjOkZLi4ujo6O4vH4+vq6EVSA4ggk4Fe7u7sPHz5cWloSkZ2dnaOjI+NFNpvd2Nj4+eefuy8FHR0d3blz5/T09PT0dGpq6ujoaH193SyM2u324uLinTt3zOrk+vrG8s3NTSOcBu9ORB4+fCgi1Wr17OzMeC0i6+vrxjZfvXolIj/++KO5/ubmprmpfpfEAKUQSMCvjGtIp6en7XbbLEGGsbq6asSYyahv9vf3La8eda8/0uWlBw8eVCqVn3/+2ayfNjc3P/vsMxF58uSJ5Y/YdfkKcAGBBPzKGCgzo+jPf/6ziBglTqVSmZqa6k6pubm5dru9sbGxsbFxcXHRbrc3Nzer1Wq1Wq1UKu/fv4/H45VKZX19fX19/eLi4vr6PXsfvDsRMWJyf3/fmChhMAovY5s9G6xUKiJiVE4j5SvglU+8bgCgqNXV1dPT0ydPnhjFx6tXr7qrjY2NjWq1OjMzIyKffvpp98QBY3zs/fv3xswIEXn48OH19XuG0QbvTkSmp6c/++yzra2t7stLT548WVpampmZefDgQU/7z87OIpGI5aYANUV4hDkwgHFBaG5uzvKcXq1Wp6enh68/blx/8O4sGfVW9zar1erDhw+Ni1LxeDwejw+5KcBbBBIQNGYgXZ+qB6iMa0hA0MTj8adPn1IYwXeokAAASqBCAgAogUACACiBQAIAKMGbQOp0Oq1Wy5NdAwDU5E0gVSqVb7/91pNdAwDU5EEgLSwsbG9vu79fAIDKPAikg4ODr776yv39AgBUxqQGAIASVLy5aiKR8LoJAIDxnZ+fj/FTKgaSjPvLhEokwl02hhK8jvr4Xd97As180fsQiiFN0ktOtEdZwTucnDB2UaFoIAGQPud61c7yA9rji/ZDHd4E0tramif7BRTXcwb3+7nbsv3XU8rvv+bkjCdX+ZddVSMVEuClgCXQMK7/jiHshOv8OxJoY5oSSH519+5dr5vgD6p1lJonX297qacTurtIkf4xqXY4BQyBBLjBPMmqdoZVUHcXqRxOsB2BBDiFk+nkCKdQIZAAO3HSdI5lONHJQUIgATbg/Ogys5/p+SAhkIDxcTb0HMk0uYuLi4uLCxGZm5ubnp7u+d9qtZrNZnsWnp2dTU9Px+O9M/gtVx4egQSMjHOfgkim8ezs7Hz99deLi4sisrq6ur+/Pzc3173Cw4cPu6ekt9vtbDY7Nzd3cXExNze3ubk5YOVREUjAsDjT+UJIkilSfHu1/fngJTe6uLj4+uuvjXJHRJaWllZXV8/OzuS3sul6DbSzs5PNZo0cisfjxot+K4+KQAJuZpzaAnleC7BgJ9PV9ufdCTRGGonIq1evFhcXzWG6ubk5I43Ozs4WFxcXFxeNf3YzaikRabfb5sJsNmu58qgIJKCvQJ7IQqgnmQLz1zQzabw0Mnz22WfGC/PaT7VaffXq1ebm5uLiYrvdnpmZ6V7fKIOq1er6+vrGxoaxsN/KoyKQAAsBO3PBYPxBg/Q5Y8I0mpmZOT09NYqearUqv8XSxcWFUTZdn+MgIuvr6xcXF/v7++YY3YCVR0IgAb8L0qkK/QRpKG/CCmlxcTGbzW5sbBhZcnZ2dnx8LCJzc3PGiJwx+874L2P5zs6OkUbd2+lZeWwEEiBCSRRKfh/KM3No7EyKx+MbGxtzc3PZbLbdbrfb7QcPHojIkydPstns6elptVq9c+fO/v7+6enpzMxMpVJpt9sXFxfd43sisr6+bq480a90pZ67d+963QQfoJeGdGNH/ftfd/79rzvuNEZZHE5Xwx0JDnXUGKdi+X/HNy4Z3tHR0Y8//ti95OPHj0dHRx8/fry6unrw4IGx8PT01KIlIt0rX03QS1RICKkADNfAXt1XmNQ/Kq7XQ2NfSZKuGQ2m6elpc6F5cejs7KznW0rXV54EgYTQ8csZB57wVyy5Y2lpaXFxcfKvGd2IQEKIcJbBkIilbouLi0Zh5HQmEUgIBc4sGAOxZHKhPBICCYFX3/rl43fxkJ9NMIn//vbSLW8bE2wEEgLLOIOkn946Pz/3ui3wPSOW6hLn841z/uB1AwBHGGcNThywV/rprZkvLj5+F+9+EiPsQoWEoGHEH07j2pJDCCQEBycIuIlYsh2BhCDgpACvmLE04eEXiURsapGPEUjwPS4yw3PGhSUZ91PR1QRPWQ0SAgk+RmEEdTCCNzkCCb7E2x5qsmsEL5wIJPgP73YobsIRvNAikOAnvMnhF5RKYyCQ4Bu8t+E7lEojIZDgA7yl4V+USsMjkKA63skIAEqlYXAvO6jLuGMYb2DnRIpvb1wCuxg3V+QmeAMQSFAUd0d1wdX2590JFCm+neQx2BgG92YdgECCiiiMXGNmEmnkGkqlfriGBLUwzu6+q+3PI/J30shlRiZxqHcjkKAQ3p+eiBTf3v3f/yGT3MdMhx4M2UEVpJEnzJG6nutJcAfDd90IJHiP2XRe6bluRCZ5hUwyEEjwGLPpPHR9jI5RO6+QSUIgwVsURoCJGeEEEjxDGgE9Qn5JiUCCB7hoBAwQ2kwikOA2LhoBNwpnJhFIcBWFETCkEF5SIpDgHtIIGEnYLikRSHAJaQSMJzyZRCDBDaQRMImQZBKBBMeRRsDkwpBJBBKcRRoBdgl8JhFIcApfNgJsF+xMIpDgCL5sBDgkwJlEIMF+FEaAo4KaSQQSbEYaAS4IZCYRSLATaQS4JniZRCABgF8FLJMIJNiG8ghwX5AyiUCCPUgjwCuBySQCCTYgjQBvBSOTCCRMijQCVBCATCKQMBHSCFCH3zOJQML4SCNANb7OJAIJYyKNADX5N5MIJIyDNAJgOwIJIyONAMX5tEgikDAa0gjwBT9mEoEEAMHku0wikDACyiPAX/yVSQQShkUaAX7ko0wikDAU0giA0wgk3Iw0AnzNL0USgQQAweeLTCKQcAPKIyAY1M8kAgmDkEYAXEMgoS/SCAgYxYskAgkAQkTlTCKQYI3yCAgqZTOJQIIF0giA+wgk9CKNgMBTs0gikAAgjBTMJGcDqdFodDodR3cBe1EeAfDKJ85tOp/Pa5rWbDaLxWImkzGX37t3L5lMikgymVxbW3OuARgVaQSEilEkqfOudyqQdF3XNO3FixetVqtUKpmB1Gq1ksnkmzdvHNovAGB4SmWSU4HUaDRSqZSIaJpWr9fN5a1WKxaLlUqlaDS6uroai8UcagBGpc5BCSCcHByy0zTNeJFOp82FnU5ndnb2/v37zWazUCj0K5UikYj5+u7du8410tcSiYRdm6pv/ZJ+ekue2rZBpdjYUQFGLw0piB11qy7x9NNbk2zihx9+MF+PfdJ2MJBarZbxortCyuVyuVxORDKZjK7r/X726urKuYYFQyKROD8/t2trH7+L27g1pdjbUUFFLw0pwB1VF9vGSMbObKdm2aVSqQ8fPshvF43M5bu7u7VazaGdYjwM1gFQgVOBlMvl6vV6uVwuFAorKysiUqvVEolEOp0ulUrlcjmfzy8vLzu0dwDASFT4WpKDQ3YHBwe1Wu3Ro0fGxaRMJmOUusfHx93L4S3KIwAGz2fcORhIItL99aNhlsNlpBEAdXDrIADAr7wduCOQwovyCIBSCCQAwO88LJIIpJCiPAKgGgIJAPBfvCqSCKQwojwCMJgnmUQghQ5pBEBNBBIAwIL7RRKBFC6URwCURSABAKy5XCQRSCFCeQRAZQQSAKAvN4skAiksKI8AKI5AAgAM4lqRRCABAJRAIIUC43UAJuFOkUQgAQCUQCAFH+URAF8gkAAAN3Nh1I5ACjjKIwB+QSABAIbidJFEIAUZ5REAHyGQAABKIJAAAMNydNSOQAosxusA+AuBBAAYgXNFEoEUTJRHAHyHQAIAKIFAAgCMxqFROwIpgBivA+BHBBIAQAkEEgBgZE6M2hFIQcN4HQCfIpAAAEogkAAA47B91I5AChTG6wD4F4EEAFACgQQAGJO9o3YEEgBACQRScHABCYCvEUgAgPHZOGpHIAEAlEAgBQTjdQD8jkACACiBQAIATMSuy0gEEgBACQRSEHABCUAAEEgAACVYB1Kj0SiXyyJSKBQSiUSpVHK3VQAAP7HlMpJ1ID179uz+/fuNRqPVap2fnzebzUajMeGeAAAYwDqQms1mJpNpNpvpdFpE0un05eWluw0DAISLdSAtLS2Vy+W9vb1cLqfr+uHhYTKZdLllGBIzGgAEg3UgFYvFqamplZWVVCp1eXn54sWLWCzmcssAAKHySc+/8/l8sVjUdb3ZbIrI3t6esTyTybjdNACAfxjzGiYZsOkNpNXVVU3Tcrnc/fv3J2sbAAAj6A0koxJigA4A4DLra0iFQqFWqxmvG43GwsKCi03CCJjRACAweiskQ7FYLJVKJycnIlKv158/f+5uqwAAoWMdSJqmFYvFx48fi8jLly9TqZS7rQIAhI51IO3u7uq6/vr161gsZpRKa2trLrcMAOAvv91A6NZ4P259DWl2dvbg4CCVSmma9ubNm//85z8TtBAAgJtZB1IulzNe6LpeKpX+8Y9/uNgkAEAYWQ/Z6bpeq9UODw+j0eijR4+q1arLzQIAhE1vIJXL5cPDw1gslsvlisXihw8fVlZWPGkZbsScbwBB0jtk980336RSqefPn6+srGia5kmbAAAh1FshVavVb7/9tlAoxGKxZDIZjUY9aRYAIGx6KyRN09bW1o6Pj58/fx6NRiuVysLCwu7urieNAwCEh/UsOxFJpVJra2vv3r1bWVn58OGDm20CAPiUnXf7vi6Xy5mzwAEAcEjfCgkAADcRSAAAJVgHUqfTcbkdAICQsw6k7e1tY3IdyaSs+tYvfCsWQJBYT2p48eJFo9HQdX1hYSGVSmUymeXlZZdbBgAIlRumfReLxVar9ezZswcPHpTLZTdbBgAIFetAqtVqpVLp3r17h4eHy8vL33///cHBwU8//bS3t+dy+wAAIWE9ZPe3v/3tyy+/LBaLsVjMXLi8vMzd7QAADrEOpD/96U+apnWnkYhkMhlXmgQACCPrQMrlcoVCIZlMmkvevHnjUosAAKFkHUjJZPLly5fO7bXRaFyvwAAAYWYdSKlUyrld5vN5TdOazWaxWGQYcDyR4tt//x+JFN9ebX/udVsAwB59b67aaDQuLy+N17qu53I5W8JD13VN0168eNFqtUqlEoE0BiOHPn4nV9ufk0kAAsM6kMrlcr1e73Q65qhasVi0ZX+NRsMovzRNq9frtmwzVHoSiEwCEBjWgVSpVN69e2d862h5eTmfz9u4S3PueDqd7rdOJBIxX9+9e9fGvfve//17IpEQkfqWGC/MJeiH/hkGvTQkOsrSDz/8YL4e+6Q96HlImqadnJyISDKZbDabdg2vtVot48WACunq6sqWfQVPpPj2/PxcRD5+FzdemEtgKZFI0D83opeGREcNY+zMtr5Tw9LS0sLCQjKZPDw8LJfLlUrFrq/EplIp4/mzrVare1o5hmSM0Zn/ZLwOQGBYB5JxF7tYLGZM/n758qVdgZTL5er1erlcLhQKKysrtmwzbMxMIo0ABIn1kJ3x1IlarSYi9+/ft3eXBwcHtVrt0aNH3IhobOYsO68bAgC2sQ6kZ8+eicjs7Ky5xN752cz2BgD0sA6kw8NDLtwBANzUd1JDo9FwuSkAgDDrO+378ePH3FwVAOAa60BaXl7O5XIuNwUAEGa9gZTP54vFoq7rzWazeznTEAAAjuoNpNXVVU3Tcrmc7bO9AQAYoHdSQyaTsXxMEXMcVJN+euvjd3GvWwEAtrG+hlSv13d2doxJDcat5zRNu7y8PDg4cLV1AIDQsA6kk5OT169fm4/pKxQKz58/39nZMR6M5GLzAABhYf09pGaz2X1fn2g02mw2p6amfvrpJ7caBgAIl75fjC0UCnt7e7qul0olI590XZ+fn3e5fQCAkLAesltbW6vVasbDkG7fvl0sFi8vL2285zcAIJA+fhcXuTXez/a9U0M0GjVnfm9vb+dyOb6KBABwjnUglcvler3e6XTMKeDFYtHFVgEAQsc6kCqVyrt37/b29kRkeXk5n8+72igAQPhYT2owaJpmPG48mUz23EkIKpj54oLvxgIIjL6z7BYWFpLJ5OHhYblcrlQqTGcAADjKOpDW1taKxWIsFnv58qWIML8OAOC0vrPsjDl1qVTKvF8DAADOsQ6kTqdzeHhoXEAyrK2tudUkAIAvffwuPvPFhTxNjPfj1oFUKBRisdjs7OwEDQMAYATWgdRsNt+9e+dyUwAAYdZ3lt3u7q7LTcEYmPkNIDD6TmrY3t7e3t42/3l+fu5KewAAIWUdSN988833339v+ehYAACcYD1kl06nXW4HAMDXfp1iNwHrCknTtGw2azzC3PDmzZtJdgMAwGDWgbS8vMyjyv3CmNcw4QcTAPBcbyDVajVP2gEACLneQHr8+LHlesyyAwA4qjeQCB4AwKhsuXAw6HlIAAC4hkAKAu7XACAACCQAgBIIJADAROz65gmBBABQAoEUEFxGAuB3BBIAQAkEEgBgfDbeuoxAAgAogUAKDi4jAfA1AgkAMCZ7HzVAIAEAlEAgBQqjdgD8i0ACAIzD9keDEkgAACUQSEHDqB0AnyKQAABKIJAAACOz/QKSEEiBxKgdAD8ikAAAo3GiPBICCQCgCAIpmBi1A+A7BBIAYAQOjdcJgRRgFEkA/IVAAgAMy7nySAgkAIAiCKQgY9QOgI8QSACAoTg6XicEUuBRJAHwCwIJAHAzp8sjIZDCgCIJgC8QSAAAJRBIoUCRBGASLozXCYEEAFAEgQQAGMSd8kgIpPBg1A6A4ggkAEBfrpVHQiCFCkUSAJURSAAAa26WR0IghQ1FEgBlEUgAAAsul0dCIIUQRRIANRFIYUQmARjM/fJICCQAQA9P0kgIpNCiSAKgGgIJAPA7r8ojIZDCjCIJgFIIJADArzwsj8TpQGo0Gp1Ox9FdYBIUSQDU8Ylzm87n85qmNZvNYrGYyWTM5ffu3UsmkyKSTCbX1tacawCGYWSSh5+JACjC81OBU4Gk67qmaS9evGi1WqVSyQykVquVTCbfvHnj0H4BAGPwPI3EuSG7RqORSqVERNO0er1uLm+1WrFYrFQqlctlRvMUwcAdABU4OGSnaZrxIp1Omws7nc7s7Oz9+/ebzWahUOhXKkUiEfP13bt3nWukryUSCbs2Vd+yc2uqCfCvZiN6aUiB7Kj61i/pp7fk6fi/2g8//GC+HvukbXMg1Wq1k5OT27dvi0ir1TIWdldIuVwul8uJSCaT0XW933aurq7sbVjwJBKJ8/NzGzdYF+8LdifY3lGBRC8NKZAdZQzWnX9h2wbHzmybAymTyRiXi3RdbzQa8ttFI3OF3d3dZDLZPccBimB2AwBvOTVkl8vldnd3y+VyvV5fWVkRkVqt9vjx43/+85+FQmF+fr7ZbC4vLzu0dwDAMJT6GOrgNaSDg4Narfbo0SPjYlImkzFK3ePj4+7lUAdFEhAqqr3fHQwkEek3NMeQnbLIJABe4dZBABBGCn70JJDQi68lAYGnYBoJgQRLZBIA9xFIsEYmAUGlZnkkBBIAhIqyaSQEEgagSALgJgIJg5BJQJCoXB4JgYQbkUlAMCieRkIgAUAYqJ9GQiBhGBRJAFxAIGEoZBLgX74oj4RAwvDIJMCP/JJGQiBhJGQS4C8+SiMhkAAgqPyVRkIgYVQUSYAv+C6NhEDCGMgkQHF+TCMhkDAeMgmA7QgkjIlMAtTk0/JICCRMgkwCVOPfNBICCRMikwB1+DqNhEDC5MgkQAV+TyMhkGALMgnwVgDSSAgkAPC7YKSREEiwC0US4InApJEQSLARmQS4LEhpJAQS7EUmAa4JWBoJgQTbkUmAC4KXRkIgwQlkEuCoQKaREEhwCJkEOCSoaSQEEpxjZBKxBNgowGkkBBIcNfPFBaUSYJdgp5EQSHABmQRMLvBpJAQS3EEmAZMIQxoJgQTXkEnAeEKSRkIgwU1kEjASY1pQSNJICCS4jEwChmREUXjSSAgkuI9MAm4UqsLIRCDBA3xFCRggnGkkBBK8wleUAEuhTSMhkOAtMgkwhW0Kw3UEEjzG8J2HIsW3Ny6BO0I4heE6AgneY/jOK1fbn5NAKgh5YWQikKAKMskT3ZkUKb692v7c2/aEEGlk+sTrBgC/MzKJN6fLrrY/j8jfSSP3GZ/AOOBNBBLUYtZJvEsRbHz2uo5AgnKMdylvV9dEim/v/u//nJ+fUyS5g49c/XANCYpi9p07ukOIOQ4uYDbdAFRIUBelkgt6SiIqJOdQGN2IQILquKqEAOBz1TAIJPgApRL8i49TwyOQ4BuUSvAdPkWNhECCn1AqwS/48DQGAgn+Q6kExfGZaTwEEnzJLJWEWIJKOCYnQSDBxxjBgzqIoskRSPA9RvDgOT4V2YJAQhAwggevcNTZiEBCcDCCBzcRRbYjkBA0jODBaRxgDiGQEEDdI3git7xtDIKkvvULJbhzCCQElnHWqEucMwgmZ3y+ST+9dX5+7nVbAovHTyDg0k9v8SQLTIhnRriDCgmhwDQ8jMH8EMMx4w4CCSFCLGFIHCSeIJAQOsQSBuDA8BCBhJD675l4nIBAFHmPQEKomWcfTkZhxl9fEQQSIMI4XihRH6uGQAJ+xzheSPDJQ00EEtCLcbyg4qOG4ggkoC8KpsDgs4UvEEjADXoKJuG85h/8yfyFQAKGRTL5BX8gnyKQgJGRTGriz+F3BBIwPpLJc933zKXz/Y5AAmxwPZmE86OT+AQQSAQSYKfu8yMnTXsR9oFHIAFOoWyaHF0XKgQS4DjLsun6f8FACIVW5Orqyv29djqdy8tLTdMs/zeRSPCQ4BtFIt787XxH8Y5SJJ+87SVFOmEYih9Oihj7HO5NhVSpVH7++ee1tTVP9g6oo+fk66NT8yRC8mtiVB4E0sLCQrPZ/Oqrr9zfNaC4wflkuY7iAvArwDUeBNLBwUG5XHZ/v4DvWJ64LU/xA9Z3gYJNgh8xqQHwmQGn+AHBMNi//3Vn7J8VUgc2cekCXa1WOzk5uX379vLysogYFVK/a0iJRMKFJgEAHKL0pIZMJpPJZIZcmSl2ABBCf/C6AQAAiHj1PSQAAHpQIQEAlKBQIHU6nVar5XUr/IQeG6DRaHQ6Ha9b4Rt0Vz/0zBjGPjUpNO3b8vYN9+7dSyaTIpJMJrmzQw9ueNFPPp/XNK3ZbBaLxe7ZNBxOlvp1Fwb0DMfSAGOfmlQJJMvbN7RarWQy+ebNG48apTRueNGPruuapr148aLVapVKJfM8wuFkqV93YUDPcCwNMMmpSZUhu4ODg+u/QKvVisVipVKpXC5TNfew7DGISKPRSKVSIqJpWr1eN5dzOFnq110Y0DMcSwNMcmpSJZAsdTqd2dnZXC43NTVVKBS8bg58w7yRfDqdNhdyOPVj2V2Q/j3DseQQz4bseu7dYCmXy+VyORHJZDK6rrvYOkUN02mhZXaOiJgXVLs/2HI49WPZXZD+PcOx5BDPAmmYezfs7u4mk0kGtU0j3fAibMzO0XW90WjIbwP95gocTpZSqZRld2FAz3AsOUSVSQ3darXa48ePz8/P0+l0oVCYn59vNpvUBBhSLpfb3d0tl8v1en1lZYXDabCe7vK6OQq53jMcS07zwZ0aarWapmn9Hi8LWOp32HA4WaJb+hnQM3Sa7XwQSACAMFB6lh0AIDwIJACAEggkAIASCCQAgBIIJACAEggkAIASCCRgNIlEIp/P5/P5Bw8e3Lt3b29vb8gfrNVq+Xz++vJWq5VIJMrl8qg/CASMindqABRnPneg0Wj89a9/nZ+fj8ViY2/t22+/nZ+fPzw87Pf8mGQyWSwWx94+4BdUSMD4jMcTXF5eisje3p5RM5m1TrlcvnfvXvcSS4eHhysrK7FYzLhNp67ru7u7ItJqtYw7SbdaLV3XO51OqVRKJBILCwvGPdaAgKFCAkZWq9VEpNPp1Gq1paUlTdMajcbe3t7BwYGIFAqFvb094yE67969E5EBEdJoNGKxWCqVyuVytVrNuI30wsJCOp3e3t5eXV0VkcvLy2azWalUotHo+fl5rVbTdd3IQiBIqJCAkZ2cnJycnDQajVarFY1GRUTX9WQy2Ww2m81mMpn88OFDJpN5/fp1rVbb29vrdDpGFXWd8YO1Wi0ajR4eHhoLnz9//vjx4577Sc/Oztbr9b29vWg0yjOzEUhUSMDIzDzodDp/+ctfjH+2Wq2TkxNj+e3btxuNhnFD6Nu3bw+4wlSpVJaWlowfjEajuq7ncjnL9XO5nKZpuq4b5deULD/yAAAA1klEQVTLly/t/8UATxFIwPjM5Lh9+/ZPP/1kJNPu7q4RLfPz88aSfjPxdF3XNM2Mt6mpKWPUrlQqvXz5cnt723yKtrERszZKJBJO/2qA+wgkYCLRaLRWqy0vL+u6ns/nNU1rNpuvX7+OxWLGlIRmsykiRunT87NG/Jj/nJ+fz2azf/zjH43BuufPnz979sy4LiUimqaVSqVGo1Gv15eWltz6/QD38PgJwDaNRuPy8tK88NNqtcyHjTabzckfMNrpdJrNZjQaZUYDAolAAgAogVl2AAAl/H8UkIZx+cNIOwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pzmap(Gd_c2d);hold on;\n",
    "plot(cos(2*pi*(0:0.01:1)),sin(2*pi*(0:0.01:1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 - Réponse harmonique \n",
    "---\n",
    "Tracez la réponse harmonique du filtre continu et superposez-la avec celle de votre filtre numérique. \n",
    "\n",
    "Préférez l'utilisation de vecteurs $f$, $p$, $z$ et des outils **plots** (comme dans les exemples du notebook du [corrige TD FILT prepa bilineaire](../FILT_prepa_bilineaire_octave_corr.ipynb))\n",
    "\n",
    "\n",
    "Faites un affichage du :\n",
    " - gain de la fonction de transfert en continu en noir (plot avec 'k')\n",
    " - de celui de celle en continu en rouge (avec 'r')  \n",
    " \n",
    "avec en abscisse les fréquences linéaires (plot) et en ordonnée le gain en échelle linéaire (abs).\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {
    "vscode": {
     "languageId": "json"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAJMmlDQ1BkZWZhdWx0X3JnYi5pY2MAAEiJlZVnUJNZF8fv8zzphUASQodQQ5EqJYCUEFoo0quoQOidUEVsiLgCK4qINEWQRQEXXJUia0UUC4uCAhZ0gywCyrpxFVFBWXDfGZ33HT+8/5l7z2/+c+bec8/5cAEgiINlwct7YlK6wNvJjhkYFMwE3yiMn5bC8fR0A9/VuxEArcR7ut/P+a4IEZFp/OW4uLxy+SmCdACg7GXWzEpPWeGjy0wPj//CZ1dYsFzgMt9Y4eh/eexLzr8s+pLj681dfhUKABwp+hsO/4b/c++KVDiC9NioyGymT3JUelaYIJKZttIJHpfL9BQkR8UmRH5T8P+V/B2lR2anr0RucsomQWx0TDrzfw41MjA0BF9n8cbrS48hRv9/z2dFX73kegDYcwAg+7564ZUAdO4CQPrRV09tua+UfAA67vAzBJn/eqiVDQ0IgALoQAYoAlWgCXSBETADlsAWOAAX4AF8QRDYAPggBiQCAcgCuWAHKABFYB84CKpALWgATaAVnAad4Dy4Aq6D2+AuGAaPgRBMgpdABN6BBQiCsBAZokEykBKkDulARhAbsoYcIDfIGwqCQqFoKAnKgHKhnVARVApVQXVQE/QLdA66At2EBqGH0Dg0A/0NfYQRmATTYQVYA9aH2TAHdoV94fVwNJwK58D58F64Aq6HT8Id8BX4NjwMC+GX8BwCECLCQJQRXYSNcBEPJBiJQgTIVqQQKUfqkVakG+lD7iFCZBb5gMKgaCgmShdliXJG+aH4qFTUVlQxqgp1AtWB6kXdQ42jRKjPaDJaHq2DtkDz0IHoaHQWugBdjm5Et6OvoYfRk+h3GAyGgWFhzDDOmCBMHGYzphhzGNOGuYwZxExg5rBYrAxWB2uF9cCGYdOxBdhK7EnsJewQdhL7HkfEKeGMcI64YFwSLg9XjmvGXcQN4aZwC3hxvDreAu+Bj8BvwpfgG/Dd+Dv4SfwCQYLAIlgRfAlxhB2ECkIr4RphjPCGSCSqEM2JXsRY4nZiBfEU8QZxnPiBRCVpk7ikEFIGaS/pOOky6SHpDZlM1iDbkoPJ6eS95CbyVfJT8nsxmpieGE8sQmybWLVYh9iQ2CsKnqJO4VA2UHIo5ZQzlDuUWXG8uIY4VzxMfKt4tfg58VHxOQmahKGEh0SiRLFEs8RNiWkqlqpBdaBGUPOpx6hXqRM0hKZK49L4tJ20Bto12iQdQ2fRefQ4ehH9Z/oAXSRJlTSW9JfMlqyWvCApZCAMDQaPkcAoYZxmjDA+SilIcaQipfZItUoNSc1Ly0nbSkdKF0q3SQ9Lf5RhyjjIxMvsl+mUeSKLktWW9ZLNkj0ie012Vo4uZynHlyuUOy33SB6W15b3lt8sf0y+X35OQVHBSSFFoVLhqsKsIkPRVjFOsUzxouKMEk3JWilWqUzpktILpiSTw0xgVjB7mSJleWVn5QzlOuUB5QUVloqfSp5Km8oTVYIqWzVKtUy1R1WkpqTmrpar1qL2SB2vzlaPUT+k3qc+r8HSCNDYrdGpMc2SZvFYOawW1pgmWdNGM1WzXvO+FkaLrRWvdVjrrjasbaIdo12tfUcH1jHVidU5rDO4Cr3KfFXSqvpVo7okXY5upm6L7rgeQ89NL0+vU++Vvpp+sP5+/T79zwYmBgkGDQaPDamGLoZ5ht2GfxtpG/GNqo3uryavdly9bXXX6tfGOsaRxkeMH5jQTNxNdpv0mHwyNTMVmLaazpipmYWa1ZiNsulsT3Yx+4Y52tzOfJv5efMPFqYW6RanLf6y1LWMt2y2nF7DWhO5pmHNhJWKVZhVnZXQmmkdan3UWmijbBNmU2/zzFbVNsK20XaKo8WJ45zkvLIzsBPYtdvNcy24W7iX7RF7J/tC+wEHqoOfQ5XDU0cVx2jHFkeRk4nTZqfLzmhnV+f9zqM8BR6f18QTuZi5bHHpdSW5+rhWuT5z03YTuHW7w+4u7gfcx9aqr01a2+kBPHgeBzyeeLI8Uz1/9cJ4eXpVez33NvTO9e7zofls9Gn2eedr51vi+9hP0y/Dr8ef4h/i3+Q/H2AfUBogDNQP3BJ4O0g2KDaoKxgb7B/cGDy3zmHdwXWTISYhBSEj61nrs9ff3CC7IWHDhY2UjWEbz4SiQwNCm0MXwzzC6sPmwnnhNeEiPpd/iP8ywjaiLGIm0iqyNHIqyiqqNGo62ir6QPRMjE1MecxsLDe2KvZ1nHNcbdx8vEf88filhICEtkRcYmjiuSRqUnxSb7JicnbyYIpOSkGKMNUi9WCqSOAqaEyD0tandaXTlz/F/gzNjF0Z45nWmdWZ77P8s85kS2QnZfdv0t60Z9NUjmPOT5tRm/mbe3KVc3fkjm/hbKnbCm0N39qzTXVb/rbJ7U7bT+wg7Ijf8VueQV5p3tudATu78xXyt+dP7HLa1VIgViAoGN1tubv2B9QPsT8M7Fm9p3LP58KIwltFBkXlRYvF/OJbPxr+WPHj0t6ovQMlpiVH9mH2Je0b2W+z/0SpRGlO6cQB9wMdZcyywrK3BzcevFluXF57iHAo45Cwwq2iq1Ktcl/lYlVM1XC1XXVbjXzNnpr5wxGHh47YHmmtVagtqv14NPbogzqnuo56jfryY5hjmceeN/g39P3E/qmpUbaxqPHT8aTjwhPeJ3qbzJqamuWbS1rgloyWmZMhJ+/+bP9zV6tua10bo63oFDiVcerFL6G/jJx2Pd1zhn2m9az62Zp2WnthB9SxqUPUGdMp7ArqGjzncq6n27K7/Ve9X4+fVz5ffUHyQslFwsX8i0uXci7NXU65PHsl+spEz8aex1cDr97v9eoduOZ67cZ1x+tX+zh9l25Y3Th/0+LmuVvsW523TW939Jv0t/9m8lv7gOlAxx2zO113ze92D64ZvDhkM3Tlnv296/d5928Prx0eHPEbeTAaMip8EPFg+mHCw9ePMh8tPN4+hh4rfCL+pPyp/NP637V+bxOaCi+M24/3P/N59niCP/Hyj7Q/Fifzn5Ofl08pTTVNG02fn3Gcufti3YvJlykvF2YL/pT4s+aV5quzf9n+1S8KFE2+Frxe+rv4jcyb42+N3/bMec49fZf4bmG+8L3M+xMf2B/6PgZ8nFrIWsQuVnzS+tT92fXz2FLi0tI/QiyQvpTNDAsAAAAJcEhZcwAACxMAAAsTAQCanBgAAAAddEVYdFNvZnR3YXJlAEdQTCBHaG9zdHNjcmlwdCA5LjI2WJButwAAFFFJREFUeJzt3TFTG+v58OGH9z01nA/gpUSF1LpALp0ClVFhThdp5kBLg7qQyQxdRKMWCnsmjeWCSRVRnJSWipNSKnA6LV9AfAH+xU40GLCD8UrchusqMmLRavc8bPTzrla7K1dXVwkAHtv/e+wVAICUBAmAIAQJgBAECYAQlhGk2WyW5/kSFgTAj2sZQfrw4cP79++XsCAAflwLD1Kz2Tw6Olr0UgD40S08SKenp7/++uuilwLAj85JDQCE8NNjr0CqVCqPvQoAlOn8/PwBcz1+kNJDV50HW1lZccmo5TPsj8KwL9+DdzMcsgMghGXsIXU6nSUsBYAfmj0kAEIQpOdoY2PjsVcBlsTW/gMRJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQhAkAEIQJABCECQAQigzSOPxeDab3X86AMyVFqRWq9Xv99vt9nA4nE+czWbNZrPf7xf/W9ayAHh6ygnSYDDIsuzw8LDX6x0fH8+nn52dbW5uHh4enp6eDgaDUpYFwJP0UymvMh6Pa7VaSinLstFoNJ9erVb7/f5wOJxMJtVqtZRlAfAklROklFKWZcWDzc3N6xNXV1cHg8FkMtne3v7SvCsrK/PHGxsbZa0SX1GpVB57FZ4jw/4oDPsSfPr0af74wW/jpQUpz/PiwfU9pOPj40ajUaTo9evXX2rS1dVVWavBfVQqlfPz88dei2dnZWXFsC+frX35HvwvgHI+Q6rVatPpNKWU57lDcwA8QDl7SI1G4+TkpNvtjkajnZ2dlNJwOGy327/99lur1ZpOp5PJZGtrq5RlAfAkrZR4uGw4HGZZNv8w6X9OL9ihXj5j/ihWVsr8vxv3ZGtfvgePeWmfIaWU6vX6N00HgDmXDgIgBEECIARBAiAEQQIgBEECIARBAiAEQQIgBEECIARBAiAEQQIgBEECIARBAiAEQQIgBEECIARBAiAEQQIgBEECIARBAiAEQQIgBEECIARBAiAEQQIgBEECIARBAiAEQQIgBEECIARBAiAEQQIgBEECIARBAiAEQQIgBEECIARBAiAEQQIgBEECIARBAiAEQQIgBEECIARBAiAEQQIgBEECIARBAiAEQQIgBEECIARBAiAEQYIluXrsFYDgBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBDKDNJ4PJ7NZren53me53mJCwLg6fmprBdqtVpZlk0mk/39/Xq9Pp9+cHCQUsrzvNFobG9vl7U4AJ6YcvaQBoNBlmWHh4e9Xu/4+Hg+fTgcppSK6cVjALhTOXtI4/G4VqullLIsG41G8+kfP35cX1/v9/sppV6vV8qyAHiSSjtkl2VZ8WBzc/P69Pfv3+/u7k6n01ar9e7duzvnXVlZmT/e2Ngoa5X4ikql8tir8OycG/ZHYtiX4NOnT/PHD34bLy1I89MWru8hpZS2traKj46azeaX5r26uiprNbiPSqVyfn7+2Gvx/KysGPbls7Uv34P/BVDOZ0i1Wm06naaU8jyvVqvz6evr6/PHd56ABwCFcvaQGo3GyclJt9sdjUY7OzsppeFw2G63z8/Pm81mt9udTCa7u7ulLAuAJ2mlxMNlw+Ewy7L5h0n/c3rBDvXyGfPHsbKSHJ1eOlv78j14zEv7DCmldP3rR/eZDgBzLh0EQAiCBEAIggTLcP3LdsCdBAmAEAQJgBAECYAQBAmWxydJ8BWCBEAIN78Ym+f5+/fvU0q7u7tra2uPsUoAPEc395BarVbxoN1uuxwqAEvz2R7SeDzOsqzT6aSU+v1+s9nc2tpqNBpHR0dfupURAJTisyBlWTa/rdH29nZxI4ksy1yoG4BF++yQ3dra2tbW1suXL4ubjtdqtVqttra25uqo8J2uUqpsbLjWN3zFzc+QOp3O77//rkAALNndp31/6d5FALAgvocEQAiCBEAId98xdjabTSaT61N8qgTAQt0dpL/85S8ppRcvXsynCBIAC3V3kM7Ozs7Pz5e8KgA8Z3d/hvTmzZvxeLzkVQHgObt7Dyml1G63iys1FFw6CICFujtI29vbjUZjyasCwHN2d5BqtdqS1wOAZ+6Lh+z6/f50Ol1fX19dXa3X6+6NBMBC3X1SQ7fbHQwGKaXpdHpxcdFut5e7VgA8O3cH6cOHD71e79WrVymlnZ2darU6HA6Xu2IAPC93B2l1dfX6j3me35gCAOW6+zOk3d3d4rTvPM/39vaS0xwAWLAvnvZdrVYHg0G1Wl1fX9/e3l7yagHw3HzxLLvidrHLXBUAnjO3nwAgBEECIARBAiAEQQIgBEECIAS3MAcgBLcwByAEtzAHIAS3MAcgBLcwByAEtzAHIISbQWq1Wvv7+4PBwFl2ACzTzSDt7u5mWdZoNIq78wHActw8qaFer6+trdVqtet35BsMBu4YC8BC3f0ZUrfbHY1Gs9lsbW2tmLK/v7/EtQLg2bk7SB8+fPj999/7/X5KaXt7u9VqLXWlAHh+vnYtuyzLptNpSqlard44xwEAyvXFL8Y2m81qtXp2dtbtdj98+JBl2ZLXDIBn5e4gdTqd/f39tbW1Xq+XUur1eoIEwEJ98UoNxReParVarVZb4voA8EzdEaR+v7+6uppl2dHR0XyiSwcBsFA3D9l1u91+v59l2eXlZUqpOHDnMkIALNrNIJ2dnb19+3Z+mK5Wq/V6veL8bwBYnM+CNB6PsywrvgxbXEComL66uupKDQAs1GdByrIsz/P54+3t7eLx5eXl9SsJAUDpPgvS2tpalmUnJyfXJ/b7/dls5lw7ABbq5ll2vV6v3W4PBoPNzc2UUnFFu+LbSACwODeDtLa2dnp6OhgMiluY7+zsOMUOgCW4+4uxjUZDhwBYpq9dXBUAlkaQAAhBkAAIQZAACEGQAAihzCCNx+PZbHbnr/I8/9KvACCVGKRWq9Xv99vt9u2r3s1ms2azWVwiDwDuVE6QBoNBlmWHh4e9Xu/4+PjGb4+OjlZXV+0hAfAV5QRpPB4XF7vLsmw0Gl3/1cnJyfr6+vwi4gBwpy/ewvxbZVlWPCguglcYj8fj8bjX6338+PEr866srMwfb2xslLVKfEWlUpk/Pv/06RHX5Gmr/Hd7Pp9Pura1swTnyZgvROXz9+pP195GHvw2XlqQ5vetuL6HdHJy8uLFi263m+f5wcHB7u7uvFvXXV1dlbUa3EelUjk/P//fz+O7nX/+Vrhia186W/uCfGVMr/9795uUE6RarVZcjDXP82q1Op++s7NT3Ap9NBo1Gg03VeLZ+W9+VlZW7PvD15UTpEajcXJy0u12R6PRzs5OSmk4HLbb7fk/TFZXV+v1einLAuBJKu2Q3enp6XA4/OWXX4qDcvV6/fpu8rt378paEABPUmlBSinZBwLgwVw6CIAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBAECYAQBAmAEAQJgBDKDNJ4PJ7NZndOz/O8xAUB8PT8VNYLtVqtLMsmk8n+/n69Xi8mzmazdrtdrVbzPK9Wq51Op6zFAfDElLOHNBgMsiw7PDzs9XrHx8fz6R8+fNjc3Dw8PHz37t3Z2VkpywLgSSpnD2k8HtdqtZRSlmWj0Wg+fWtrq3hw56E8AJgr7ZBdlmXFg83NzRsTh8Ph0dHR7u7ul+ZdWVmZP97Y2ChrlfiKSqXy2KvwTBn55TPmS/Dp06f54we/jZcWpPlpC9f3kFJK3W734uKi1+vNi3Xb1dVVWavBfVQqlfPz88dei+dl/q8uI79ktvble/C/AMoJUq1WG4/HKaXi5IX59H6/X9SolKUA8ISVE6RGo3FyctLtdkej0c7OTkppOBy22+03b97ked5qtYqnvXv3rpTFAfD0lHbI7vT0dDgc/vLLL8WhuXq9bjcZgPsrLUgppfnXjwDgW7l0EAAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAIQgSACEIEgAhCBIAISwjSOPxeDabLWFB3NOnT58eexWeI8P+KAz7D2ThQWq1Wv1+v91uD4fDRS8LgB/XTwt99cFgkGXZ4eFhnucHBwf1en2hiwPgx7XYPaTxeFyr1VJKWZaNRqOFLguAH9rCD9llWVY82NzcXPSyAPhxLfaQXUopz/PiwVf2kCqVyqJXg+s2NjaM+ZJtbGwUD4z8ktnafyCLDVKtVhuPxymlPM+r1eqdzzk/P1/oOgDwQ1i5urpa6AKazebm5uZoNNrZ2Wk0GgtdFgA/roUHKaU0HA6zLJt/mAQAty0jSADwPy3wLLvZbDY/o6Fw45INt6/g8K1P4P6M9qMznkvwgK3a3+Wevv8t/X/O8v//+te/lrvSc3//+9///e9/v3r1qvix1WrleX5yclIcvrvx4wOesKDVfhpevnw5HA7/8Y9//Oc//3n16pXRfnTGs3Q3NvL0oK3a3+X+vvMt/T6zpKvF+OMf/7ixsfG3v/2t+PGf//znn//856urq+l0+qc//enGjw94woJW+2m4MURG+9EZz9LdHskHbNX+Lvf3nW/p95nl6upqUYfsTk9Pf/311/mPNy7ZcPsKDt/6hAWt9tOQ5/na2trBwUG3253NZkb70RnP0t3YyNODtmp/l/v7zrf0+8ySlnn7iRuXbLh9BYdvfQJfMpvNXrx40Wg0fv755729vWS0AzCe5bq9kacHbdX+Lg/2/aN9+wllBmk4HHa73X6/f+dvb1yy4fYVHL71CVx3ffAbjUan06nX6zs7O5eXl8loB2A8y3V7I08P2qr9XR7s+0f79hPKDFK9Xu90Otvb27d/VavVptNp+u8lG278+IAnlLjaT8P1wT85Obl+sw+j/eiMZ+lubOTpQVu1v8uDff9o3zn4C/weUrfbTSl1Op3ixxuXbLh9BYdvfcKCVvsJGI/He3t7W1tbk8mk0Whsb28b7UdnPMt1eyNPD9qq/V3u7zvf0u8zy6LOsrvTx48fp9Ppl358wBP4CqMdjfEsXSlbtb/Lg33/aN/40ZUaAAhheWfZAcBXCBIAIQgSACEIEgAhCBIAIQgSPFnFlQwfey3gvn567BWAZatUKtcvXLa/v19c5PGJ6Xa7Z2dnWZblef7u3Tv3ViA+30Pi2alUKufn54+9Fgv3+vXrf/3rX+nWF+whLIfsII3H436/32q1ivfufr//+vXrly9fFj+mlLrd7suXL5vNZr/f7/f74/F4/qvix9tzjcfjg4ODg4ODSqWyt7dX3CKheJ3Xr1/Pr0F8Y67ZbFbM0mw2bxxtu8/rz+V5vru7u8gxg/I5ZMdzNM/J+vr69vb25eXl0dFRr9er1+tFnE5PT1NKe3t7RTkmk8lvv/2WUmo2m1tbW1mWTSaT4hWm02mWZbfnyrLs7Ozs7du3h4eHe3t7xfPnr/OHP/xha2srz/Mbc11eXq6urp6fnw+Hw8FgMD+WeJ/Xr9fr8//ALMuKy7v1+/3RaPT27dvlDCx8D0HiOZrfhnl1dbV4UK1Wizf0wWBQrVaLflSr1el0enFxsb29vba2llLa2tq68wVvz5VlWXFJ45TSixcvUkofP35sNBrF6xSFuD1XrVYbDAb9fr9arV4/yHaf179tPB4fHx+fnp4WC4XgBInn6PrOxG15nn/8+LF4vL6+PplM7vOGfmOu20+4uLiYh3C+63NjrkajkWVZ0aQsy3q93v1f/05bW1tqxI/CZ0jwmfX19bW1tU6n0+l0fv7559XV1Wq1Oi/B/EjdXHFvsdtz3X7lWq02f53iI6Lbc/X7/TzPO53O6enp2dnZV9bqPv8tWZa5pQI/EHtI8Jnt7e3BYNBqtYoPit6+fVuv15vN5sXFxfwGl8XRs263e3FxUbTh9ly30/XmzZt2u31wcDC/QVlxgO7GXAcHB+PxeDQavXnz5itrdfv1b5tMJsfHx+/evSttdGCRnPYNdxiPx5eXl9eP7A2HwyzL3r9/n1LqdDqz2WwymVSr1esHxG7PddtwOFxdXb3+zacbcxWvfOM59399+HEJEnwD3+mBxXHIDr7B/KwEoHT/B2qIN1DG2hHrAAAAAElFTkSuQmCC",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Df = Fe/100;        % résolution fréquentielle d'affichage\n",
    "f = -Fe: Df: 2*Fe ; % vecteur des fréqeunces [Hz]\n",
    "\n",
    "%%VOTRE CODE ftilde = , p= , z=\n",
    "f_tilde = f';      % fréquences normalisées [1<-> Fe]\n",
    "p = exp(f);        % équivalents harmoniques de p\n",
    "z = exp(f_tilde);  % équivalents harmoniques de z\n",
    "\n",
    "\n",
    "% fonction de transfert continue\n",
    "Gc = @(s,wn,zeta) wn^2./(s.^2 + 2*zeta*wn*s + wn^2)  ;\n",
    "\n",
    "%%VOTRE CODE G_continu=..., Gd_math=...\n",
    "% utilisez la fonction anonyme Gd qui prend l'équivalent de z en argument\n",
    "% pour calculer la fonction de transfert pour chaques fréquences\n",
    "G_continu = Gc(f,1,1);\n",
    "Gd_maths = Gd(f);\n",
    "\n",
    "%%VOTRE affichage ici \n",
    "plot(f,G_continu,'k');hold on;\n",
    "plot(f,Gd_maths,'r');grid on;\n",
    "\n",
    "%%VOS BONS LABELS\n",
    "xlabel(\"Frequences en ?\");\n",
    "ylabel (\"Gain  en ? \");"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vous devez constater que :\n",
    " - le spectre discret est périodique\n",
    " - les fréquences de résonnance ne collent pas !\n",
    " \n",
    "C'est l'effet de la fameuse compression des fréquences (\"frequency wrapping\") induites par la méthode de la transformation bilinéaire. \n",
    "\n",
    "Faites un vrai **tracé de bode** en changeant les fréquences de calcul, car l'échelle log sur l'axe des fréquences ne permet pas d'afficher les fréquences négatives ou nulles : On peut prendre des **fréquences de 100Hz à 10 kHz** par exemple.\n",
    "\n",
    "On utilise `semilogx` à la place de plot pour avoir les abscisses en échelle log.\n",
    "On peut caluler le gain en dB avec une fonction anonyme `dB_de`.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [],
   "source": [
    "Df = Fe/1000;        % résolution fréquentielle plus fine\n",
    "f = 100: Df: 10*Fe ; % fréquences positives car échelle log !\n",
    "\n",
    "dB_de = @(gains_complexes) 20*log10(gains_complexes);\n",
    "\n",
    "%% COPIER/COLLER de la cellule précédente\n",
    "%% PUIS modif de l'affichage pour passer en \n",
    "% échelle log sur x (semilogx)\n",
    "% dB sur ordonnées (utiliser dB_de)\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si on cherche les fréquences de résonnance en continu et en discret avec la fonction max. On voit que l'on n'a pas le $\\omega_n$ voulu...\n",
    "\n",
    "On ne l'a pas en continu, car on a confondu pulsation naturelle $\\omega_n$ et pulsation de raisonnance $\\omega_r$\n",
    "\n",
    "En discret, il va falloir comprendre le phénomène de compression dû à la bilinéaire."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "freq_naturelle =  0.15915\n",
      "freq_raisonnance =  0.0027658\n",
      "freq_r_continu_mesuree =  4900\n",
      "freq_r_discret_mesuree =  4900\n"
     ]
    }
   ],
   "source": [
    "freq_naturelle = wn/2/pi    % on pense que la fréquence wr est wn !!!\n",
    "wr = wn * sqrt(1-2*zeta** 2);\n",
    "freq_raisonnance = wr/2/pi\n",
    "[gain_max, indice]=max(abs(G_continu));\n",
    "freq_r_continu_mesuree = f(indice) % on à la fréquence du max à Df près\n",
    "[gain_max, indice]=max(abs(Gd_maths(1:indice))) ;\n",
    "freq_r_discret_mesuree = f(indice) % en discret ça colle pas\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  2.1 - Compression fréquencielle \n",
    "\n",
    "Regardez la [suite de la vidéo](https://youtu.be/uBZOh7vhjHU?t=587) sur l'effet de compression de la méthode bilinéaire.\n",
    "\n",
    "Dans la transformation $p = \\frac{2}{Te}\\frac{z-1}{z+1}$ utilisée, remplacez p par $i\\omega_c$ d'un côté et z par son expression $e^{T_e\\,i\\,\\omega}$. On obtient un nombre imaginaire des deux côtés, mais avec une pulsation continue $\\omega_c$ et une pulsation pour le discret $\\omega$ qui sont différentes. \n",
    "\n",
    "\n",
    "Vous retrouvez ainsi la relation de compression entre pulsation en coninue $\\omega_c$ et pulsation équivalente par la transormation bilinéaire $\\omega$ :\n",
    "$$ \\omega_c = \\frac{2}{T_e} \\tan\\left(\\frac{T_e}{2}\\omega\\right)$$\n",
    "\n",
    "On a donc égalité entre les expressions  \n",
    "$G_d(z\\equiv e^{i\\omega\\,T_e}) = G(p = \\frac{2}{Te}\\frac{z-1}{z+1}) = G(p = \\frac{2}{Te}\\frac{e^{i\\omega\\,T_e}-1}{e^{i\\omega\\,T_e}+1}) = G\\left(i\\, \\frac{2}{T_e} \\tan\\left(\\frac{T_e}{2}\\omega\\right)\\right) = G(i\\, \\omega_c)$  et donc   \n",
    "$G_d(z\\equiv e^{i\\omega\\,T_e}) = G(p\\equiv i\\omega_c) = G(\\frac{2}{T_e} \\tan\\left(i\\,\\frac{T_e}{2}\\omega\\right))$ \n",
    "\n",
    "Ainsi, si le polynome $G(p)$ a une résonnance pour $\\omega_c=\\omega_n$ en continu, en discret $G\\left(\\frac{2}{Te}\\frac{z-1}{z+1}\\right)$ a une résonnace pour $\\frac{2}{T_e} tan(\\frac{Te}{2}\\omega)=w_n$.\n",
    "\n",
    "Retrouvez donc analytiquement cette fréquence de résonnance en discret et vérifiez que c'est bien cela qui est observé sur la réponse harmonique."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fr =  0.0027658\n",
      "freq_r_discret =  0.0027661\n",
      "freq_r_discret_mesuree =  4900\n"
     ]
    }
   ],
   "source": [
    "wr = wn * sqrt(1-2*zeta** 2);\n",
    "fr = wr/2/pi\n",
    "\n",
    "%% VOTRE CODE wr_discret=\n",
    "% donnant la pulsation wrdiscrete, une fois wr compressée\n",
    "% faut-il tan ou atan ?\n",
    "wr_discret = tan(wr) ;\n",
    "\n",
    "%% on vérifie si c'est pareil\n",
    "freq_r_discret = wr_discret/2/pi\n",
    "freq_r_discret_mesuree"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 - Compensation fréquentielle \n",
    "\n",
    "Pour terminer, effectuez une précompensation de $\\omega_n$ (en anglais \"prewarp\"):\n",
    "- on cherche un $\\omega_n'$ tel que la résonnance obtenue avec $\\omega_c=\\omega_r'$ donne une fois compressée par la transformation  une résonnance en discret $\\frac{2}{T_e} tan(\\frac{Te}{2}\\omega)=w_r'$ lorsque $w = \\omega_r$.\n",
    "\n",
    "> On prend en continu un $\\omega_r'$ plus grand qui, une fois compressé, donnera le $\\omega_r$ recherché.  \n",
    "\n",
    "> On ne fait pas la confusion entre la pulsation naturelle $\\omega_n$ et la pulsation de raisonnance $\\omega_r = \\omega_n \\sqrt{1-2\\zeta^2}$\n",
    "\n",
    "Il faut donc tout reprendre dans cette cellule de calcul :\n",
    " - le calcul du $\\omega_n$ (naturel) en fonction du $\\omega_r$ voulu (raisonnance)\n",
    " - le calcul du $\\omega_n'$ (prewarp) qui compense la compression (il doit donc être plus grand)\n",
    " - le calcul des coefficients $a_0$ à $b_2$ avc ce wn'\n",
    " - definir une fonction de transfert Gd_prewarp avec ces gains\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "metadata": {},
   "outputs": [],
   "source": [
    "clear all; close all; clc;\n",
    "%% Parametres du filtre continu désiré\n",
    "wr = 2000 ;% pulsation de raisonnance voulue [rad/s]\n",
    "zeta = 0.25; % facteur d'amortissement [1]\n",
    "\n",
    "%% VOTRE calcul de wn=  en fonctionde wr\n",
    "wn = wr/sqrt(1-2*zeta^2) ;      % pulsation propre [rad/s]\n",
    "\n",
    "%% COPIER/COLLER de la toute première cellule \n",
    "%% MAIS avec votre calcul du wn'\n",
    "%% wn_wrap = \n",
    "wn_wrap =  atan(wn);\n",
    "\n",
    "%% ET calcul  des a0,a1,...,b2 AVEC wn_wrap\n",
    "%% PUIS redefinir une fonction anonyme Gd_warp avec ces valeurs\n",
    "Gd_warp = @(z) (b0 + b1 .* z.^(-1) + b2 .* z.^(-2)) ./ (a0 + a1*z.^(-1) + a2*z.^(-2) ) ;\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "On refait le tracé de bode et on vérifie les maxs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [],
   "source": [
    "%% COPIER / COLLER \n",
    "% afficher le continu voulu en noir \"k\"\n",
    "% le continu prewarp en bleu \"b\"\n",
    "% le discret en rouge \"r\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les fréquences de résonnance rouge et noire doivent coller !\n",
    "\n",
    "Mesurez les fréquences de raisonnance et vérifier que cela colle avec $\\omega_r$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fr =  318.31\n",
      "freq_r_continu_mesuree =  42\n",
      "freq_r_discret_mesuree =  42\n"
     ]
    }
   ],
   "source": [
    "%% VOTRE CODE freq_r_continu_mesuree = \n",
    "freq_r_continu_mesuree = 42; \n",
    "%% VOTRE CODE freq_r_discret_mesuree = \n",
    "freq_r_discret_mesuree = 42; % en discret\n",
    "\n",
    "%% ça doit coller !\n",
    "fr = wr/2/pi\n",
    "freq_r_continu_mesuree\n",
    "freq_r_discret_mesuree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 - Application : filtrage \"overtone\"\n",
    "---\n",
    "\n",
    "Vous êtes enfin prêtes pour synthétiser un filtre dans le [tp2_overtone](../tp2_overtone/tp2_overtone.ipynb) pour le même effet d'overtone que la chanteuse Ana Maria Hefele.\n",
    "\n",
    "Explication sur \n",
    "[![](https://i.pinimg.com/originals/fa/fe/6d/fafe6d300d6fa568aed43b31fd3bc67d.jpg)](https://www.youtube.com/watch?v=UHTF1-IhuC0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Octave",
   "language": "octave",
   "name": "octave"
  },
  "language_info": {
   "file_extension": ".m",
   "help_links": [
    {
     "text": "GNU Octave",
     "url": "https://www.gnu.org/software/octave/support.html"
    },
    {
     "text": "Octave Kernel",
     "url": "https://github.com/Calysto/octave_kernel"
    },
    {
     "text": "MetaKernel Magics",
     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "octave",
   "version": "4.2.2"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}