Numerical interpolation and Curve Fitting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical Interpolation and Curve Fitting\n",
    "_By Dhruv Jain_\n",
    "\n",
    "### **Objective: Use given data points to determine a value of an intermediate point**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Key libraries: Numpy(for mathematical procedures) and matplotlib(to create plots)\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt \n",
    "import copy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lagrange Interpolating Polynomial Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lagrange Interpolating Polynomial Function:\n",
    "def lagrange_interp(xi,fi,x_interp):\n",
    "    \"\"\"Dhruv Jain, 18 Oct, 2021\n",
    "    Obj: Use Lagrange Form of the Interpolating Polynomial to find f(x_i) for given x_i using given (x,y) data points  \n",
    "    Args:\n",
    "        xi: (n x 1) vector, x-coordinates of the given data set\n",
    "        fi: (n x 1) vector, y-coordinates of the given data set\n",
    "        x_interp: (n_interp x 1) vector, given x-coordinate points for which we interpolate and finds its y-coordinates\n",
    "    Output:\n",
    "        f_interp: (n_iterp x n) vector, interpolated y-ccordinates of x_interp using Lagrange form of Interpolating Polynomial\n",
    "    \"\"\"\n",
    "    n = len(xi)\n",
    "    coeff = np.zeros(n)\n",
    "    coeff = copy.copy(fi) # coeff contains all the f(xi) terms \n",
    "\n",
    "    # Step 1: Create fi/denominator(Ln,i); independent of the points we want to interploate for\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            if i==j:\n",
    "                continue\n",
    "            coeff[i] = coeff[i]/(xi[i]-xi[j])\n",
    "    \n",
    "    # Step 2: Multiply the numerator(Ln,i); dependent on x_interp, thus needs to be called len(x_interp) times\n",
    "    def lag_interp_num_coeff(xval):\n",
    "        \"\"\"Dhruv Jain, 18 Oct, 2021, MSAAE Purdue\n",
    "        Obj: Calculate and multiply the numerators of Ln,i based of x_interp[i] value to get f_interp[i]\n",
    "        Args:\n",
    "            xval: float, value for which the interpolated value is calculated\n",
    "        Output:\n",
    "            sum: float, interpolated y-ccordinates of x_interp using Lagrange form of Interpolating Polynomial\n",
    "        \"\"\"\n",
    "        sum = 0.0;\n",
    "        for i in range(n):\n",
    "            temp_c = coeff[i]\n",
    "            for j in range(n):\n",
    "                if i==j:\n",
    "                    continue\n",
    "                temp_c = temp_c*(xval-xi[j])\n",
    "            sum = sum + temp_c\n",
    "        return sum\n",
    "    \n",
    "    # Initialize and calculate y-coord for each x_interp element\n",
    "    f_interp = np.zeros(len(x_interp))\n",
    "    for i in range(len(x_interp)):\n",
    "        f_interp[i] = lag_interp_num_coeff(x_interp[i])   \n",
    "    \n",
    "    return f_interp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Newton Interpolating Polynomial Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Newton Interpolating Polynomial Function: \n",
    "def newton_interp(xi,fi,x_interp):\n",
    "    \"\"\"Dhruv Jain, 18 Oct, 2021\n",
    "    Obj: Use Newton Form of the Interpolating Polynomial to find f(x_i) for given x_i using given (x,y) data points  \n",
    "    Args:\n",
    "        xi: (n x 1) vector, x-coordinates of the given data set\n",
    "        fi: (n x 1) vector, y-coordinates of the given data set\n",
    "        x_interp: (n_interp x 1) vector, given x-coordinate points for which we interpolate and finds its y-coordinates\n",
    "    Output:\n",
    "        f_interp: (n_iterp x n) vector, interpolated y-ccordinates of x_interp using Lagrange form of Interpolating Polynomial\n",
    "        coeff: (n-1 x 1) vector, coefficients of Newton form of Interpolating Polynomial: f[x0], f[x0,x1], f[x0,x1,x2].......,f[x0,x1,....,xn]\n",
    "    \"\"\"\n",
    "    n = len(xi)\n",
    "    coeff = np.zeros(n)\n",
    "    coeff[0] = copy.copy(fi[0]) # coeff contains all the f[x0], f[x0,x1]....f[x0,x1,.....,xn] terms \n",
    "    temp = copy.copy(fi)\n",
    "    \n",
    "    # Step 1: Compute f[x0], f[x0,x1] ... f[x0,x1,....xn] coeffecient values; independent of the points we want to interploate for\n",
    "    for i in range(1,n):\n",
    "        for j in range(n-i):\n",
    "            temp[j] = (temp[j+1]-temp[j])/(xi[j+i]-xi[j])\n",
    "        coeff[i] = temp[0]\n",
    "    \n",
    "    # Step 2: Multiply the numerator (x-x0)(x-x1)...(x-xn); dependent on x_interp, thus needs to be called len(x_interp) times\n",
    "    def newton_interp_num_coeff(xval):\n",
    "        \"\"\"Dhruv Jain, 18 Oct, 2021, MSAAE Purdue\n",
    "        Obj: Calculate and multiply (x-x0)(x-x1)...(x-xn) based on x_interp[i] value to get f_interp[i]\n",
    "        Args:\n",
    "            xval: float, value for which the interpolated value is calculated\n",
    "        Output:\n",
    "            sum: (n_iterp x n) vector, interpolated y-ccordinates of x_interp using Lagrange form of Interpolating Polynomial\n",
    "        \"\"\"\n",
    "        sum = coeff[0]\n",
    "        for i in range(1,n):\n",
    "            temp = coeff[i]*(xval-xi[i-1])\n",
    "            for j in range(1,i):\n",
    "                temp = temp*(xval-xi[i-j-1])\n",
    "            sum = sum + temp\n",
    "        return sum\n",
    "    \n",
    "    # Initialize and calculate y-coord for each x_interp element\n",
    "    f_interp = np.zeros(len(x_interp))\n",
    "    for i in range(len(x_interp)):\n",
    "        f_interp[i] = newton_interp_num_coeff(x_interp[i])   \n",
    "    \n",
    "    return f_interp, coeff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f95d1e222b0>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABN4UlEQVR4nO3dd1zV5fvH8dcFgogDc4QbNPcCBMSd+G1ou77OSDMzstR2Niw1yzKzZfWtbFlGaqWZabYlM0UZgnuLinsrKipw/f7gwA/xMAQOHzjcz8fjPDp85vs+5Ln4rPsWVcUwDMMwcnKxOoBhGIZROpkCYRiGYdhlCoRhGIZhlykQhmEYhl2mQBiGYRh2mQJhGIZh2GUKhGEUkIhEisjwQq7bSESSRcS1uHNZQUQSReQ6q3MYjmUKhGGpK/miKcoXdEnL2S5V3a2qVVQ1zQH7UhE5YytAe0XkrZIsRCIyQUS+Lqn9GSXHFAij3HCWv95z4aeqVYBrgQHAMIvzGE7AFAij1BCRoSKyTESmishxEdkpIn1s8yYB3YH3bX8pv2+b3lJEfheRYyKyWUT6Z9veDBH5UER+FpEzQKht2ke2dU6LyN8i4pNtnS4iEi0iJ23/7ZJL1mtE5C8ROSoiR0QkQkSq2+bNBBoBP9myjhERX9tf+hVsy9QTkQW23NtE5IFs254gIt+KyFe2jOtFJKggn6GqbgP+Bfyzbe8WEYkXkRMislxE2meb94ztqOO07fP7T7bP7pVsy/UUkSQ7n0Nv4HlggK2tCQXJaZQNpkAYpU0IsBmoBUwBPhMRUdWxwD/AKNupmlEiUhn4HfgGuBoYBPxPRNpk297dwCSgKrDMNi0MeNm2j3ggAkBEagCLgGlATeAtYJGI1LSTU4DXgHpAK6AhMAFAVQcDu4FbbVmn2Fl/FpBkW78v8Grml7PNbcBsoDqwAHg/rw8tK5RISzIK6Tbbzx2Az4EHbW36GFggIhVFpAUwCghW1arAjUBiQfaTSVV/AV4F5tja6ncl6xulmykQRmmzS1U/sZ2r/xKoC3jnsuwtQKKqfqGqqaoaB8wl4ws304+q+q+qpqtqim3aIlVdqqrngbFAZxFpCNwMbFXVmbbtzQI2Abfm3LGqblPV31X1vKoeJqOYXFuQBtr21Q14RlVTVDUe+BQYnG2xZar6s+1zmAnk98UbZztK2ghEAv+zTX8A+FhVV6pqmqp+CZwHOgFpQEWgtYi4qWqiqm4vSBuM8sEUCKO0OZD5RlXP2t5WyWVZHyDEdurkhIicIOPooE62ZfbYWS9rmqomA8fI+Eu+HrArx7K7gPo5NyAiV4vIbNvpmVPA12QckRREPeCYqp7OYz8Hsr0/C3hknp7KRQcyPqcBZByFVbZN9wGezPEZNQTq2U5HPUbGkc8hW3vqFbANRjlgCoRRluTsengP8LeqVs/2qqKqD+WxDmR8QQIgIlWAGsA+28snx7KNgL12tvGabdvtVbUacA8Zp53y2m+mfUANEalagP0UmGb4FlgBjLNN3gNMyvEZedqOjlDVb1S1GxntVuB123pnAM9sm89edC/bdVFyG6WXKRBGWXIQaJLt54VAcxEZLCJutlewiLTKZzs3iUg3EXEn41rESlXdA/xs297dIlJBRAYArW37yakqkAycEJH6wNP5ZM1i29dy4DUR8bBdNL4f27WQYjAZCBeROsAnwAgRCZEMlUXkZhGpKiItRKSXiFQEUoBzZJx2goxrMzeJSA3bdh7LY38HAV8RMd8nTsb8Qo2y5F2gr+0Op2m2UzQ3AAPJ+Kv8ABl/AVfMZzvfAOPJOLUUSMZpKVT1KBnXNZ4EjgJjgFtU9YidbbxExmmdk2Rc2J6XY/5rwAu20zpP2Vl/EOBry/0DMF5Vf88nd4Go6lrgb+BpVY0h4zrE+8BxMi5eD7UtWpGMYnKEjM/uajLuSIKM6x4JZFy0/g2Yk8cuv7P996iIxBVHG4zSQcyAQUZ5IiIzgCRVfcHqLIZR2pkjCMMwDMMuUyAMwzAMu8wpJsMwDMMucwRhGIZh2JXXgzdlTq1atdTX17dQ6545c4bKlSvnv6ATMW12fuWtvWDafKViY2OPqGpte/OcqkD4+voSExNTqHUjIyPp2bNn8QYq5UybnV95ay+YNl8pEcnZe0AWc4rJMAzDsMsUCMMwDMMuUyAMwzAMu0yBMAzDMOwyBcIwDMOwy2EFwtZL5SoRSbANmfiSnWVERKbZhlxcYxv9KnNeb9sQiNtE5FlH5TTKh4iICHx9fYmNjcXX15eIiOLqONUwnJcjb3M9D/RS1WQRcQOWichiVY3KtkwfoJntFQJ8SMYAMK7AB8D1ZAzLGC0iC1R1gwPzGk4qIiKC8PBwzp7NGH9o165dhIeHAxAWFmZlNMMo1RxWIDSjD49k249utlfOfj1uB76yLRslItVFpC4Z3SBvU9UdACIy27asKRBGntJSUzm0dzvH9+3gzKGdpJ09QdKfP/Pif1txPhW8kxN48BY/9hw9x+cfvEHf/95JRQ/P/DdsGOWQQ/tish0JxAJNgQ9U9Zkc8xcCk1V1me3nP4FnyCgQvVV1uG36YCBEVUfZ2Uc4EA7g7e0dOHv27EJlTU5OpkqV3Ea2dE7O0OaLKWc4u2c1lY+vp8657TRJT8RTzhd8fXUlURqwv2JjTni1xq1+AJW8rnZg4pLlDL/jK2XafGVCQ0NjVTXI3jyHPkltG3DdX0SqAz+ISFtVXZdtEbG3Wh7T7e1jOjAdICgoSAv7NKF5+rLsOHHkAJv/+ooq2xbQ/PwG3CSNc+pOontT1lx1Oy5Xt6LS1Y2pXrcJlb1q0bnbtWzdugN39wpMenkCH7//Dj51vGjb2Jv+/wmg8rENtEuJ46pDkXAIdrk0ZF/93jToMZiGzfysbm6RlNXfcVGYNhefEulqQ1VPiEgk0BvIXiCSyDY+MNCAjBG23HOZbpRTmp7OxpW/krLsA9omLydE0kh0aUhM/TC82t1Es8BetHK3P5DcuPEv/f81CNeKbN1zmL1HzzD00Ql0tl2D0PR0dmyI5lD8Yqru/pOQ3Z/iEvEJW12bcqLNENr1GY5HpfLVv49hOKxAiEht4KKtOFQCruP/B0TPtAAYZbvGEAKcVNX9InIYaCYijckYyH0gcLejshqlV1pqKqt/+Ryv1R/TOm0bx6lKXJ3+1O46hCZtO+Hrkv+NeJkXoseOHQuAj48PkyZNuuQCtbi40KRtCE3ahgATOLR3JzsiZ+K9/XuC14zj+JqprG7Ql5a3j+Gq2nUd0lbDKG0ceQRRF/jSdh3CBfhWVReKyAgAVf2IjEHibyJjnNyzwH22eakiMgr4FXAFPlfV9Q7MapQymp5OwpJv8fr3VYLSd7HbpT4r27yI380j6OR55edaw8LCCAsLIzIyksTExHyXv7p+Y64OG4emv8C6FYu4+O//CNnzBWffn8WKRmG0+e/zVKtesxAtM4yyw5F3Ma0BAuxM/yjbewVG5rL+z2QUEKOc2b42igsLnsD/4nqSpC6xwW8S0Ps+Grm6lngWcXGhbddboeut7NoYy9FFL9F5z6ecfGcWUS1GEdT3KSq4uZd4LsMoCeZJaqPUOHP6BFEfjsDn+z54X9zDytYv4P1cAoE3D8fFguKQk0+rQDo8tYDtdy1md8UWdNr8OrtfC2bjyl+tjmYYDmEKhFEqrFu2gNNvBtLp4Cxia96C6+gYQvo/jVsuF56tdE37LrR95k9Wd56GZ3oyrRb3Z+V793Lm9AmroxlGsTIFwrDU+ZSzRH04grZ/DOaCS0U23fQ9IY/MxKumt9XR8iQuLgTceC/VnoojynsQwUd+5MRbHdkQ9YvV0Qyj2JgCYVhmz9YE9k7pTKeDs1hZ8w5qP7mSlh2vtzrWFfGs4kWnhz5i801zAGi5eCArPh5J6sULFiczjKIzBcKwRPzv31D96xu5Kv0o8d0/JmT0l1SqXNXqWIXWKuRGqj+xiuiat9J5/9dsmRLK4X2JVscyjCIxBcIoUelpaaz47En8/32IgxXqkXLfn/j/Z6DVsYpF5arVCXlkJjGBU/C9sBWX6T1Yt2yB1bEMo9BMgTBKTMrZZBLeuo3Oez4l2qs3DZ5cSl2fFlbHKnZBtz7I4YGLSXapSqvfh7ByzmSrIxlGoZgCYZSI44f3k/jWdfgl/0tU86cIenQWHoV44K2s8GkVSM3HlrG2cidCNr7Gyg/uN9cljDLHFAjD4fbu2Ejy/3rR+OI24ju/Q6e7X0QK0EVGWVel2lW0e2IhUd6DCDn8PevfvInTJ49ZHcswCsz5/5UalkrcGEPFr3pTVU+x86Zv6NB7qNWRSpRrhQp0eugjVrYZR5tzsRx47waOH95vdSzDKBBTIAyH2b42Cq85d6IIJwf+RMuQG6yOZJmQfk+yrseHNLyYyKkPr+Ng0narIxlGvkyBMBxia/w/1Jp7FxdwI+Wen/Bp2SH/lZyc/38GsqP3TGqkHUU/vZE929ZaHckw8mQKhFHstsRF4j2/P2fxJO3en2nYtJ3VkUqN1p37cPCu73HnPJ5f38yujbFWRzKMXJkCYRSrHetWUmfB3ZyWqsiwxdRr3NLqSKVOU79uJN/9E+m4UHnOXezaHG91JMOwyxQIo9gkbVtHte/7k0JF5N6fqNOomdWRSq1Gzf05O2g+AJ6z7mDP1gRrAxmGHaZAGMXiYNJ2XL++A1fSODdwLvV8ne8BuOLm08KfMwPm4UoaFSPuMNckjFLHYQVCRBqKyBIR2Sgi60XkUTvLPC0i8bbXOhFJE5EatnmJIrLWNi/GUTmNojtx5AApn99GVU3myB2zzAXpK+DTKpCT/efhRipuX9/Job07rY5kGFkceQSRCjypqq2ATsBIEWmdfQFVfUNV/VXVH3gO+FtVsz9JFGqbH+TAnEYRpJw7w/6P76JO2kF23fg5zfy7Wx2pzGncOpijd86iiiZz9rPbOHn0oNWRDANwYIFQ1f2qGmd7fxrYCNTPY5VBwCxH5TGKX3paGuv+N5hWF9ezruNk2nS5yepIZVZTv24kXj+demn72P/RHZw7c9rqSIaBZAwL7eCdiPgCS4G2qnrKznxPIAlomnkEISI7geOAAh+r6vRcth0OhAN4e3sHzp49u1AZk5OTqVLFefsGsqeobU6J+ZLeyfP42etuPAMGFGMyxyntv+eTW/7h1r1vElshgFOdn8e1gluRtlfa2+sIps1XJjQ0NDbXszSq6tAXUAWIBe7KY5kBwE85ptWz/fdqIAHokd++AgMDtbCWLFlS6HXLqqK0edXcd1THV9OV7wzS9LS04gvlYGXh9xw1Z0rGZ/v2wCJ/tmWhvcXNtPnKADGay3eqQ+9iEhE3YC4Qoarz8lh0IDlOL6nqPtt/DwE/AB0dldO4MhtWLMY/4SXWVuxAwENflIuO90pSSP+niWpwPx1P/EzUzBetjmOUY468i0mAz4CNqvpWHst5AdcCP2abVllEqma+B24A1jkqq1FwB/Zso86v4ex3rYPPQ9/j5l7R6khOKWTYVGKr9qLzzveJW/yF1XGMcsqRf/p1BQYDvbLdynqTiIwQkRHZlrsT+E1Vz2Sb5g0sE5EEYBWwSFXNaPAWSzl3htNfDsRdL6IDvqFa9ZpWR3Ja4uJCm4e/ZlOFVrSOepotcX9bHckohyo4asOqugyQAiw3A5iRY9oOwM8hwYxC0fR01n58P8GpW1nd9QMCWvhbHcnpeVSqTO0HvufYR6HUWDCEA7X+ME+nGyXKnDw2CmTVd1MIPrGYFQ2HE3DDPVbHKTdqejfgQv9ZVNTznPmyv7n91ShRpkAY+doU/QcdNkwhvlInQoZOsTpOuePbKogd175D49SdrP/4PjQ93epIRjlhCoSRp5NHD1J90QgOu9SicXgELq6uVkcql/x6DWSl74MEnfqdlXNeszqOUU6YAmHkStPT2fnZvdTQY5y57VO8rqpldaRyLWTIq6z27ELQpqmsX/6z1XGMcsAUCCNXK2e9jP/ZFcS1eIJmAT2sjlPuubi60vTBCPa51qXubw9yYM82qyMZTs4UCMOuzTF/EbjlXVZ7diVk4PNWxzFsqnrVQAd8g7te5NSXgzifctbqSIYTMwXCuMzJY4eptjCcI1KDJsO/NE9KlzI+LfzZ0uUNmqduYfXnl/WibxjFxvzLNy6h6els++IBaukxTt0yHa8ata2OZNjR4cbBRF3dn06HviXu15lWxzGclCkQxiViF04n8PQSYhqH0yKol9VxjDx0uP89tlRoTtMVY9i3c5PVcQwnZAqEkeXA7q00j3uJTW6tCQ6baHUcIx/uFT2oEpZx9JAcMZgL51MsTmQ4G1MgDADSUlM5+vV9uGo6VQd9TgU3d6sjGQVQr3FLtnWeQvPULcR99ojVcQwnYwqEAUD0rIm0ubCWDX7PU79JK6vjGFegw42Diardj06H5rD6t6+tjmM4EVMgDLavjaLDtvdZXbkbQXeMtjqOUQgB909jm+s1+C5/lsP7Eq2OYzgJUyDKufMpZ3H5IZxTUhXfoZ+aW1rLqIoenrj1/xwPPc/+r4aRnpZmdSTDCZhvg3IububzNE7fxd4er3NV7bpWxzGKwKeFP2vaPkP7lFhWzXnV6jiGE3DkiHINRWSJiGwUkfUictkTPSLSU0ROZhtQaFy2eb1FZLOIbBORZx2Vszw7tX8LwUlfEu11I369BlodxygGHf/7BKs9uxCw+V12rl9pdRyjjHPkEUQq8KSqtgI6ASNFpLWd5f5RVX/bayKAiLgCHwB9gNbAoFzWNQohIiKCZk2b0HLTuxzTKmz0ut7qSEYxERcXfIZ+ymk8Of/NEGKiV+Hr60tERITV0YwyyGEFQlX3q2qc7f1pYCNQv4CrdwS2qeoOVb0AzAZud0zS8iUiIoLw8HCGda5Bc0li1D/VeGj0Y+YLxIks/j2SB5d60dLtAM3P/MuuXbsIDw83v2PjipXINQgR8QUCAHvHvJ1FJEFEFotIG9u0+sCebMskUfDiYuRh7NixtGxUg6ea7OAvDeK7yA2cPXuWsWPHWh3NKCZjx45lXuQ6Pt7XjNuI5OYuLczv2CgUUVXH7kCkCvA3MElV5+WYVw1IV9VkEbkJeFdVm4lIP+BGVR1uW24w0FFVL7sHU0TCgXAAb2/vwNmzZxcqZ3JyMlWqVCnUumVJTPQqbk7+mtoc59drxnPg8MmseYGBgRYmKxnl4fccGxub8SY9lTvPzMBTU4ioHA6u7uZ37KSK0ubQ0NBYVQ2yO1NVHfYC3IBfgScKuHwiUAvoDPyabfpzwHP5rR8YGKiFtWTJkkKvW5a8Pqyr6vhqOqhXG506daoCCqiPj4/V0UpEefg9+/j4ZP1e3530jKaO89KZDwea37ETK0qbgRjN5TvVkXcxCfAZsFFV38plmTq25RCRjmSc8joKRAPNRKSxiLgDA4EFjspaXuzZmsCjDTbx08mmzPprfdZ0T09PJk2aZGEyozhNmjQJT09PAC5WrM07u5pyT+2tvPxYmMXJjLLGkdcgugKDgV7ZbmO9SURGiMgI2zJ9gXUikgBMAwbailoqMIqMo4+NwLequt7eToyC0fR0Tn03mvPixrE2w/Hx8QHAx8eH6dOnExZmvjycRVhYGNOnT8/6HX+87DhbaUjvExGcPHbY4nRGWVLBURtW1WWA5LPM+8D7ucz7GTAD7xaTmAX/I/hCAivbvMi9/R7i3uEPERkZSWJiotXRDAcICwsjLCyMyMhItmzbzraEZVSfdxvxXz5M8OPfWR3PKCPMk9TlwPHD+2kaPzmjG++7Hrc6jmGBpn7diGk0jOCTv5kO/YwCMwWiHNj29WNU0bNUvHMaLq6uVscxLBI0eBLbXZvQaPnzHD+83+o4RhlgCoSTW7dsAcEnfyGm/j00bh1sdRzDQm7uFZE7P6KqJrPjq4etjmOUAaZAOLGUc2fw+vMZ9oo3AfeYztsMaNI2hFjfBwg8/Rfxf8yyOo5RypkC4cTiI16koe7jWOjreHiWrweHjNwF3v0SO118qbfseU6fPGZ1HKMUMwXCSe3aFEeHPTOIqXYd7XrcaXUcoxRxr+jBxVvepaYeZ8NX5qYFI3emQDghTU/n9A+PkyIe+N79jtVxjFKoeYeeRNcZQMjR+WyI+sXqOEYpZQqEE4r75Qvano9nY6tHqFWnodVxjFKq/eAp7BNvqvz2BCnnzlgdxyiFTIFwMmdOn6DhqlfY7tqEoP8+ZXUcoxTzrOLF0dDXaZS+l9Vfm55ejcuZAuFk1nzzIldzjIs3TsG1gsMelDecRLsedxLt1ZugpK/Ysc6MQGdcyhQIJ7J7SzyB+yKI9upNy45mlDijYJoNfpdTUoW0+SNJS021Oo5RipgC4SQ0PZ2Tcx8nRSrSeNBUq+MYZUj1WnVIDB5Hs9StRM8xvfoa/88UCCex+revaHc+jg0tRpsL08YV69BnGPGenWm/5QP2JW62Oo5RSpgC4QTOJp+kXtTL7HDxJaivuTBtXDlxcaHOwPcAODTnETQ93eJERmlgCoQTSPjmRepwhPM3TqGCm7vVcYwyqk6jZqxp9jD+56KI/+Mbq+MYpYApEGXcnq0JBO79mmivG2gVcqPVcYwyLrD/c+xw8aXe8nEknzpudRzDYo4ccrShiCwRkY0isl5EHrWzTJiIrLG9louIX7Z5iSKy1jYSXYyjcpZlmp7O8blPcB53Gg960+o4hhNwc6/IhT5vUluPsS7iOavjGBZz5BFEKvCkqrYCOgEjRaR1jmV2AteqanvgZWB6jvmhquqvqkEOzFlmJfw1h/YpMaxvMZJadRpZHcdwEi2DryO61m0EHZjD9jXLrY5jWMhhBUJV96tqnO39aTLGlq6fY5nlqpp5HBsFNHBUHmdz4XwKNf+dyC6XBgT2HWN1HMPJtAx7k1NShdQFj5GelmZ1HMMioqqO34mIL7AUaKuqp3JZ5imgpaoOt/28EzgOKPCxquY8ushcLxwIB/D29g6cPXt2oTImJydTpUrZ6RI7OWEutxz/ih8bjcWrScfCbaOMtbk4lLc2F6W9Jzf9ye0HprGg1gNUa3tLMSdznPL2O4aitTk0NDQ217M0qurQF1AFiAXuymOZUDKOMGpmm1bP9t+rgQSgR377CgwM1MJasmRJodctaUcO7NFT47w1/rX/FGk7ZanNxaW8tbko7U1PS9O1r/bQk+Pr6OH9u4svlIOVt9+xatHaDMRoLt+pDr2LSUTcgLlAhKrOy2WZ9sCnwO2qejRzuqrus/33EPADULg/k53QtjnP4cEFqt9pnpg2HEdcXKj233fx0AskfvOY1XEMCzjyLiYBPgM2qupbuSzTCJgHDFbVLdmmVxaRqpnvgRuAdY7KWpZsXxtF0NGfiPX+Lz4t/K2OYzi5Rs39iW00lKBTf7B26Y9WxzFKmCOPILoCg4FetltV40XkJhEZISIjbMuMA2oC/8txO6s3sExEEoBVwCJVLfejmmh6Oik/Pc0pqUKrgWaMaaNkBNw9kSSpQ/Ulz3I+5azVcYwS5LD+oFV1GSD5LDMcGG5n+g7A7/I1yrf432cScGENK1uPJaRGbavjGOWER6XKHO/5Gu2W3EfU7FfoNNT8cVJemCepy4iUc2fwjnqFnS4+BN75mNVxjHKm3bV3EVe5O347P2H/LtOZX3lhCkQZsfrbSdTTQ5wJfdn0t2RYot6AtwHY/+2TFicxSoopEGXAkX278NvxKas9u9C2++1WxzHKqTqNmpHQ+AE6nPmHNZFzrY5jlABTIMqAnd8+QwVSqX3XG1ZHMcq5gIEvsEfqcdXfL5gL1uWAKRCl3NbVSwk+sZi4ugNp0LSt1XGMcq6ihyfHe7xMQ91H3JxXrI5jOJgpEKWYpqeT+vMzHMWLNoPMP0ajdGgf2pfVlbvht+NTDuzeanUcw4FMgSjFYhd/RquLG9je7nGqetWwOo5hZKnT/y0EZd+cJ6yOYjiQKRCl1Lkzp2kQPZntrk0IvH201XEM4xJ1fVoQ3/h+OpxZytqlP1gdx3AQUyBKqfg5L2cMI3rdq7hWcNjzjIZRaB0GjiNJ6lJ9yfPmgrWTKlCBEBEXEQkQkZtFpJeIeDs6WHl2MGk7frtmEFelB60797E6jmHYVdHDk2M9XrFdsJ5kdRzDAfIsECJyjYhMB7YBk4FBwMPA7yISJSL3iYg5Cilme74dgyvp1Olrems1Srf2oX1Z7dkVvx2fcGDPNqvjGMUsvy/3V4CvgWtU9UZVvUdV+2rGEKG3AV5kdMhnFJNNMX8SdOoP4hrcQz3fFlbHMYx8efd/23bB+nGroxjFLM8CoaqDVHWpbVCJnPMOqeo7qvql4+KVL+lpabj88iyHuYr2AydYHccwCqSebwvifYfRIXmp6RLcyeR59VNEeuQ1X1WXFm+c8i124ccEp24h2n8SwVWrWx3HMAosYOA4kqb8iFfkc1wIuRH3ih5WRzKKQX63xzxtZ5qS0RV3A8C12BOVU2eTT+Kz+g22VGhO4K0PWR3HMK6IR6XKHO3xMn5/DydqziQ6DXnZ6khGMcjvFNOt2V/A64AbsB+4owTylRsJcyZyNcdIv+FVXFxN3TXKHr/Qfqz27EL77R9zMGm71XGMYlDQ21z/IyKRwMvAW6raSVV/ymedhiKyREQ2ish6EXnUzjIiItNEZJuIrBGRDtnm9RaRzbZ5z15hu8qUA7u3ErD7S2Kr9qJlx+utjmMYhebd721cSCdpzlNWRzGKQX63ud4sIsuBp4Cxqhqqqr8XcNupwJOq2groBIwUkdY5lukDNLO9woEPbft1BT6wzW8NDLKzrtNI+u4ZAOr3m2JxEsMomnqNW7La5z4CT//Fun/z/BvSKAPyO4L4iYxrDanAMyKyIPsrrxVVdb+qxtnenwY2AvVzLHY78JVmiAKqi0hdoCOwTVV3qOoFYLZtWaezadXvBJ3+k/iGQ6jTqJnVcQyjyAIGjmefeFPlz+e5eOG81XGMIsjvInVocexERHyBAGBljln1gT3Zfk6yTbM3PSSXbYeTcfSBt7c3kZGRhcqYnJxc6HULKz09jbr/PMMhvYrTDa8r8f1b0Warlbc2W9Xekw3u4/Y9k1n44dNU8burRPdd3n7H4Lg251kgVPXvou5ARKoAc4HHVPVUztn2dpvH9Msnqk4HpgMEBQVpz549C5UzMjKSwq5bWNE//o9Wup3oDq9x/Y03lei+wZo2W628tdmq9mqPHiRM+ZOex74lpeVT1KrTqMT2Xd5+x+C4Nud3DeInEblVRNzszGsiIhNFZFge67uRURwiVHWenUWSgIbZfm4A7MtjutM4c/rE/9/WesuDVscxjGIlLi7U6Ps27lxk5ywzhnVZld81iAeA7sAmEYkWkZ9F5C8R2QF8DMSq6uf2VhQRAT4DNqrqW7lsfwEwxHY3UyfgpKruB6KBZiLSWETcgYG2ZZ3G2szbWm98zdzWajilhk3bEdvgHoJP/sbGlb9aHccohPxOMR0AxgBjbNcR6gLngC2qml//vl3J6KdprYjE26Y9DzSybfsj4GfgJjI6AzwL3Geblyoio4BfyXgY73NVXX+ljSut9u/ajP+er4ip9h+Cgq+zOo5hOIzfoIkcmLqQir+OIbVDKBXc3K2OZFyBAg80oKqJQOIVLL8M+9cSsi+jwMhc5v1MRgFxOvu+f5bqCA36vW51FMNwKM8qXmzq9CIdoh5l5dw3CRn4nNWRjCtguuouYZtW/kbg6b+Ib2RuazXKh4AbhrC2YgCtNr3H0YNJVscxroApECUoPS0N19+f5xA18Bswzuo4hlEixMWFane9QyVNYfsse927GaVVfncxldy9aeVA7E8f0Sx1K7s7jMGzipfVcQyjxPi08Ce27kA6nviZzTF/WR3HKKD8jiDmZ74RkbmOjeLczpw+gW98xm2tHW4OtzqOYZS4tndP4hA1cFn8NGmpqVbHMQogvwKR/SJzE0cGcXZrZo+nNsfNba1GuVWl2lXsDnqOZmnbiPnhHavjGAWQX4HQXN4bVyBp2zoCk74m2usGWprbWo1yLPCm4ax3b0eL9W9z4sgBq+MY+civQPiJyCkROQ20t70/JSKnRSRntxlGLo7OfYILuNF4wFSroxiGpcTFBc873qaKnmXzrDFWxzHykd+AQa6qWk1Vq6pqBdv7zJ+rlVTIsizhr9n4nVvJumYPUquej9VxDMNyjVsHE+Pdj+AjC9ga/4/VcYw8mNtcHeh8yllq/jOB3S716dDPPCBkGJla3/0ax8SL9IVPkp6WZnUcIxemQDhQ3JxXaKD7OXntK2YQd8PIplr1muz0H0OL1M3E/vi+1XGMXJgC4SAH9mzDb8enrK7cjXbXlmx/+IZRFgTd9hAb3VpzzZqpnDx22Oo4hh2mQDjI3m+fwoV0vPvl1pGtYZRv4uKC+61v4qWn2fTNM1bHMewwBcIB1v+7iMDTS1jdaCj1fFtYHccwSq1r2nchpvadBB2ex/a1UVbHMXIwBaKYpV68gOefz7Of2gQMmmB1HMMo9VrePYVTUpULC55A09OtjmNkYwpEMYv5/g0apydyoPOLeHhWsTqOYZR6XjVqs7Xdk7S6uJ7YhR9bHcfIxmEFQkQ+F5FDIrIul/lPi0i87bVORNJEpIZtXqKIrLXNi3FUxuJ2eF8ibTe9xxqPIPyvH2x1HMMoM4LuGM2WCs3xjZvM6ZPHrI5j2DjyCGIG0Du3mar6hqr6q6o/8Bzwt6pm/z8j1DY/yIEZi9Xubx7FjVRq9p+GuJiDM8MoKBdXV+TmN6mhJ1n/jXlmqLRw2LeYqi4FCvqnwCBglqOylIQ1S74nMDmSON/7qd+kjdVxDKPMaRbQg+iatxJ04FsSN5aZEwdOzfI/c0XEk4wjjezdiSvwm4jEikip7xs75WwyNZc+n/HE9KDxVscxjDKr+aApJIsnZ3543FywLgUkY1hoB21cxBdYqKpt81hmAHCPqt6abVo9Vd0nIlcDvwOjbUck9tYPB8IBvL29A2fPnl2orMnJyVSpUriLyimxX9H79Fx+bPwSXj7+hdqGFYrS5rKqvLW5LLb31LqF3HbkE36s+wReLa694vXLYpuLqihtDg0Njc31VL6qOuwF+ALr8lnmB+DuPOZPAJ4qyP4CAwO1sJYsWVKo9RI3rdbz467SVW/1LfS+rVLYNpdl5a3NZbG9qRcv6taJAXpwvK8mnzp+xeuXxTYXVVHaDMRoLt+plp5iEhEv4Frgx2zTKotI1cz3wA2A3TuhrKbp6STPfYQU8aDJ3W9bHccwnIJrhQqk9p7C1RxjzTdjrY5TrjnyNtdZwAqghYgkicj9IjJCREZkW+xO4DdVPZNtmjewTEQSgFXAIlX9xVE5iyJmwYe0uZDAxjZPUNO7gdVxDMNptAy+jujqfQjcN4tdm+OtjlNuVXDUhlV1UAGWmUHG7bDZp+0A/ByTqvgcObCb5vGT2OjWmuC7Hrc6jmE4nSaDppLyYTCn5j2OPvOnuXXcAuYTL6TdM0fioReo3O9DM8a0YThATe8GbGgxmnbn41j921dWxymXTIEohLhfZtDhzFLimjxIo+b+VscxDKcV1Pcpdrj4Uj9qIsmnjlsdp9wxBeIKnTx6kEZR49nmeg1Bg8ZZHccwnFoFN3cu9HmT2nqMdV+bLsFLmikQV2jLV49QXU8ht7+Pm3tFq+MYhtNrGXwd0bVuJ/jgt2xLWGZ1nHLFFIgrsCZyLsEnfyG6wRCuad/F6jiGUW60vOdNjosXuuBR0lJTrY5TbpgCUUAnjx2mTuRT7HJpQMA9r1odxzDKFa+rarGr4ziapW0j+rvXrY5TbpgCUUBbv3iQq/QkF279EI9Kla2OYxjlTofe95HgEUy7TdM4mLTd6jjlgikQBRC76FOCTv9JjO8DNAvoYXUcwyiXxMWF2gPex4V09s16xOo45YIpEPk4tHcnTaPHsblCC4LvednqOIZRrtVr3JL4a0YQcGYZq3/72uo4Ts8UiDxoejoHZt6Pm6biOeATKri5Wx3JMMq9oIEvsNPFh3rLx5lnIxzMFIg8rJwzmfYpsaxt8xQNm5X63j8Mo1xwc6/I+T5v4c1R1kU8a3Ucp2YKRC62JSyjw6Y3SagUQse+T1kdxzCMbFoGX8fKmncQfGCOeTbCgUyBsOP0yWN4zL+fE1KNRsO+NJ2EGUYplPlsBAseIfXiBavjOCXzzZeDpqez5dNh1Ek/xNE+H3FV7bpWRzIMww6vq2qxu9NLNE3bTvQ3L1kdxymZApHDqrlvEXh6CdFNHqJVyI1WxzEMIw8deg8lrnIPOuz42Iwb4QCmQGSzJS4S/3WTWeMRSIi5pdUwyoRGgz/gnFTk3PcPmW44ipkjR5T7XEQOiYjd4UJFpKeInBSReNtrXLZ5vUVks4hsExGH3qYQERGBr68vsSv/odr8IRzBi4b3R5gxHgyjjKhVpxFbA8bS8uIGXh4WSmxsLL6+vkRERFgdrcxz5BHEDKB3Psv8o6r+ttdEABFxBT4A+gCtgUEi0toRASMiIggPD2ff3r1cf/YHqsk5+v3ows+//eWI3RmG4SCbT1Xjz+TGjGm8CZeLp9m1axfh4eGmSBSRwwqEqi4FjhVi1Y7ANlXdoaoXgNnA7cUazmbs2LGcPXuW6cP9aC/beXBFHVYkbGPsWDNQumGUJS+8+AL3f5NEGi50P5cxhP3Zs2fNv+UiElV13MZFfIGFqtrWzryewFwgCdgHPKWq60WkL9BbVYfblhsMhKjqqFz2EQ6EA3h7ewfOnj27wPliY2MhLYWBZz5ntUcIm90Ds+YFBgbmsaZzSE5OpkqVKlbHKFHlrc3lpb2xsbEA1DqzgXvTv+cLl34cq9wKMP+W8xMaGhqrqkF2Z6qqw16AL7Aul3nVgCq29zcBW23v+wGfZltuMPBeQfYXGBioV8LHx0cBrVGtsr4xZYoCCqiPj88VbaesWrJkidURSlx5a3N5aW/mv2VBNHZcsJ4c562N6tYw/5YLAIjRXL5TLbuLSVVPqWqy7f3PgJuI1CLjiKJhtkUbkHGEUewmTZqEp6cnx06dyXoYztPTk0mTJjlid4ZhOEjmv2VF+bvSTbiRxhd3N+SVl1+xOlqZZlmBEJE6IiK29x1tWY4C0UAzEWksIu7AQGCBIzKEhYUxffp0fHx8APDx8WH69OmEhYU5YneGYThI9n/L6W5VeXVHC3pV2UkzjwNWRyvTHHmb6yxgBdBCRJJE5H4RGSEiI2yL9AXWiUgCMA0YaDviSQVGAb8CG4FvVXW9o3KGhYWRmJhIYGAgiYmJpjgYRhmV/d/ySzP+Zm3FDrRbN4WkbXbvtDcKoIKjNqyqg/KZ/z7wfi7zfgZ+dkQuwzCcn4urK1cP/pTUT7uRPOcB0p75B9cKDvu6c1rmSWrDMJySd4Nr2NJhPC0vbmCV6aupUJy+pF68eJGkpCRSUlLyXM7Ly4uNGzeWUKrSwbTZ+RW2vR4eHjRo0AA3NzcHpCo5gbeEE7dlEYHbP2DHulto0jbE6khlitMXiKSkJKpWrYqvry+2a+J2nT59mqpVq5ZgMuuZNju/wrRXVTl69ChJSUk0btzYQclKhri40Pje6Zz6oBPMC+d80xVU9PC0OlaZ4fSnmFJSUqhZs2aexcEwjP8nItSsWTPfo+6y4qradUnq/jpN0hOJ+/Jpq+OUKU5fIABTHAzjCjnbvxn//wxkZY3bCNkXwdqlP1odp8woFwXCMAyj/f3/Y49rA+r+9QhHDyZZHadMMAWiBDhLXzg9e/YkJiYmz2Xmz5/Phg0bHJ4lMTGRtm0v6+LrsmW++eabrJ9jYmJ4+uniOcXg6+tLu3bt8PPz44YbbuDAgdwfyIqMjOSWW24plv1eqc8++4yvvvoqz2UmTJjA1KlTSyiRdSpVrkraXZ9RVc+Q9MVQ0tPSrI5U6pkCUQapKunp6VbHsKswBSLVQYO85CwQQUFBvPHGG8W2/SVLlpCQkEBQUBCvvvpqsW23ON1///0MGTLE6hilRpO2IcS3fgq/lGhWzTZd6uTH6e9iyu6ln9azYd8pu/PS0tJwLcQgQa3rVWP8rW2ueL2ffvqJV155hQsXLlCzZk0iIiLw9vbm8OHD3H333Rw9epTg4GB++eUXYmNjSU5Opk+fPoSGhrJixQrmz5/P5MmTiY6O5ty5c/Tt25eXXsq419vX15d7772Xn376iYsXL/Ldd9/RsmXLy7a9ePFi4uLiqFWrFl9//TXTpk3jwoULhISE8L///S/Pz6NKlSo8+uijLFy4kEqVKvHjjz+yfft2FixYwN9//80rr7zC3LlzARg5ciSHDx/G09OTTz75hJYtWzJ06FBq1KjB6tWr6dChA1WrVmX79u3s3buXPXv2MGbMGB544AFUlTFjxrB48WJEhBdeeIEBAwZckiUxMZHBgwdz5swZAN5//326dOnCs88+y8aNG/H39+fee+8lICCAyZMn88svv3Ds2DGGDRvGjh078PT0ZPr06bRv354JEyawe/duduzYwe7du3nsscd45JFH8vxd9ujRg2nTppGSksJDDz1ETEwMFSpU4K233iI0NDRrufT0dFq0aMHy5cupXbs26enpNG/enKioKJ566imqVatGTEwMBw4cYMqUKfTt2zfX9kdGRjJ+/Hi8vb2Jj4/nrrvuol27drz77rucO3eO+fPnc8011/Dqq69Ss2ZNnnrqKT755BOmT5/OhQsXaNq0KTNnzsTTs/zd0dOx3xhWT/2bDlveYVtCL5r6dbM6UqlljiAs0q1bN6Kioli9ejUDBw5kypQpALz00kv06tWLuLg47rzzTnbv3p21zubNmxkyZAirV6/Gx8eHSZMmERMTw5o1a/j7779Zs2ZN1rK1atUiLi6Ohx56KOv0Qc5t79mzB4CNGzcyZ84c/v33X+Lj43F1dc13oJUzZ87QqVMnEhIS6NGjB5988gldunThtttu44033iA+Pp5rrrmG8PBw3nvvPWJjY5k6dSoPP/xw1ja2bNnCH3/8wZtvvgnAmjVrWLRoEStWrGDixIns27ePefPmER8fT0JCAn/88QdPP/00+/fvvyTL1Vdfze+//05cXBxz5szJ+kKfPHky3bt3Jz4+nscff/ySdcaPH09AQABr1qzh1VdfveSv7E2bNvHrr7+yatUqXnrpJS5evJjnZ7Fw4ULatWvHBx98AMDatWuZNWsW99577yV3Arm4uHDPPfdkfbZ//PEHfn5+1KpVC4D9+/ezbNkyFi5cyLPPZgykmFf7ExISePfdd1m7di0zZ85ky5YtrFq1iuHDh/Pee+9dlvOuu+4iOjqahIQEWrVqxWeffZZnu5yVuLjQeNgXnBAvKs5/gORTx62OVGqVqyOIvP7SL+n745OSkhgwYAD79+/nwoULWfebL1u2jB9++AGA3r17c9VVV2Wt4+PjQ6dOnbJ+/vbbb5k+fTqpqans37+fDRs20L59eyDjywAy+sKfN2+e3W1Xr14dgD///JPY2FiCg4MBOHfuHFdffXWe+d3d3bPOqwcGBvL7779ftkxycjLLly+nX79+WdPOnz+f9b5fv36XHKXcfvvtVKpUiUqVKhEaGsqqVatYtmwZgwYNwtXVFW9vb6699lqio6Oz2gkZD0OOGjUqq7ht2bIlz+yZn0XmEU6vXr04evQoJ0+eBODmm2+mYsWKVKxYkauvvpqDBw/SoEGDy7YRGhqKq6sr7du355VXXuG+++5j9OjRALRs2RIfH5/LsgwbNozbb7+dxx57jM8//5z77rsva94dd9yBi4sLrVu35uDBg1k57bW/WrVqBAcHU7duXQCuueYabrjhBgDatWvHkiVLLsu7bt06XnjhBU6cOEFycjI33nhjvp+Ts6peqw7rb/iAlr/eTfz0oXR44oesHp2N/1euCkRpMnr0aJ544gluu+02IiMjmTBhAkDmGBh2Va5cOev9zp07mTp1KtHR0Vx11VUMHTr0kr9WK1asCICrq2vWOf7ctq2q3Hvvvbz22msFzu/m5pZ1K2T2fWSXnp5O9erViY+Pz7c9cPmtlSKS5+eR6e2338bb25uEhATS09Px8PDIdx17283cf+ZnB7m3DTKuQWT+9Z/bNnNq2LAh3t7e/PXXX6xcufKSI7Xs+83cVl7bzL68i4tL1s8uLi52Mw8dOpT58+fj5+fHjBkziIyMzDevM2vT5SZWbBtN5x3TiJr9Kp3ufsHqSKWOKZkWOXnyJPXr1wfgyy+/zJrerVs3vv32WwB+++03jh+3f/h76tQpKleujJeXFwcPHmTx4sX57jPntk+cOAHAf/7zH77//nsOHToEwLFjx9i1a1eh2lW1alVOnz4NQLVq1WjcuDHfffcdkPFll5CQkOu6P/74IykpKRw9epTIyEiCg4Pp0aMHc+bMIS0tjcOHD7N06VI6dux4yXonT56kbt26uLi4MHPmTNJsd6dkz5JTjx49sr6cIyMjqVWrFtWqVStUm+1tc8uWLezevZsWLVpcttzw4cO555576N+/f77XvQrS/oI6ffo0devW5eLFi2asZptO97zEas8uBG5+i00rf7M6TqljCkQJOHv2LA0aNMh6vfXWW0yYMIF+/frRvXv3S/4KHT9+PL/99hsdOnRg8eLF1K1b1+6pLz8/PwICAmjTpg3Dhg2ja9eu+ebIue06depQtWpVWrduzSuvvMINN9xA+/btuf766y87z19QAwcO5I033iAgIIDt27cTERHBZ599hp+fH23atOHHH3N/SKljx47cfPPNdOrUiRdffJF69epx55130r59e/z8/OjVqxdTpkyhTp06l6z38MMP8+WXX9KpUye2bNmSdWTSvn17KlSogJ+fH2+//fYl60yYMIGYmBjat2/Ps88+e0mRLqyHH36YtLQ02rVrx4ABA5gxY8Ylf+Vnuu2220hOTr7k9FJuCtL+gnr55ZcJCQnh+uuvp2XLloXahrMRFxeaPDCTgy61qbH4QY4c2GN1pNIlt6HmyuLL3pCjGzZsyHfIPVXVU6dOFWg5R0tJSdGLFy+qqury5cvVz8/PYdtu165dsW27qMaPH69vvPGGw/dTGn7P0dHR2q1btxLZV1HaW9B/O6VNYYbf3L42Ss+Oq6XrJnXTixfOF38oB3PUkKPmGkQps3v3bvr37096ejru7u588sknDtv2tGnTim3bRsFMnjyZDz/80JziKWWatA0hetsEguOfJ+rTR+j00EdWRyoVHFYgRORz4BbgkKpe9siriIQBz9h+TAYeUtUE27xE4DSQBqSqapCjcpY2zZo1Y/Xq1SWy7dzOz1sh8yK9s3v22WezbmE1SpfgO0ayMimWTgdnseqH1nS8M+/nX8oDR16DmAH0zmP+TuBaVW0PvAxMzzE/VFX9y1NxMAzDWoEPfsTaigH4x09g48pfrY5jOYcVCFVdChzLY/5yVc28RScKuPxGc8MwjBJUwc2dRg9+x0EXb7wXD2df4marI1lKtAD3bhd64yK+wEJ7p5hyLPcU0FJVh9t+3gkcBxT4WFVzHl1kXzccCAfw9vYOnD179iXzvby8aNq0ab5ZC9vVRllm2uz8itLebdu2ZT08WJYkJycXuYPM5CN76L52DIelJjs6TcG9lA8yVJQ2h4aGxuZ2psbyi9QiEgrcD2TvEKWrqu4TkauB30Vkk+2I5DK24jEdICgoSHv27HnJ/I0bNxboCenyNtIYmDaXB0Vpr4eHBwEBAcWcyPEiIyPJ+T1QGGuvrkKrP4dxNv4tWj2xCDf3y29ZLi2Kq805WfochIi0Bz4FblfVo5nTVXWf7b+HgB+Awj0ZVEqICE8++WTWz1OnTi3yRdn4+Hh+/vnnIib7f76+vhw5ciTPZWbMmMG+ffuKbZ+5KUj32Dnbv2DBAiZPnlws+3d1dcXf35+2bdvSr18/zp49m+uyM2bMYNSoUcWy3ys1btw4/vjjjzyXGTFiBN9//30JJXIu7XrcSVy7F2mfEk38B0PQUtqDsiNZViBEpBEwDxisqluyTa8sIlUz3wM3AOusSVk8KlasyLx58/L9Ar4SxV0gCqIwBcJRXXnnbP9tt91WbHcHVapUifj4eNatW4e7uzsffVQ6b3mcOHEi1113ndUxnFrHvk+wouEDBJ/8hajPHs9/BSfjyNtcZwE9gVoikgSMB9wAVPUjYBxQE/ifrQ+czNtZvYEfbNMqAN+o6i/FEmrxs3Bgrd1ZldJSwbUQH0eddtAn779cK1SoQHh4OG+//TaTJl3aB/3hw4cZMWJEVq+t77zzDl27dqVdu3b8888/eHl5UatWLd5++22GDBnC4MGDGTx4MOPGjePcuXMsW7aM5557juuvv75Yuq9OTEykT58+dOvWjeXLl1O/fn1+/PFHFi1aRExMDGFhYVSqVIkVK1awYcMGnnjiCZKTk6lVqxYzZsygbt269OzZky5duvDvv/9y2223sXbtWjw8PFi/fj0HDx7krbfe4pZbbsm3e2yAVatW8dhjj3Hu3DkqVarEF198QePGjS9r/7lz54iJieH9999n165dDBs2jMOHD1O7dm2++OILGjVqxNChQ/Hw8GDNmjWXdKmdl+7du7NmzZpcuwfPdPr0adq3b8+WLVtwc3Pj1KlTtG/fnq1bt3L99dcTEhLCkiVLOHHiBJ999hndu3fPtf0zZsxg/vz5pKWlsW7dOp588kkuXLjAzJkzqVixIj///DM1atRg6NCh3HLLLfTt25eJEyfy008/ce7cObp06cLHH3/sdMOGWqXTfVNY9f4BOu+dwco59QgZ8Ez+KzkJR97FNEhV66qqm6o2UNXPVPUjW3FAVYer6lW2W1mzbmdV1R2q6md7tVFVpxjVY+TIkURERFx20e/RRx/l8ccfJzo6mrlz5zJ8+HAAunbtyr///sv69etp0qQJ//zzDwBRUVF06dKFiRMnMmDAAOLj4xkwYECxdl+9detWRo4cyfr166levTpz586lb9++BAUFERERQXx8PBUqVGD06NF8//33xMbGMmzYMMaOHZu1jRMnTvD3339nnVpLTEzk77//ZtGiRYwYMYKUlJR8u8eGjF5Rly5dyurVq5k4cSLPP/887u7ul7U/u1GjRjFkyBDWrFlDWFjYJQXx4MGDl3WpnZvU1FQWL15Mu3bt8vx8IaPfp549e7Jo0SIAZs+ezX//+1/c3NyytrVq1SreeeedrHE78mr/unXr+Oabb1i1ahVjx47F09OT1atX07lzZ7sjxI0aNYro6GjWrVvHuXPnWLhwYZ5tMwpOXFzo8NDnxHt2JnjDa8T89LHVkUqM5RepS1Qef+mfc/DFy2rVqjFkyBCmTZtGpUqVsqb/8ccfl4zAdurUKU6fPk337t1ZunQpPj4+PPTQQ0yfPp29e/dSo0YNu3crFEf31ZkaN26Mv78/kNGVd2Ji4mXLbN68mXXr1nH99dcDGXfLZHY9DVz2pd2/f39cXFxo1qwZTZo0YdOmTSxbtizf7rFPnjzJvffey9atWxGRfIsbwIoVK7K6OB88eDBjxozJmnfzzTdf1qV2TufOnctqf/fu3bn//vsJCQnJ9fPNNHz4cKZMmcIdd9zBF198cclT8Nm7X8/8PPNqf2hoKFWrVqVq1ap4eXlx6623AhldeWcf9yPTkiVLmDJlCmfPnuXYsWO0adMmax2j6Cq4udNi5HdsfPcm/GOeJa6COx365N+XVllXvgqExR577DE6dOhwSSdt6enprFix4pKiARm9eH7wwQfs3r2bSZMm8cMPP/D999/TvXt3u9u2d7vylXZfnSnn8ufOnbO7vzZt2rBixQq72yiurrxffPFFQkND+eGHH0hMTCzUnRrZ922vS+2cMq9BZJfX55upa9euWUdKaWlpl4yZfSXdr+fMmV9X3ikpKTz88MPExMTQsGFDJkyYcNmRmFF0lSpXxXf0T2yd1od2UU8SX6Ei/tffbXUshzK9uZagGjVq0L9//0tG8rrhhht4//33s37O/GJq2LAhR44cYevWrTRp0oRu3boxderUrAKRsytrR3RfnVP2fbZo0YLDhw9nFYiLFy+yfv36XNf97rvvSE9PZ/v27ezYsYMWLVoUqHvs7N2iz5gxw26WnLp06ULm8zARERF061b0ISUL+vkOGTKEQYMGFain1oJ2D56fzGJQq1YtkpOTzV1LDlS5anUajFrETrdraL1sNAlLvrM6kkOZAlHCnnzyyUvuZpo2bVpWt9OtW7e+5I6ZkJAQmjdvDmSc6ti7d2/Wl11oaCgbNmzA39+fOXPmOKT76pyGDh3KiBEj8Pf3Jy0tje+//55nnnkGPz8//P39Wb58ea7rtmjRgmuvvZY+ffrw0Ucf4eHhUaDusceMGcNzzz1H165ds8Z5sNf+7KZNm8YXX3xB+/btmTlzJu+++26R217QzzcsLIzjx48zaNCgfLdZ0O7B81O9enUeeOAB2rVrxx133JE1MqDhGFW9auD98M/sruBDq8gRrP61+P+tlRq5dfNaFl/O0N13SSqpNt9777363Xfflci+8uPoNn/33Xd6zz33OHQfV8J09+04J44e0o0vh2jqOC9dOW9aiewzN6a7b8Mo5UaPHs3ixYtL/PkUwxpeNWrj9tivbHjvDjomvEDUuZNON2ypKRCGw2W/duDM3nvvPasjGCXMs4oXzR9fRNx7/em05Q2iPtxD8AMf4FrBOb5azTUIwzCMIqjo4Un7x+YRdXV/Oh2czbo3byL5lP2x5MsaUyAMwzCKqIKbO50e/oSVrcfS5mw0h9/p6RRdhZsCYRiGUUxC+o9h438+p2b6IarMCGX1b19bHalITIEwDMMoRu163MnpIX9wyLUOActHEvW/cC6cL5sPLpoCUQIOHjzI3XffTZMmTQgMDKRz58788MMPAMTExOTZeV5xK0j31JGRkXk+02AYRt7qN2lDw6eXsbLWf+l0aA57pnRiW8Iyq2NdMVMgcoiIiMDX1xcXFxd8fX2znnQtLFXljjvuoEePHuzYsYPY2Fhmz55NUlISAEFBQUybNq04ohcbUyAMo+gqengSMupz4rt+SLW0E/jOu5UVnzxKyrkzVkcrMFMgsomIiCA8PJxdu3ahquzatYvw8PAiFYm//voLd3d3RowYkTXNx8cnq5O2zMFx0tPT8fX15cSJE1nLNW3alIMHD3L48GH++9//EhwcTHBwMP/++y+Q8XTvsGHD6NmzJ02aNMm10HzxxRc0b96ca6+9NmtdgMWLFxMSEkJAQADXXXcdBw8eJDExkY8++oi3334bf39//vnnH3766afLljMMo2D8r78b90eiibuqN533zuDYlADiFn9RJgYgMgUim7Fjx142etjZs2cv6cb6Sq1fv54OHTrku5yLiwu333571qmnlStX4uvri7e3d65dgkP+XXnv37+f8ePH8++///L7779f0nNsp06diIqKYvXq1QwcOJApU6bg6+vLiBEjePzxx4mPj6d79+5069btsuUMwyg4rxq16fjYLNb95yvOSyU6rHyMTa91Zd2yBaW6UDhywKDPgVuAQ6ra1s58Ad4FbgLOAkNVNc42r7dtnivwqaoWz1iS+cgctKeg0wtj5MiRLFu2DHd3d6Kjoy+ZN2DAACZOnMh9993H7Nmzs7rMzq1LcMi/K++VK1fSs2dPateunbWPzC6l9+3bx/Dhw9m/fz8XLlygcePGdjMnJSUxYMCAfJczDCNvbbvfTlrnm1k1fxpN1r1LrT8GsyWyOWeCR9E2dGCpG/fakUcQM4DeeczvAzSzvcKBDwFExBX4wDa/NTBIRFo7MGeWRo0aXdH0gmjTpg1xcXFZP3/wwQf8+eefHD58+LJlO3fuzLZt2zh8+DDz58/PGkMgs0vw+Ph44uPj2bt3b9bYFQXpyju3kcWefvppRo0axdq1a/n4449z7SJ69OjRBVrOMIz8uVaoQMe+T1DlmQ2sbPMildNOErDiEU692oyoD0ewNf4f0rN1TJmXzGumsbGxxXLNNCdHjii3FDiWxyK3A1/Z+ouKAqqLSF2gI7BNM0aWuwDMti3rcJMmTcLT0/OSaZ6enpcNE3olevXqRUpKCh9++GHWtJynsTKJCHfeeSdPPPEErVq1ombNmkDuXYIXREhICJGRkRw9epSLFy/y3Xf/3z3xqVOnsrrSzt47ac6utLN3ue2IXmINozzyqFSZkH5P4f38OuK7fcTuyu0IPPAtzebfwrGXmxD9dj9WzpnMpug/OHn04GWnoiIiInj44YdwTTuD2/nDxXLNNCcrOwypD+zJ9nOSbZq96SElESgsLAzIuBaxe/duGjVqxKRJk7KmF4aIMH/+fB5//HGmTJlC7dq1qVy5Mq+//rrd5QcMGEBwcPAl/RdNmzaNkSNH0r59e1JTU+nRo8cl3YLnpW7dukyYMIHOnTtTt25dOnTokNVt9nPPPUe/fv2oX78+nTp1YufOnQDceuut9O3blx9//JH33nuPCRMm2F3OMIyiq+Dmjv91g+C6QRw/vJ9t/85Dtv/JNSejqHHyN9iYsVyKunFCvEgXVxTok5bMyccFuMCR8xE8yv9fMy3Kd1Z2ogUY1avQGxfxBRbmcg1iEfCaqi6z/fwnMAZoAtyoqsNt0wcDHVV1dC77CCfjFBXe3t6BmQPFZPLy8qJp06b5Zk1LS8PV1bXgjXMCps3Oryjt3bZt22XDqpYFycnJdoflLWs0PZ2U00e4cGQ7rsn7cT9/DI/UUwjpCOmcuiCcoQpnXKpQoWYTDpytiLhknE4ODAws8H5CQ0NjVTXI3jwrjyCSgIbZfm4A7APcc5lul6pOB6YDBAUFac4hKTdu3FigsaZPO3hM6tLItNn5FaW9Hh4eBAQEFHMix4uMjCzU0LRlja+vL7t27QJg6tSpPD3mKSDjNnp748gXhpW3uS4AhkiGTsBJVd0PRAPNRKSxiLgDA23LGoZhGDaOuGaakyNvc50F9ARqiUgSMB5wA1DVj4CfybjFdRsZt7neZ5uXKiKjgF/JuM31c1XNfbDjAlDVXO/kMQzjco489WwUj+zXTCHjyKGo10xzcliBUNU8B+W1DXU3Mpd5P5NRQIrMw8ODo0ePUrNmTVMkDKMAVJWjR4/i4eFhdRQjH2FhYYSFhREZGVlsp5Wyc45hj/LQoEEDkpKS7D53kF1KSkq5+wdh2uz8CtteDw+PSx64NMonpy8Qbm5uBXryNzIyskxekCsK02bnV97aaxQv0xeTYRiGYZcpEIZhGIZdpkAYhmEYdjn0SeqSJiKHgV2FXL0WcKQY45QFps3Or7y1F0ybr5SPqta2N8OpCkRRiEhMbo+bOyvTZudX3toLps3FyZxiMgzDMOwyBcIwDMOwyxSI/zfd6gAWMG12fuWtvWDaXGzMNQjDMAzDLnMEYRiGYdhlCoRhGIZhV7kvECLSW0Q2i8g2EXnW6jwlQUQ+F5FDIrLO6iwlQUQaisgSEdkoIutF5FGrMzmaiHiIyCoRSbC1+SWrM5UUEXEVkdUistDqLCVBRBJFZK2IxItITLFuuzxfgxARV2ALcD0ZI9xFA4NUdYOlwRxMRHoAycBX9oaDdTYiUheoq6pxIlIViAXucObfs2T0bV9ZVZNFxA1YBjyqqlEWR3M4EXkCCAKqqeotVudxNBFJBIJUtdgfDizvRxAdgW2qukNVLwCzgdstzuRwqroUOGZ1jpKiqvtVNc72/jQZw8DXtzaVY2mGZNuPbraX0/81KCINgJuBT63O4gzKe4GoD+zJ9nMSTv7FUd6JiC8QAKy0OIrD2U61xAOHgN9V1enbDLwDjAHSLc5RkhT4TURiRSS8ODdc3guEvSHmnP6vrPJKRKoAc4HHVPWU1XkcTVXTVNUfaAB0FBGnPp0oIrcAh1Q11uosJayrqnYA+gAjbaeQi0V5LxBJQMNsPzcA9lmUxXAg23n4uUCEqs6zOk9JUtUTQCTQ29okDtcVuM12Tn420EtEvrY2kuOp6j7bfw8BP5Bx6rxYlPcCEQ00E5HGIuIODAQWWJzJKGa2C7afARtV9S2r85QEEaktItVt7ysB1wGbLA3lYKr6nKo2UFVfMv4t/6Wq91gcy6FEpLLtxgtEpDJwA1BsdyeW6wKhqqnAKOBXMi5cfquq661N5XgiMgtYAbQQkSQRud/qTA7WFRhMxl+U8bbXTVaHcrC6wBIRWUPGH0K/q2q5uO2znPEGlolIArAKWKSqvxTXxsv1ba6GYRhG7sr1EYRhGIaRO1MgDMMwDLtMgTAMwzDsMgXCMAzDsMsUCMMwDMMuUyAMw0FEpJKI/G3rFLKg64wSkfscmcswCsrc5moYDiIiI4EKqvruFazjCfyrqgGOS2YYBWOOIAzjColIsIissY25UNk23oK9fo7CgB9t6/S0HU18KyJbRGSyiITZxmxYKyLXAKjqWSBRRIqtuwTDKKwKVgcwjLJGVaNFZAHwClAJ+FpVL+newNZ1SxNVTcw22Q9oRUZX6zuAT1W1o20Ao9HAY7blYoDuZDwZaxiWMQXCMApnIhldWKQAj9iZXws4kWNatKruBxCR7cBvtulrgdBsyx0CWhZnWMMoDHOKyTAKpwZQBagKeNiZf87O9PPZ3qdn+zmdS/9Y87CtbxiWMgXCMApnOvAiEAG8nnOmqh4HXEXEXvHIT3OKsUdOwygsUyAM4wqJyBAgVVW/ASYDwSLSy86ivwHdCrGLrsAfRYhoGMXC3OZqGA4iIgHAE6o62JHrGIajmCMIw3AQVV1NxpgMBX5QjoyL2y86KJJhXBFzBGEYhmHYZY4gDMMwDLtMgTAMwzDsMgXCMAzDsMsUCMMwDMMuUyAMwzAMu/4PF+gGDfOpNVkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Sample data\n",
    "xdat = np.array([0,1,2,3,4,5])\n",
    "fdat = np.array([1.5,2.5,3,2.5,1.5,1])\n",
    "x_cont_dat = np.linspace(min(xdat),max(xdat),1000)\n",
    "\n",
    "# Calculate interpolating polynomial \n",
    "f_cont_dat_lagrange =lagrange_interp(xdat, fdat, x_cont_dat)\n",
    "f_cont_dat_newton, coeff = newton_interp(xdat, fdat, x_cont_dat)\n",
    "\n",
    "# Plot\n",
    "plt.figure(1)\n",
    "plt.title('Interpolation Result')\n",
    "plt.scatter(xdat,fdat,color = 'black',label='Given data')\n",
    "plt.plot(x_cont_dat,f_cont_dat_lagrange,label='Lagrange Interpolation Polynomial')\n",
    "plt.plot(x_cont_dat,f_cont_dat_newton,label='Newton Interpolation Polynomial')\n",
    "plt.ylabel('F (N)')\n",
    "plt.xlabel('x (m)')\n",
    "plt.grid()\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}