{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Von Karman turbulence spectrum\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Total phase variance for Kolmogorov spectrum\n", "\n", "We compute the Kolmogorov spectrum total variance over a 10m telescope and we compare it with Noll('76) formula ($\\Delta_{1}$ = 1.029 ($\\frac{D}{r0})^{5/3}$)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Variance computed using the VonKarmanPsd class: 2213\n", "Variance from Noll's formula: 2216\n" ] } ], "source": [ "import numpy as np\n", "from scipy.special import jv\n", "from arte.atmo.von_karman_psd import VonKarmanPsd\n", "\n", "'Compute the total variance of Kolmogorov over a 10m telescope:'\n", "R = 5\n", "r0 = 0.1\n", "L0 = np.inf\n", "freqs = np.logspace(-8, 4, 1000)\n", "bess = jv(1, 2*np.pi*R*freqs)\n", "\n", "psd = VonKarmanPsd(r0, L0)\n", "psd_piston_removed = psd.spatial_psd(freqs) * (1 - (bess/(np.pi*R*freqs))**2)\n", "var_in_square_rad = np.trapz(psd_piston_removed*2*np.pi*freqs, freqs)\n", "\n", "noll_var_in_square_rad = 1.029*(2*R/r0)**(5./3)\n", "\n", "print(\"Variance computed using the VonKarmanPsd class: %d\" %(var_in_square_rad))\n", "print(\"Variance from Noll's formula: %d\" %(noll_var_in_square_rad))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Kolmogorov and Von Karman spectra\n", "\n", "Plot spectra for different outer scale values" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEQCAYAAACX5IJuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd1xV5R/A8c9zLxsEAUVRUFBwIYpbc6/cK3eWVlqWZe6RmplpWjmyXJlZ+jO1NFfmyFGuHLhy4EIcoGwZyr73Pr8/UAK5KDIE9Hm/Xrzinvs9534PqV+e85zzfYSUEkVRFEXJLU1BJ6AoiqI8H1RBURRFUfKEKiiKoihKnlAFRVEURckTqqAoiqIoeUIVFEVRFCVPmBR0AgWpRIkS0s3NraDTUBRFKVJOnjwZIaUs+ej2F7qguLm5ceLEiYJOQ1EUpUgRQtw0tl1d8lIURVHyhCooiqIoSp5QBUVRFEXJEy/0HIqiFLSUlBSCgoJITEws6FQUJRMLCwtcXFwwNTXNVrwqKIpSgIKCgihWrBhubm4IIQo6HUVJI6UkMjKSoKAg3N3ds7WPuuSlKAUoMTERR0dHVUyUQkcIgaOj41ONnlVBUZQCpoqJUlg97Z9NVVByQEpJ0IiRRP74EzIlpaDTUZQcu3HjBtWrV892/N9//03nzp3z7POTkpLo27cvHh4eNGjQgBs3bmR735MnT+Lt7Y2HhwcffvghxtZ2unTpEo0aNcLc3Jw5c+bkWd6Kcaqg5IBMSMCQmEDYF19wvWcv4tXDkYry1HQ6HT/88AP29vb4+/szatQoJkyYkO3933vvPZYtW8bVq1e5evUqO3fuzBTj4ODAN998w9ixY/MydSULqqDkgMbKCtelS3FZ+C36+/e4+drr3JkwEV1UVEGnpig5FhAQQK1atfD19SUxMZE333wTb29vatWqxV9//ZUpftq0aQwaNIiXX34ZNzc3Nm7cyPjx4/H29qZ9+/akGBm9t2jRgkmTJtG8eXMWLFjAli1bGDRoEAC9evVi7969RkcajwoODiY2NpZGjRohhGDgwIFs3rw5U5yTkxP16tXLdJfSjRs3qFKlCkOGDKF69eoMGDCAPXv20LhxYzw9PTl+/Hh2f2xKOuourxwSQlCsTRusX3qJiKXfEbttGyKbt9YpijGf/n4BvzuxeXrMamVs+aSL1xPjLl++TL9+/fjxxx/x8fFh7ty5AJw7d45Lly7x8ssvc+XKlUz7Xbt2jb/++gs/Pz8aNWrEb7/9xpdffkmPHj34448/6N69e6Z9oqOj2b9/PwA//vgjrq6uAJiYmGBnZ0dkZCSRkZH07dvXaK5///03t2/fxsXFJW2bi4sLt2/ffvIPJB1/f3/Wr1/PsmXLqFevHmvWrOHQoUNs3bqVzz//3GiBUh5PFZRc0lhZ4TR6FCXeH4bG3ByZksKdiR/h8PprWPr4FHR6ivJE4eHhdOvWjd9++w0vr9Tic+jQIYYPHw5AlSpVKF++vNGC0qFDB0xNTfH29kav19O+fXsAvL29s5wPSV8ojI1GhBBUrlyZM2fOZJlzVvs9DXd3d7y9vQHw8vKidevWCCEem7vyeKqg5EBCsh6DlFib//fj05ibA5B86xbxvr7E/vEHxXv3ouTo0ZjY2xdUqkoRkp2RRH6ws7PD1dWVw4cPpxWU7Fx2AjB/8Odeo9Fgamqa9o+6RqNBp9MZ3cfa2jrtexcXFwIDA3FxcUGn0xETE4ODgwOXL19+7AjFxcWFoKCgtG1BQUGUKVMmWzk/mvvDfNOfS1a5K4+n5lBy4Jt9V2kzbz/bzwVn+otnXrEiFbZvx+GNN4jeuImA9h2I+vVXpMFQQNkqyuOZmZmxefNmVq1axZo1awBo1qwZP//8MwBXrlzh1q1bVK5cOc8/u2vXrqxcuRKADRs20KpVqwwjFGNfxYsXx9nZmWLFinH06FGklKxatYpu3brleX7K0ymyBUUIUUEI8YMQYsMj262FECeFEHl3b+Mj2lQthb2VGcN+PsXAFce5HhGX4X2tjTWlJk7AfeNGzDw9iFq7DlRBUQoxa2trtm3bxvz589myZQvDhg1Dr9fj7e1N3759+emnnzL8Rp9XBg8eTGRkJB4eHsybN4/Zs2dne98lS5YwZMgQPDw8qFixIh06dABg6dKlLF26FICQkBBcXFyYN28eM2bMwMXFhdjYvJ2nUv4jsju0fRaEECuAzkCYlLJ6uu3tgQWAFlgupZyd7r0NUspe6V5PB+KAC1LKbY/7vLp168qcroei0xv439GbzPvzCkk6A0ObV2BYCw8szbQZ4qSU6O/excTREX1MDBHLllFi6FC0trY5+lzl+XLx4kWqVq1a0GkoSpaM/RkVQpyUUtZ9NLawjVB+Atqn3yCE0AKLgA5ANaC/EKKasZ2FEG0APyA0f9MEE62GNxu7s3dsczrVcObbff60nb+fPX4ZP1oIgYmjIwBxR45w98efuNahI9GbN2f7OrWiKEpRUKgKipTyAHD3kc31AX8pZYCUMhlYB2R1sbQl0BB4FXhbCJHv5+dUzIL5fX1Y905DLE21DFl1giErfQm8G58p1rZ9e9zW/4qpS1mCJ37EzddfJ9HInTOKoihFUaEqKFkoCwSmex0ElBVCOAohlgK1hBAfAUgpJ0spRwJrgO+llJkmLoQQ7wghTgghToSHh+dZkg0rOLJ9RFMmdazCP9ciaTNvP9/uvUqSTp8hztLLC7e1ayn92XSSr/oT/s03eZaDoihKQSoKtw0bu7lcSikjgXeN7SCl/Cmrg0kplwHLIHUOJS8SfMhUq+GdZhXpUrMMM7ZdZO7uK/x2KohPu1WneaWSaXFCo8G+d2+KtWmDTE59mjj51i0Sz5+nWIcOqlmgoihFUlEYoQQBruleuwB3CiiXbHG2s2TRgNqseqs+QggGrTjOe6tPcic6IUOcib09pqWcALi7ejW3R48hcPBgkgKuF0TaiqIouVIUCoov4CmEcBdCmAH9gK0FnFO2NKtUkp0jmzL25Ur8dTmMNvP2s3T/NZJ1mW8hLjVhAqU+nkLCufMEdOtG2Lz5GOIzz8MoiqIUVoWqoAgh1gJHgMpCiCAhxGAppQ74ANgFXAR+lVJeKMg8n4a5iZYPWnmye1RzXqpYgtk7LtHxm4McuRaZIU5otTgMGEDFHdux69iRyGXLiPj++wLKWlFyp0WLFuT0lvyCNGvWLDw8PKhcuTK7du3K9n53796lbdu2eHp60rZtW6KyaBTbvn17ihcv/tglAN544w02bNiQafvKlSvx9PTE09Mz7WHQ7Jg3bx7VqlWjRo0atG7dmps3b6a9p9Vq8fHxwcfHh65du2b7mFmSUr6wX3Xq1JHP2h6/ENnki72y/IRt8sO1p2RoTILRuDhfX6mLjZVSSpng5yeTbt58lmkqz4ifn19Bp5AvmjdvLn19fXO0r06ny+NsniwlJUVeuHBB1qhRQyYmJsqAgABZoUKFbOcybtw4OWvWLCmllLNmzZLjx483Grdnzx65detW2alTpyyPNWjQILl+/foM2yIjI6W7u7uMjIyUd+/ele7u7vLu3bvZym3fvn0yLi5OSinl4sWLZZ8+fdLes7a2fuL+xv6MAiekkX9TC9UI5UXQumopdo9qzoetPdlxLoTWc/ez4tB1dPqMl8Gs6tZFW6wYACHTPyOgcxfCFy7CkJRUEGkrz6kJEyawePHitNfTpk1j7ty5SCkZN24c1atXx9vbm19++QVI7aPVokULevXqRZUqVRgwYMBjn6cyGAwMGjSIKVOmAKlrmNStWxcvLy8++eSTtDg3NzemT59OkyZNWL9+PS1atGDUqFE0a9aMqlWr4uvryyuvvIKnp2fasQC6d+9OnTp18PLyYtmyZWnbbWxsmDx5MjVr1qRhw4aEhmZ+NG3atGm88847vPzyywwcOJAtW7bQr18/zM3NcXd3x8PDI9tt7NO34R80aFCWnYpbt25NsQd/r5/Grl27aNu2LQ4ODtjb29O2bVuj678Y07JlS6ysrABo2LBhhh5oea0o3OX13LEw1TK6bSV61CrLJ1svMH2bH+tPBjGjuxd1yjtkii/79XzCvviCiIULidm6ldJTJmPTrFkBZK7kqx0TIeRc3h6ztDd0yLqdSb9+/Rg5ciTDhg0D4Ndff2Xnzp1s3LiRM2fO8O+//xIREUG9evVo9uDP3OnTp7lw4QJlypShcePGHD58mCZNmmQ6tk6nY8CAAVSvXp3JkycDMHPmTBwcHNDr9bRu3ZqzZ89So0YNACwsLDh06BCQ2j7FzMyMAwcOsGDBArp168bJkydxcHCgYsWKjBo1CkdHR1asWIGDgwMJCQnUq1ePnj174ujoSFxcHA0bNmTmzJmMHz+e77//PkMheujkyZMcOnQIS0tLPvjgAxo2bJj2XvqW+E2bNuXevXuZ9p8zZw5t2rQhNDQUZ2dnAJydnQkLC3vy/5uncPv27bQ2/4/m1rdvXy5fvpxpn9GjRzNw4MAM23744Ye0FjUAiYmJ1K1bFxMTEyZOnGh0uYGnoQpKAXIvYc3KN+ux83wI07f50XPJEXrXcWFihyo42vzXN8m0VCnKzptH8d69CZn+GYHvDMVl0UKKtW5dgNkrz4NatWoRFhbGnTt3CA8Px97ennLlyjF//nz69++PVqulVKlSNG/eHF9fX2xtbalfv37aWiQ+Pj7cuHHDaEEZOnQoffr0SSsmkFqwli1bhk6nIzg4GD8/v7SC8mh34YfX9L29vfHy8kr7B7tChQoEBgbi6OjIN998w6ZNmwAIDAzk6tWrODo6YmZmljZPUadOHXbv3m30/Lt27YqlpSXw+Jb4Bw8ezOZPNH88LreHo8cnWb16NSdOnEhbiwbg1q1blClThoCAAFq1aoW3tzcVK1bMcZ6qoBQwIQQdvJ1pVqkk3+y7yg8Hr/OnXyjj2lWmf/1yaDX/PZNi3agRFbZsJnrjJmyaNwcg8fIVzN3dEGZmBXQGSp55zEgiP/Xq1YsNGzYQEhJCv379gMe3r0/fJFKr1WbZ6v2ll17ir7/+YsyYMVhYWHD9+nXmzJmDr68v9vb2vPHGGyQmJqbFp29rn/5z0reWf/hap9Px999/s2fPHo4cOYKVlRUtWrRIO176VvqPy9FYK/2H0rfEf9IIpVSpUgQHB+Ps7ExwcDBOTk5GPy+nXFxc+PvvvzPk1qJFCyB7I5Q9e/Ywc+ZM9u/fn+Fn+fD8KlSoQIsWLTh9+nSuCoqaQykkrM1N+KhDVXaMaEpV52JM2XyeHosPczYoOkOcMDPDvl9fhIkJhvh4bg0eTED3HsQdPVpAmStFXb9+/Vi3bh0bNmygV6/UPqvNmjXjl19+Qa/XEx4ezoEDB6hfv/5THXfw4MF07NiR3r17o9PpiI2NxdraGjs7O0JDQ9mxY0eu8o6JicHe3h4rKysuXbrE0Vz+HejatSvr1q0jKSmJ69evc/Xq1bRzPnjwoNFW+m3atEnb9+GdVytXrszzVvrt2rXjzz//JCoqiqioKP7880/atWsHpI5QjOX2sJicPn2aoUOHsnXr1gyFLioqiqQHc7IREREcPnyYatWMtknMNlVQChnPUsVY+3ZDFvTzITgmkW6LDjN50zmi45MzxWqsrCgzcwYyJYVbb7zJ7dFjSAnN22u3yvPPy8uLe/fuUbZs2bTLSj169KBGjRrUrFmTVq1a8eWXX1K6dOmnPvbo0aOpXbs2r7/+etr69F5eXrz11ls0btw4V3m3b98enU5HjRo1+PjjjzPMf+SEl5cXffr0oVq1arRv355Fixah1WqfvCMwceJEdu/ejaenJ7t372bixIkAnDhxgiFDhqTFNW3alN69e7N3715cXFyyvDV56NChuLi44OLiQqNGjXBwcODjjz+mXr161KtXj6lTp+LgkHm+1Zhx48Zx//59evfuneH24IsXL1K3bl1q1qxJy5YtmThxYq4LSqFqX/+s5aZ9/bNwLzGF+buvsvLIDewsTZnYoQq9arug0WRszWJITCTy++VEfv89wsQE961bMEu33rZSeKn29UphV5Tb1yvpFLMwZWqXavz+QRMqlLBm/Iaz9P7uCH53Mi4QpLGwoOTwD6iw7XcchgzGtGxZADVaURTlmVIFpQioVsaWX4c24qteNbgREUfnbw/y6e8XiE1MyRBnVq4cJYcNQwhBclAQ19q1486EiegiIgooc0VRXiSqoBQRGo2gd11X9o1pwasNyvHTPzdoPXc/m0/fNnpHjomDAw4DBxKzfTvXOnTk7po1SL3eyJEVRVHyhiooRYydlSkzunuz5f3GlLGzYOQvZ+j//VGuhma8pVFjZYXT6FFU2LIZi+pehE7/jBv9+iOTM0/uK4qi5AVVUIqoGi7F2TisMTN7VOdi8D06LDjIrO0XiUvKeL+9eYUKlFuxgrLz5mLTvHna8yqGhARjh1UURckxVVCKMK1GMKBBefaNac4rtcvy3YEA2szbz45zwRkugwkhsO3YkZIfvA9A/KlT+LdqTdT69UhD5lb6iqIoOaEKynPA0cacL3vV5Lf3GlHcyoz3fj7FoB99uR4RZzRea2uLWcUKhHw8lZv9XyXRz+8ZZ6w871T7euPt67PTLl61ry+iXwXRvj6/pej0csWhAFl96k7pOWm7nLvrkkxIztyC22AwyKhNm+TllxpLv6rVZOicuQWQraLa12f2PLevz067eNW+Xik0TLQa3mzszt4xzenoXZpv9vnTZt5+9vhlbN8thKB49+5U3LEd+379MCnhCPz3C4byYlDt659t+/qcUu3r85kQogIwGbCTUvZ6sM0aWAwkA39LKX8uwBQLlJOtBV/3q0XfeuWYuuU8Q1adoE1VJz7p4oWrg1VanNbWltJTP057HbvtD6J/+YVSUz/GolKlgkj9hfXF8S+4dPdSnh6zikMVJtSfkOX7qn39s21fn9N28ap9fQ4IIVYAnYEwKWX1dNvbAwsALbBcSjlbShkADBZCpL/Y+AqwQUr5uxDiF+CFLSgPNaroyPYRTVlx6DoL9l6lzbz9fNDSg3eaV8DcxEifIo0g6epVrvd4BYeBAynx/vtobawzxynPBdW+/tm2r89pu/jH5aba12ftJ2AhsOrhBiGEFlgEtAWCAF8hxFYppbGZZBfg4QpF6im+B0y1GoY2r0hXnzJ8ts2PubuvsPH0bT7t6kWzSiUzxNp16oT1Sy8RPm8ed3/8kdjt2yk97ROKtWxZQNm/OB43kshPqn19qmfRvj6n7eJV+/ockFIeAO4+srk+4C+lDJBSJgPrgKx6QweRWlQgi3MTQrwjhDghhDgRHh6eF2kXGc52liweUIdVb6W25B644jjDfj5JcEzGZ1JM7O1x/uwz3NatRevogExJMXY45Tmh2tenyu/29blpF19U2tcXthGKMWWBwHSvg4AGQghHYCZQSwjxkZRyFrARWCiE6AT8buxgUsplwDJI7Tacr5kXUs0qlWTnyKZ8fyCAb/f58/flcEa09uStJu6Yav+rw5Y+PrivXw+a1G2RP/yAPvYeJd4diubBZQKl6Muqff2RI0eoWbMmQoi09vWXLj3dHM/o0aOJiYnh9ddf5+eff05rX1+hQoU8aV+/dOlSatSoQeXKlfO0fb2JiclTt6/v06cPP/zwA+XKlWP9+vVAavv6pUuXsnz5ci5evMjQoUPRaDQYDIbHtosfOnQoI0eOBMDV1ZUjR46kta8Hcty+HqBcuXJs3br1qfLJrkLXvl4I4QZseziHIoToDbSTUg558Pp1oL6UcnhuP6uwt69/FgLvxvPp737suRiKp5MN07tVp1FFR6OxwZ9MI/qXXzAtU4ZSkydh06pV2mUFJWdU+3qlsHve2tcHAa7pXrsAdwool+eOq4MVywfVZfnAuiSk6On//VFGrjtN2L3ETLHOn06j/P9WobG2Juj9Dwh69z2Sg24XQNaKohRGRaGg+AKeQgh3IYQZ0A/YWsA5PXfaVCvF7lHNGd7Kg+3nQmg9Zz8rDl1Hp8/YmsWqXj3cN/6G0/jxxJ88iT6LJ4IVRXnxFKqCIoRYCxwBKgshgoQQg6WUOuADYBdwEfhVSnmhIPN8XlmaaRnzcmV2jWqGT7niTN/mR5eFhzl5M2PREKamOL71Jh5//42ld+rd3eELF3H/wIGCSFtRlEKiUBUUKWV/KaWzlNJUSukipfzhwfbtUspKUsqKUsqZBZ3n8869hDWr3qrPkgG1iY5PpueSfxi/4V8i7ydliHv4fIohIYHY7dsJfGcoQcM/JOWOuiKpKC+iQlVQlMJDCEEHb2f2jG7O0OYV2HjqNq3m7ufnYzfRGzLeyKGxtMR98yZKjhrF/YMHudapMxHLvldrryjKC0YVFOWxrM1N+KhDVXaMaEpV52JM3nSeVxYf5mxQdIY4jZkZJYa+Q8U/tmHd+CUili5Fd/fRR4oURXmeqYKiZItnqWKsfbshC/r5cCcmkW6LDjNl8zli4jM+9GhatiyuCxdSYetWTEuXRkpJxHfLSMmit5FSsG7cuEH16tWNvhccHIyHhwe1a9c2+pT4Q5GRkbRs2RIbGxs++OCDDO9NnjwZV1dXbGxsjB7/5Zdf5ubNm9SpUwcfHx+8vLxYunRp7k5KKTCqoCjZJoSgm09Z9o5pzhsvubHm2C1azv2bX08EYnjkMpiZS1kAkv39iVi0iIAOHbm7ciUyixYYSuFy7949unfvzhdffMGgQYPo1asXKVl0TLCwsOCzzz5jzpw5md7r0qVLlh17d+7cSbt27XB2duaff/7hzJkzHDt2jNmzZ3NHzcMVSaqgKE/N1sKUT7p4sW14U9xLWDN+w1n6fHcEvzuxmWLNPT2p8PtWLOvUJnTWbK6/0pP4kycLIGvlSQICAqhVqxa+vr7079+fCRMm0LNnT0aMGEHXrl15++23je5nbW1NkyZNsLCwyPRew4YN056+f9TOnTvp0KEDZmZmaf2lkpKSMKhVRIusotB6RSmkqpWxZf3QRmw4FcTsHZfosvAQAxuVZ3TbShSzME2LMytfHtfvvuPenj2Efj6L22PG4vHnrrT17ZX/3Hx9YKZtxTq0x+HVVzEkJBD4ztBM79v16EHxV3qgi4ri9ocjMrxX/n+rMsUbc/nyZfr168ePP/6Ij48P27Zty/D++++//xRn8WR6vZ7Lly+ntfoIDAykU6dO+Pv789VXX6U1LVSKFjVCUXJFoxH0qevKvjHN6VfPlZ/+uUGrufvZcuZ25nXt27al4h/bcF26BGFmhkxOJnrTZqReNYYuSOHh4XTr1o3Vq1fj4+PzTD7z2LFjNGjQIO21q6srZ8+exd/fn5UrVxpdEEsp/NQIRckTxa3MmNnDmz51Xfl4y3lGrDvDuuOBTO/mhWepYmlxGisrLKpUASB2506CP/qIqNWrKf3JVCwfrIvxInvciEJjafnY903s7bM9IknPzs4OV1dXDh8+jJeX1xPjN23axKeffgrA8uXLqVs3U0unJ9qxYwft27fPtL1MmTJ4eXlx8ODBtM7HStGhRihKnqrpWpxNwxozo3t1/IJj6bDgILN2XCQuKfNkvG2XLpSZOwddWBg3+vYjeOon6FQrl2fOzMyMzZs3s2rVKtasWfPE+B49eqS1SM9JMQHYu3cvrVu3BlLX9khISF1CISoqisOHD1O5cuUcHVcpWKqgKHlOqxG81rA8+8Y0p0etsny3P4A28/az41xwpstgdp06UWHHdhwGDiT6t9+4M3ZcAWb+4rK2tmbbtm3Mnz+fLVu2PPX+bm5ujB49mp9++gkXFxf8/FLXvxs/fjwuLi7Ex8fj4uLCtGnTCA8Px8LCAltbWyC1m22DBg2oWbMmzZs3Z+zYsXh7e+fp+SnPRqFrX/8sqfb1z8aJG3eZsvk8l0Lu0axSST7t6oV7iczLCidevgLSgEWVKuiiotAFB2ORy/UZCrsXsX396tWrCQoKYuLEiQWdipINT9O+Xs2hKPmurpsD24Y3YdWRm8zbfYV28w/wbvMKDGvpgYXpfwsYWVSulPZ95NKl3P3fauz796fkiA/RPvhtVin6XnvttYJOQckn6pKX8kyYaDW81cSdfWOa08G7NN/s86ft/P3svWj8bp4Sw4Zh368vUWvWcK1jJ2K2bHnsOueKohQ8VVCUZ8rJ1oIF/Wqx5u0GmJtoGbzyBENWniDwbnyGOK2dHaWnTsVt/XpMy5ThzoSJRHz7bQFlrShKdqiCohSIlyqWYPuHTfmoQxX+uRZB2/n7WbjvKkm6jM+kWFb3wm3dWkp/+inFe/YEICU0FP39uIJIO1+okZdSWD3tn83npqAIITRCiJlCiG+FEIMKOh/lycxMNAxtXpE9o5vTqooTc/68QoevD3LwaniGOKHRYN+3D6ZlU/uDBU+aTECnTsTu2FHk/zG2sLAgMjKyyJ+H8vyRUhIZGWm0pU5WCvVdXkKIFUBnIExKWT3d9vbAAkALLJdSzhZC9AC6AXeBP6SUe590fHWXV+Gy/0o407Ze4HpEHJ28nZnSuSrOdpaZ4hLOnCF4+nSS/C5i/dJLlJoyBfMK7gWQce6lpKQQFBREYmJiQaeiKJlYWFjg4uKCqalphu1Z3eVV2AtKM+A+sOphQRFCaIErQFsgiNQ15/sDXYEoKeV3QogNUsonPmarCkrhk6TTs2x/AAv/8kerEYxo7clbTdwx1WYcTEu9nqi16wj/+msMSUm4LvwWm+bNCyhrRXmxZFVQCvUlLynlAVJHHOnVB/yllAFSymRgHakjkyDg4WPWWTaHEkK8I4Q4IYQ4ER4enlWYUkDMTbQMb+3JntHNeamiI7N2XKLjgoMcDYjMECe0WhxeG0DFHdux79MHyzp1ANBHR6vLR4pSQAp1QclCWSAw3eugB9s2Au2EEN8CB7LaWUq5TEpZV0pZt2TJkvmbqZJjrg5WLB9Uj+UD65KQoqffsqOMXHeasHsZLw2ZlCxJ6Y+noLWxQSYnc2PAawS9+x7JgYFZHFlRlPxSFAuKMLJNSinjpZSDpZTDpZSLnnlWSr5oU60Uu0c1Z3grD7afC6H1nP38ePg6Or2RNTOEoHjPnsT7+hLQuQvhixZhSEp69kkryguqKBaUIMA13WsXQC3v9p1ZheUAACAASURBVByzNNMy5uXK7BzZFJ9yxfn0dz+6LjzMyZsZG0kKU1Mc33qTCtv/wKZVSyK+XUhAl64kB90uoMwV5cVSqCflAYQQbsC2dJPyJqROyrcGbpM6Kf+qlPLC0x47x5PyyfGgf/ibb7oBkxAZt4n0g6lHtz1mv7w4ljA2kCv6pJTsOB/C9N/9CIlNpE9dFyZ2qIqDdebFuu4fPkz0+g2UnfMVwsQEmZKCeORuFUVRnl5RvctrLdACKAGEAp9IKX8QQnQEvib1tuEVUsqZOTl+jgvKzklwtCheVcvHoqYxBa0paM3AxCz1vxm+TMHEHMxswMIWzG3Bwu6//9o4QTFnsHUGi+JPLIhxSTq+2XuVHw5dx9rchPHtK9O/Xjk0GuP76aOjud6zF8X79cVx0CC1WqSi5EKRLCj5LccF5cZhCDkLGX52D75P25buvUe3PXa/Jx0rp/vlZQ5GjmXQgz75ka8U0CWl+z4Rku9DYgwkxoIuAaNMLFMLi6MHlKgEJTyhRGVwrgFmGbsUXwm9x8ebz3Ps+l1qutgxo7s33i52mQ6ZEhZGyKfTub93L2YVK1L644+xbtggU5yiKE+mCooR6jmUAqZPSS0sidFwPwzu3YF7IRB7B2JvQ6Q/RPj/V3iEBpy8wKUOuDaEiq2gWCmklGw5c4cZf1wkMi6JAQ3KMe7lKthZZb68de/vvwmdMZOUoCBsO3XCedbnaNRoRVGeiiooRqiCUgQYDBATCGEX4fZJCPKF26cgKSb1/dLe4NEGKncitkRN5u2+yqojN7C3MmNihyr0rO2S6TKYITGRyO+XkxRwDZf584HUuRnxnM47KUpeUwXFCFVQiiiDAULPgf8e8N8LgcfAoIPi5aF6T/xLd2D8gRRO3Yqmbnl7PutenarOmddTeVhEkm/c4Pb4CZSaOAGr2rUL4IQUpWhRBcUIVVCeEwnRcOkPOL8BAvaD1CPL1uN4ie6MOOdOeKJgUCM3RrX1pJhF5stg8SdOcHvceHTBwdj16IHT2DGYODoWwIkoStGgCooRqqA8h+6Hw7n1cOIHiPTHYGHPoWLtmRD0EnqbMkzuVJWuNctkurxliI8nYskSIn/8CY2VFU5jxmDft08BnYSiFG6qoBiR04Ky9tJa/rzxJ0CGf5jEg1tpRdottWT5nrH9njY+x/sZmSvI8ec8iDcRJphqTTHVmGKmMfvve60ZphpTTDQmmGnNsDaxxsbMhmJmxbAxTf1vMbNiWJlY5e0chpRw/QD4LodLf2AQGvaYtmB2bDtKuXvzWXcvPJyKZdot6do1Qj6bgXklT0pPmpR3+SjKc0QVFCNyW1BkultpH/05Gnvv4TZp5Dbg7Mbn9HMy7Zfuda4/R0r0Uk+KIYVkfTIphhRS9CnopI7ssjSxpKRlSZysnHCycqKUVSnK25bH3c4ddzt37C3ss32sTKJuwpGFyFOrQJfEHuqzIKUHjZu04MNWnlibm2Q8VykhJQVhZka8ry8x2/7AadRItMWL5zwHRXmOqIJihLrklb8M0oDOoEsrMkn6JOJT4rmXco/7yff/+2/yPSISIgiLD8vwlWxITjuWvbk91UpUo2aJmniX9KZGyRrYmmWeaH+s++FwbCmG48vQJMXyu74hP1sMYFDXtrSvXtroCCnyp58I+2oOWltbnMaOwa5HD4SmKHYsUpS8owqKEaqgFF56g57guGCux1wnICaAa9HXOBdxjmvR15BINEKDdwlvGpdtTJMyTfAq4YVGZPMf+oRo+Odb9EcWgy6R33RNOVrubT58pRVuJawzhSdevkzI9M9IOHkSSx8fSn8yFYuqVfP4jBWl6FAFxQhVUIqe+8n3uRB5Ad8QXw7fPsyFyAtIJE6WTrR3b09H945Uc6yWvfmY++EYDs7D4LscqdezRrYlruEY3mpbBwtTbYZQKSUxm7cQ9tVXlHh3KA4DB+bTGSpK4ZergiKEcMjGZxiklNE5Sa6gqIJS9N1NvMvh24f58+afHLp9CJ1BR3nb8rzi+Qo9PHpkb+4l5jYJez7H/Nwa7kkL/mfWB6/uY2np5ZopVB8Tg8baGmFiQuzOXciUZGw7d1YPRSovlNwWlERSW8Q/7m+NVkpZLucpPnuqoDxfYpJi2HtrL1uvbeVk6ElMNaa87PYy/Sr3o2bJmk/+Rz/Uj6gtE7C/c4CbBie2l36Pzn2H4uqY+TIYQOB7w7j/119Y1a9P6akfY+7hkQ9npSiFT24LymkpZa3cxhQ2qqA8v/yj/Pn1yq9svbaVuJQ4fEr6MMR7CE1dmj5xriXl8m7ubZ2AQ9w1TsjKBNT6iG6du2Ju8shlMIOB6PUbCJ83D31cHA6DBlJy2DA01sYLkKI8L3JbUCyklIm5jSlsVEF5/sWnxLPZfzMrL6zkTtwdPIp7MNh7MO3d2mOiMcl6R72O6H9WoPn7c2z1UewyaYld189pWKNaplBdVBRhc+cSs+E3XJYspljLlvl4RopS8HJcUIQQbYE+wCIp5RkhxDtSymX5lOczpQrKiyPFkMLO6ztZcX4F/tH+uNm6McxnGO3c2j1+xJJ0j1tbPqO03w8kSVN2lxxEo1cn4eyQuUV+kr9/2mWvmG1/YOFVDXN39/w6JUUpMLkpKJuAN4EpwHagl5RyWL5kmQtCiO5AJ8CJ1OL355P2UQXlxWOQBv4K/IuFpxfiH+1PJftKDK81nOYuzR87x5IUeoU7v4zG/e5BAmQZLvlMom3XAZhqMxcjQ2Ii/m3aYoiJwWHwW5QYOhSNpWV+npaiPFNZFZTs3LgfLqWMllKOBV4G6uV5dlkQQqwQQoQJIc4/sr29EOKyEMJfCDERQEq5WUr5NvAG0PdZ5agULRqhoXW51mzosoHZTWeTqEtk+L7hvLb9NXxDfLPcz7xUJdw/3EZYl9VYmwo6/vsBvrPac/rf05k/w8KCCps3YduxA5FLvyOgcxfu7duXn6elKIVCdkYo3aSUW9K9Hi6l/DbfM0v9rGbAfWBVujXltaSuKd8WCCJ1Tfn+Ukq/B+/PBX6WUp560vHVCEVJMaSw1X8rS/5dQmh8KC1dWzK6zmjc7Nyy3kmXxNUtX1L23EK0Us/+Ev3xGTANJ4fMHYrjjh8n9LPPSLoWQMUd2zErXz7/TkZRnpEi+2CjEMIN2JauoDQCpkkp2z14/dGD0NkPvnZLKfc85njvAO8AlCtXrs7NmzfzL3mlyEjUJbL64mqWn1tOki6JvlX68m6NdylukXX/roTIQK6vHUu1iJ0ES0cu1xhPk25vY/Lo3WApKcSfOIF1o0YA3Nu7F+smTdCYm+frOSlKfsmTgiKEqAtMBsoDJqQ+lyKllDXyKlEjn+lGxoLSC2gvpRzy4PXrQANSRy2DSB2xnJFSLn3SsdUIRXlUREIEi88s5rerv2Ftas3QGkPpX6U/Ztqslwm+/e8+UraNxS3lGv+aeGPSZR5eNesbjU0KCCCgYydMy5ej9JSPsWnaJL9ORVHyTV4VlMvAOOAcYHi4XUqZb7/mGykovYF2jxSU+lLK4U97bFVQlKz4R/kz9+RcDt0+hIuNC2PrjaWVa6ssJ+6lXsf537+l3Jk5WMkE/inZG+8Bs3Cwz9xk4v7hw4R+NoPkGzco9vLLlPpoIqbOzvl9SoqSZ3IzKZ9euJRyq5TyupTy5sOvPMoxu4KA9D0xXEh9il9R8oyHvQdL2izhuzbfYWFiwci/RvLunncJiAkwGi+0Jnh3H4V2xCkuOHWkecQ6UhbU4dDmZRj0hgyxNo0b4751CyVHjuT+gQNc79MHQ3Ky0eMqSlHytCOU1kB/YC+Q9HC7lHJj3qeW9pluZByhmJB6eas1cJvUS1yvSikvPO2x1QhFyY4UQwq/Xv6VRacXkaBLYEDVAbxb811szGyy3OfWv3+j/3007rpr/Gvqg2W3eVSqXidTXHLQbZKuXKZYq1ZIKUk8fwFL7+r5eTqKkmt5dclrNVAFuMB/l7yklPKtPMky8+etBVoAJYBQ4BMp5Q9CiI7A14AWWCGlnJmT46uCojyNyIRIvj39LRuvbsTBwoFRdUbRpWKXLB+MlHodZ7d8jfvZeVjIRI6W6k/NV2dgV9x4w8rYHTu4PWo0tp064TRhPKZOTvl5OoqSY3lVUM5JKb3zNLMCpAqKkhPnI84z6/gszoafpUbJGkyqPwmvEl5Zxt+LvIP/z2OodXc7IThyvc5kGnZ6M9NCXYbERCKXfU/k998jzMwo+eFw7AcMQJg8pkWMohSAvCoo3wPzHz7zUdSpgqLklEEa+P3a78w/OZ+7iXfp4dmDD2t9iKNl5mdRHgo4tQexfRzuugDOmNWmWI/5VKzqkyku+eZNQmbMJO7gQWxatMB16ZL8PBVFeWp5VVAuAhWB66TOoeT7bcP5SRUUJbfuJ9/nu7PfsdpvNZYmlrxf6336Vu6bZeNJgy6FMxvn4Om3AHOZzDHn1/B5dTrFbDM+7yKl5N7u3WgsLbFp2hRDcjKG+/cxccjO0kSKkr/yqqAYfcy3AO70yhOqoCh5JSAmgC+Of8E/d/6hkn0lpjScQi2nrFdziA4NJGDdGGpH7SKYEgTWn0q99q9nuV59xJIlRP60EqdRIyneuzdCqzUapyjPQpF9Uj4/qYKi5CUpJXtu7eFL3y8JiQuha8WujK4z+rGXwa4e34XJrvG4629wxrwe9r3mUd4z84A/6do1QqZ/RvyxY1hUr07pT6Zi6f3cTGcqRUxu10M5JaWsnduYwkYVFCU/xKfEs+zsMlb6rcRSa8nw2sPpU6kPWo3xUYVel8KpDV9Q9dJCzGQKvi4DqfXqp1hZ22aIk1IS+8d2Qr+YjT4iEqfx43F8841ncEaKklFuC0oCcPVxIYCdWgJYUf4TEBPArGOzOBp8lKoOVZnUYBI+Tpkn4R+KDLnFjbWjqROzmzs4EfzSNGq36Z/pMpj+/n0ivv0W2y5dsazuhSEhAWFunuXlMkXJa7ktKNlpkaqXUgblJLmCogqKkt+klPx580++9P2SsPgwenj0YGSdkThYZD25funodiz+nICb4RZnLBpQovd8XCpmfVvyncmTSQ64TumpH2NRtWp+nIaiZKDmUIxQBUV5VuJT4ln671L+5/c/rEytGFF7BD09e2Z5GUyXnMTJ9bOpfmUxJug5Ve4NavWfhoVV5qfzozdtJuyrr9BHR2M/YAAlPxyOtlix/D4l5QWmCooRqqAoz9q16Gt8fuxzjoccp5pjNaY0mIJ3yawn18NvX+fWutHUubeP26IUEU0/o2arzOvH6WNiCPv6a6LX/YK2hCMu8+ZhVe+ZrYWnvGBUQTFCFRSlIEgp2XljJ1/5fkVEQgSveL7CiNojsLcw3pIF4MLBLdj8NYnyhiBOWzWmdN/5OJevnCku4dx5wubMocwXszEtXRppMKi5FSXP5VlBEUKUBJBShudRbgVGFRSlIN1Pvs+Sf5fw88WfsTGzSbsMllVvsOSkRE6tm0GNgGUIJGcqvE3tfh9jbm58vXopJYFD3sa8cmVKvj8MjbV1fp6O8gLJVft6kWqaECICuARcEUKECyGm5nWiivKisDGzYVy9cazvsh6P4h5MPzKdAX8M4EKE8cbZZuYWNBw0g9jBh7lkU59G1xcROrsO5w5sMRovk5MxKV2KuytWcK1TZ2J37uJFviKh5L/s3uU1CugIvCOlvP5gWwVgCbBTSjk/X7PMJ2qEohQWUkr+uP4Hc3zncDfxLr0r9ebD2h9iZ26X5T5n/1qPw4EpuMgQTtq0wLX/1ziVdc8UF3/qNCHTp5N06RLWjRvj/PnnmJZSnYyVnMvtbcOngbZSyohHtpcE/pRSZt1johBTBUUpbO4l32PxmcWsubQGWzNbRtUZRXeP7lleBktMiOP02k+pfXMFOrSc83yPOn0+wtQs43r1Uqcjas1aotevx23dWnX5S8mV3BaU8w8XuHqa9wo7VVCUwury3cvMPDaT02GnqVGiBpMbTqaaY7Us4+8E+BG+fgQ1E45zXVOehLZfUK1Rh0xxDyfpDcnJBA17H/tX+1OsVav8PBXlOZTbJYAftz5poVm7VAjRXQjxvRBiixDi5YLOR1FyqrJDZVa2X8nMJjMJuh9E/z/6M+PoDGKSYozGl6lQjRrjdnH6pUVYGuKptqsfvvN7Ex5yK0Pcwzu+dGFhpIQEEzTsfQLffY/koCL1TLJSSGV3hKIH4tJvSve9hZTSNNeJCLEC6AyEpR/xCCHaAwtIXZ1xuZRydjaOZQ/MkVIOflycGqEoRUFsciwLTy/kl8u/UNy8OKPrjKZrxa4IIYzGJ9yP5d+1U6gdtJpEzPGr8iF1e43FxDTjX1OZksLdVasIX7QY9HpKvDsUxyFDEKa5/uusPOcK/XMoQohmwH1gVbr147Wkrh/fFggidf34/qQWl1mPHOItKWXYg/3mAj9LKU897jNVQVGKEr9IP2YencnZiLPUdqrN5IaTqWRfKcv4wCtniPltJNWTTuOvrYiuw1dUqds6U1xKSAihs2aTcucObuvWqtb4yhPldg6lHhAopQx58Hog0BO4AXwqpbybR0m6AdvSFZRGwDQpZbsHrz8CkFI+Wkwe7i+A2cBuKeWeJ32eKihKUWOQBjZd3cTXp77mXvI9BlQdwDCfYVibGp9klwYDp3f+iMvxGThxl+P2nfF8dQ72JZ0zHzsuDo21NbqoKMK+/IqSHw7H1DlznKLkdg7lOx7MlTwYScwGVgGxwLK8StKIskBgutdBD7ZlZTjQBuglhHjXWIAQ4h0hxAkhxInw8CL/bKbygtEIDT0r9eT37r/T3aM7//P7H103dWXH9R1GnzERGg21Ow7GavQpjpZ+ldp3tyMW1eXYhnkY9PqMx35w51fi2bPEbt/OtU6diVy+HJlcaKZJlUIuuyOUf6WUNR98vwgIl1JOe/D6jJQy657cT5NM5hFKb6CdlHLIg9evA/WllMPz4vPUCEUp6s6Gn2XG0RlcvHuRBs4NmNRgEhXsKmQZf8PPl/jNI6mWfJ4rJpUQnefh6dM0U1xy0G1CP/+c+/v2YVaxIqWnTsW6Qf38PBWlCMntCEUrhHi4SHZrYF+694wvnp03ggDXdK9dgDv5+HmKUqTUKFmDtZ3WMrnBZPwi/Oi5tSdfn/ya+JR4o/Fu1epRdeJBTtSehaMulIqbunB04ZvE3M3wiBlmLmVxXbwIlyWLkYmJRK1e/SxORynisjtCmUzqk/IRQDmgtpRSCiE8gJVSysZ5kkzmEYoJqZPyrYHbpE7KvyqlNN6b4impEYryPIlMiGTeyXlsvbYVZ2tnJtSbQKtyrbK8GywmKoJLayZQN+w3ooUt13wmUK/re5maSRoSEzEkJGBib09SQABxhw5h/+qrCJP8/F1SKcxyfZeXEKIh4Ezqk/FxD7ZVAmyedDdVNo+/FmgBlABCgU+klD8IIToCX5N6Z9cKKeXM3H7WQ6qgKM+jk6EnmXF0Bv7R/jQp24RJ9SfhauuaZbz/v4fQ/z6ayrrL+JlWx7L7PNy9GhiNDVuwgMglSzGvUoXSU6diVbtINslQcim3d3lZAO8CHsA54AcppS7Ps3zGVEFRnlcphhTWXFzD4jOL0Rl0DPYezFvV38LCxMJovEGv5+Tmb/A4N5diMg7f0n2p/urnFLPLuLKklJJ7f+4mdNYsdCEh2L3yCk5jx2DikPUKlMrzJ7cF5RcgBTgIdABuSilH5HmWz5gqKMrzLjQulDkn5rDzxk5cbFz4qMFHNHNplmV8dEQIV9aMpW7kNiKEPbfqTaJOh8GZL4PFxRGxdCmRP/6E41tv4TR6VH6filKI5LagnJNSej/43gQ4LqWsnfdpPluqoCgviqPBR5l5dCY3Ym/QyrUVE+pPoIxNmSzjr5z8C832MXjor3He3IdiryygfOXMN3MmXbuGSanSaG2sSThzBrQmWHoXydZ+ylPI7V1eKQ+/eR4udSnKi6ahc0M2dt3IiNojOBJ8hG6bu7H83HJS9ClG4yvVaYn7R8c5VvUjyiVdwXlNK44sG078/Yy9xMwrVkRrk/r8Stj8r7nRpw/B06ahj47O93NSCp+c9PISgCUQ/+B7KaW0zbcM85EaoSgvouD7wXzh+wV7b+3FzdaNSQ0m0ahMoyzjI0MDCVgzhnoxuwihJHcafUKttgMyXQbT37tH+LffErX6Z7R2djiNHYtdj+5qCeLnUKHv5VUQVEFRXmQHgw4y6/gsAu8F0s6tHePqjqOUdaks4y8e24nFrvG4G27yr0U9HHsvwKWiV6a4xEuXCPl0OgmnT1N23lxsO3bMz9NQCoAqKEaogqK86JL0Saw4t4Ll55ZjojFhmM8wXq36KqYa4x2HdclJnFg/G+8rizFBz6lyb1Cr/zQsrGwyxEmDgXu791CsTWuEVkv86dOYe3igLVbsWZyWks9UQTFCFRRFSRUYG8is47M4ePsgHsU9mNJwCnVK1ckyPvz2dW6tG0Wde39xW5QioukMarbqYzTWkJSEf6vWoBGUGj8e286ds3zYUikacjspryjKc8zV1pVFrRfxdcuviUuJ442dbzD50GQiEiKMxpcs606dMZu50HoVemFCzQNvc/rLjgTfvJwpVmNujuvSJZiWKs2dceO5NegNkvz98/uUlAKgRihqhKIoGcSnxPP9ue/56cJPWGotGV57OH0q9UGrMb5OSnJSIifXfUbNgO8RSM5UeJs6/aZiZp7xIUqp1xO9fgNh8+djiIujwtatmFdwfxanpOQxdcnLCFVQFCVrATEBfH7sc44FH6OqQ1UmN5xMzZI1s4wPuXWV4F9GUivuELc0ZYlpMQvvZt0yxenu3iV2+w4cXhsApD7LYlahgroMVoSogmKEKiiK8nhSSnbd2MVXvl8RlhBGT8+ejKw9kuIWxbPc59+/1uN4YAouMoQTxVpSrt98nMoaH4kk37xJQOcuWDVoQOkpkzFzc8unM1HykiooRqiCoijZE5cSx5IzS1h9cTU2ZjaMrD2SVzxfQSOMT8MmJsRxZu00at38ER1aznkOo06fiZiamWeIkzodUWvWEv7NN8ikJBzfHoLjO++gsTDec0wpHFRBMUIVFEV5OleirjDz6ExOhZ2iRokaTG44mWqO1bKMvxPgR/j6EdRMOM51TXkS2n5BtUYdMsWlhIUR9tUcYn//HbPy5XHfugWNubmRIyqFgSooRqiCoihPT0rJ7wG/M/fEXKKToulTqQ/Daw/H1sx4wwxpMHBmzxqc/5lGacLxtWuHe/+5lCiduaV+3LHjJPr54fjmGwDooqIwsbfPz9NRckAVFCNUQVGUnItNjuXbU9/y65VfKW5enDF1x9ClQpcsJ9cT7sfy79op1A5aTaIwx6/KCOr1Gos2i4W64o4eI/Dddykx9B0cBg9GY2aWn6ejPIXnvqAIITTAZ4AtcEJKufJJ+6iCoii55xfpx8yjMzkbcZbaTrWZ3HAylewrZRl/68oZYn8bSfWk0/hrK6Lr8BVV6rbOFJcSEkLorNnc27ULs/LlKfXxx9g0yZPFYZVcKtQPNgohVgghwoQQ5x/Z3l4IcVkI4S+EmPiEw3QDypLaGTkov3JVFCWjao7V+F/H/zGt0TSuxVyjz+99+Mr3K+JS4ozGl6vkg9eEfZysPxdbfRRVtr3C8QUDiIoIyRBnWro0Lgu+xnX5cgAChwwh+OOP8/18lJwrFCMUIUQz4D6wKt168lpS15NvS2qB8AX6k7oU8KxHDvHWg68oKeV3QogNUspeT/pcNUJRlLwVnRjN16e+5rerv+Fk6cTYemNp79Y+y8tg92OjOP/zR9QN+YX7woor1cdQt8cINNqMD1EakpO5u2IF2uL22PfrizQYQK9HmBrvOabkr0J/yUsI4QZsS1dQGgHTpJTtHrz+CEBK+Wgxebj/a0CylPJXIcQvUsq+T/pMVVAUJX+cDT/LjKMzuHj3Ig2cGzCpwSQq2FXIMv76heMkbB5FtZTzXDapjLbLPDxqNskyPnrjJiJX/EDpqVOxrl8/P05BeYxCfckrC2WBwHSvgx5sy8pGoJ0Q4lvgQFZBQoh3hBAnhBAnwsPD8yZTRVEyqFGyBms7rWVSg0n4RfjRc2tPFpxaQHxKvNF4d6/6VP3oICdqfU4JXQjuGztzbOGbxEQZ7yVmUrIkMiGRWwMHcXvceHTq73KhUJhHKL2BdlLKIQ9evw7Ul1IOz6vPVCMURcl/EQkRzD85n63XtuJs7cyE+hNo5doqy8tgMVERXPp5HHXDNxEtbAnwmUDdru9lXtc+IYGIZcu4u/wHhLk5pT+egl23zK1elLxXFEcoQUD6G9VdgDsFlIuiKDlUwrIEM5vM5Kf2P2Ftas3Iv0by/t73CYwNNBpvZ1+CBh/8yPVXthFhUpp6ZyZxcVZTrl84liFOY2mJ04gRVPh9K5Y1a6IplvocTGH5JflFVJhHKCakTsq3Bm6TOin/qpTyQl59phqhKMqzlWJIYc3FNSw+sxidQccQ7yG85f0W5lrjT8Ub9HpObP6GSufmYCPjOVG6L9UHzMLGNuPDjlLKtBFP+MJFpATfwWnsWPVQZD4p1CMUIcRa4AhQWQgRJIQYLKXUAR8Au4CLwK95WUwURXn2TDWmDPIaxNbuW2lVrhWL/11M983dORBkfNpTo9VSv+co+OAkpxw70TB0LfHzanPij+Wpd3o9kOHymcFAzJatXGvfgah1v2SIU/JXoRmhFAQ1QlGUgnU0+Cgzj87kRuwNWrm2YkL9CZSxKZNl/OUT+9DuGIuH/hrnzGth+8rXlK/skykuyd+fkOmfEX/8OBbe3jjPmIFF5awftlSeTqEeoSiK8mJq6NyQjV03MqL2CI4EH6Hb5m4sP7ecFH2K0fjKdVvh/tFxjlX9iPJJl3Fe04ojyz4k/n5MhjhzDw/KrfyJMl99hS4sDAz6Z3E6Lzw1QlEjFEUpFO7cv8MXx79gX+A+3GzdmNxwMg2dG2YZoCifbQAAGCFJREFUHxkaSMCaMdSL2UUIJbnTaBq12r6a6W4wmZyMeNAHLHTWbMwrV8aue7dMcUr2qRGKoiiFWhmbMixotYBFrRehM+h4+8+3Gbd/HKFxoUbjHUu5Um/Ur/i1/4UEjRW1j7zP2a/a/b+9O4+OqsgXOP79dWchAglICCB7FJCdsOOIZ4SngE92AQ1uT5Eni4AjjCDMM44CyjIDzphxZBEdQEERWV1mQAcXiIQQJiwDQgImQEhYAiF7uuv9kWYmy22cQCfpJL/POTmHrv7de6vqnObXVbdvFafji95qvZZMnDk5ZMXFcfallzj16GNkHz1a5m2qbnSEoiMUpbxOdn42Kw+uZEXcCnxsPkzsMpHwtuH42qyXWsnLzWHfR6/T8VgkPjiIafY/hIVHUCOgZpE443RyeeOnpCxahOPKFeqODaf+lCnYa9Uqj2ZVGV6/9EpF0ISilHdLvJLIvB/m8e3pb7mjzh3M6T2Hbg26uY1POZ1A4ofP0y39K05LA873fY3O/UaXiHOkpZGyZAlXtm0ndNtWfENCyrIZVY4mFAuaUJTyfsYYdibu5I0f3uBsxlmG3D6E57s9T3BAsNtjDn6zicCvZtHMeZr9t/yChmN+T6PmbUrEOdLSsNepgzGGlIWLqDNyBP63316WzakSNKFY0ISiVOWRmZfJsrhlrDq0igB7AM91fY7RrUdjt9kt43Ozs4hZ9yqd45cBEBv6DN0e/j/8/EvuV5976hQJo8fgzMig3pNPEDxhAraaNUvEqQKaUCxoQlGq8om/HM+8qHlEnY2i7a1tmd17Np3rd3Ybn/zTMc6um0ZYxnf8ZGvMlXvn06FvyTW/8i9eJGXRYi5/8gk+jRrRYOZMat9/n9s1x6ozTSgWNKEoVTkZY/ji5Bcs3LuQlKwURrYaybSu06hTo47bYw58tZ7gXXNobM6xr/a9NH3494Q0blkiLjMmhuRXfoszPZ3Qz7Zj87deFqY604RiQROKUpVbRl4GkbGRrDmyhtp+tZnWdRrDWw3HJtZPRGRnXiX2gwjCflpFPnbiWk+k26iZ+PoVTRomP5+806fxa94cZ04Ol1avoe7YcGw1Sk6XVUeaUCxoQlGqajh26Rhz98wlJiWGTvU7MafXHNrWa+s2/nT8IS6sn0qn7L0k2JqTdf8C2vUeaBl75csvOT1lKr5NmtBg9kvUvvfesmpGpaEJxYImFKWqDmMMW+K3sDh6MWk5aYxpM4bJYZMJ9Au0jnc6if3bWhp9H0FDUtkbNICWjywmuGHTErEZe6JIfvVVck+coFa/fjR46SX8mlxvv7+qTROKBU0oSlU9V3Kv8IeYP7D+2Hrq+Nfhhe4vMDh0sNub65lXL3Ng7W/odno12eLPkbZT6T5yOnYfnyJxJjeXi++/T+pbkQR06kTz91aVQ2u8kyYUC5pQlKq6Dl84zNw9c/nH+X/QNaQrc3rPoVXdVm7jTx2NJf2TqXTIieW4/XbyBy3izu79SsTlnT2LMysL/9BQ8i9cIPvIP6l19y/KsileR9fyUkpVK+3qteMvD/yFl/u8zInLJxi1ZRQL9y4kIy/DMr55my60f/Er9vVYTKDjEq23jCBq6aOknU8uEufbqBH+oaEAXHz3XRLHjSNp6jTykpOtTlutVNoRioiEArOBIGPMQ66ymkAkkAt8bYxZc71z6AhFqerhUvYllsYsZcOPGwgJCGFGjxkMaDHA7TTY1SuXOLhmFt2T15EuNTnW8QV6DJuCzV70IUpnTg4XVqzgwp/fAbud+pMmcuvjjyO+1muOVRVeNUIRkZUikiIiB4uVDxSRoyJyXERmXu8cxph4Y8zTxYpHAB8bY54Bhni42kqpSqpujbpE3BXB6gdWUy+gHjN2zWD8X8eTcDnBMr5WYF16T3ibxNGfc9a3Ob3iIjg2/y6OH/iuSJzN35/6EycSum0rNXv2JGXhIlJ+9/vyaJJXqpARiojcA1wF3i+0h7ydgj3k7wOSKNhD/hHADswvdoqnjDEpruM+LjRCmQV8ZoyJFZG1xpjw69VDRyhKVT8Op4N1R9fxx/1/JMuRxZPtn+SZjs9wi+8tlvHG6SR6858IjX2DOuYK0fVHcOfYBQTVLbmWWPrOndRo3x7fBg3ITUrC5u+PT/36Zd2kcudVIxRjzC7gYrHinsBx18gjF/gQGGqMiTPGPFjsL8XNqZOAJq5/6/0hpVQJdpud8LbhbB6+mUEtBrE8bjnDNg1jx087sPqCLTYbPYZNwmdqDNH1h9M99RPyl3Zl76dvldivvna/fvg2aABAcsQrnBj0ABff/wsmP79c2lbRvOk/3cZAYqHXSa4ySyJST0TeBsJcIxOAT4CRIvInYIub48aLSLSIRKempnqo6kqpyiY4IJh5fefx7oB3qelbk2lfTWPSjkkkXkm0jA+qG0yvye+SMGIrF3wa0CP2JY7M70vC4b2W8Q1mv0RA586cmzePhFGjydy/vyyb4xUq7Ka8iLQAthaa8hoFDDDGjHO9fgzoaYx5rqzqoFNeSimAPGcea4+sJTI2knxnPuM6juOpjk/hb7dex8vpcBC9cSmtDy6mlskkuuEYOoydT63AukXijDGkf/EF5+a/Tv65czResoTAgQPKo0llyqumvNxIAgo/otoEOFNBdVFKVSO+Nl+eaP8Em4dtpl+zfkQeiGT4puF8k/SNZbzNbqfnQ7+CSdHE1HuA3uc+IPN3Xdm3bXmRaTARIXDgQEK3bSN48mRq3dMXKHiWxTgc5dK28uRNIxQfCm7K9wdOU3BTPtwYc8jdOW6WjlCUUlZ2n9nNvKh5nLxykv7N+vNijxdpVKuR2/ij0TuwfzaDOxwniPMPI2jkEpq17mIZ68zNJWHwEGyBgTR8+WUCOrQvq2aUGa8aoYjIB8BuoI2IJInI08aYfGAy8AVwBFhflslEKaXc6XNbHzYM2cDUrlP57vR3DN00lOVxy8lz5FnGt+nen5azfiDqzpk0zzlKwzX92P3OFDKvXi4RK76+BE+eRN7Zs5wcNYqzr7yC43LJuMqo0j7Y6Ak6QlFK/ZwzV8/wxg9vsDNxJy2DWjK712x6NerlNv58ciIJH/yKHpe/JJn6nL0rgi7/FY7Yin5/d6Snk/rmH7i0Zg32OnVovmY1/i1L7s/ijXQtLwuaUJRS/6ldSbuYHzWfpKtJDGoxiOk9phNyS4jb+CO7P6PGX1+kpfMUBwJ6EjxqKY1D25WIyz58mEvr1tPw/36D2O04rmZgr+Xd2w9rQrGgCUUpVRrZ+dmsPLiSFXEr8LH5MLHLRMLbhuNrs15qJS83h30fvU7HY5H44CCm+f8Q9kgENQKsE0b+pUvE//eDBA1+kODnnsNeq1ZZNueGedU9FKWUqoxq+NRgYpeJfDr0U7o26Mqi6EWM2TqGfef2Wcb7+vnTe+zLZIzfw8HAu+nz0zucX9CVAzs/sowXm43a993Hxff/QvygB7i8dZvlw5beSkcoOkJRSt0AYww7f9rJ63tfJzkjmSG3D+H5bs8THFBySZZr4nZtIujrWTRznmZ/zbtpNGYJDZuVXFI/Ky6O5IhXyD50iFt696bp23/yqu2HdcrLgiYUpdTNyszL5J1/vMN7h98jwB7Ac12fY3Tr0dhtdsv4nOxMYta9Rpf4ZRiEA6HP0O3h3+DnXzRhGIeDtI8+IvvQYRq9+tt/lYnd+rzlSROKBU0oSilPib8cz7w984hKjqLtrW2Z03sOnep3cht/9tRRktc9T1jmd5yyNSG933w63O1+kfSc+HgSn51AyPQXqH3ffW6X3i8Peg9FKaXKUGhQKMvuX8aCexZwPus8j25/lIjvI0jLTrOMb9S8DWG/3s6Be5bhY/Lo8LfH2Ld4OKlnTlrGm+xsbAEBnJ4ylcTx/0vuqVNl2JoboyMUHaEopTzsau5VIg9EsvbIWmr71WZa12kMbzUcm1h/h8/OvMr+DyLo+tMq8rET13oi3UbNxNev6FpiJj+fS2vWkPrmHzB5eQRPeJbgCRPKo0lF6JSXBU0oSqmydOzSMebumUtMSgyd6ndiTq85tK3X1m386RMHOf/RNDpn7yXB1oKs+9+gXe+BJeLyzqWQsmABvo0aEjJ9elk2wZImFAuaUJRSZc0Yw5b4LSyOXkxaThpj2oxhcthkAv0CreOdTmL/toZG30fQkPPsDRpAy0cWE9ywqWWs2Gxk7N7NxdVraDBrFn5N3O764TF6D0UppSqAiDDk9iFsGb6F0a1Hs+7oOgZvHMyWE1vcbugVdv9jBE6PYfdtj9M57W/4vd2TqHWv4yi2Ude15VzyziaT8f33xD/4IOff/jPO3NxyaVtxOkLREYpSqhwdvnCY1/a8Rtz5OLo16MbsXrNpVbfksyjXnDoaS/onU+mQE8tx++04Bi2iTfd+JeLyzpzh3OtvkP7ll/i1aEHDiAhq9na/5tjN0BGKUkp5gXb12rH6gdW83OdljqcdZ9SWUSzau4iMvAzL+OZtutD+xa+I7rGIIMdFWm0ZwQ9vPkra+eQicb633UaTN5fSdNk7GOMkNyG+PJpThI5QdISilKogl7IvsTRmKRt+3EBIQAgzes5gQPMBbp8xSb98kUNrZ9E9eT3pUpMfO75A92FTsBV72NGZk4P4+CB2O5c3bSL//AVuffwxxNd6zbHS0hGKUkp5mbo16hJxVwSrH1hNvYB6zPj7DMb/dTwJlxMs42sH3UrvCX8mcfTnJPs2o2dcBD/Ov4vjB74rEmfz9//XE/UZu/eQsnAhCSNGkLl3b5m2p9KOUEQkFJgNBBljHipUXhPYBbxsjNl6vXPoCEUp5S0cTgfrjq7jj/v/SJYjiyfbP8n4TuMJ8AmwjDdOJ9Gb/0Ro7BvUMVeIrj+CO8cuIKhuybXE0nfu5Nxrc8k7c4agoUMImTEDn2D3a479HK8aoYjIShFJEZGDxcoHishRETkuIjOvdw5jTLwx5mmLt14E1nuyvkopVdbsNjvhbcPZPHwzg1oMYnnccoZ+OpQdP+1w+2uwHsMm4TM1huj6w+me+gl5S7uxd1NkkX3tAWr360fotq3Ue/Z/ubL9M7IPHy6TNlTICEVE7gGuAu8X2lPeTsGe8vcBSRTsKf8IYAfmFzvFU8aYFNdxH18boYjIfwHBQA3gvI5QlFKVVXRyNHOj5nI87Th9G/dlVq9ZNK1d8lmUa44f+BbnludpnX+Mw34dCRj2e1q261EiLu/cOXwbNLipunndg40i0gLYWiih9AEijDEDXK9nARhjiieT4ucpnFDmAjWBdkAWMNwY43R3rCYUpZQ3y3PmsfbIWiJjI8l35jOu4zie6vgU/nZ/y3inw0H0xqW0PriYWiaT6IZj6DB2PrUC63q0Xl415eVGYyCx0OskV5klEaknIm8DYYWSz2xjzDRgLbDMKpmIyHgRiRaR6NTUVM+2QCmlPMjX5ssT7Z9g87DN3NvsXiIPRDJ803C+SfrGMt5mt9PzoV9hJkUTc+sD9D73AZm/68q+7StKTIOVBW8aoYwCBhhjxrlePwb0NMY8V1Z10BGKUqoy2X1mN/Oi5nHyykn6N+vPiz1epFGtRm7j/xm9A5/PZnCH4wRx/mEEjVxCs9ZdbroelWGEkgQUniBsApypoLoopZTX6XNbHzYM2cCUsCl8d/o7hm4ayvK45eQ58izj7+zen5azfiDqzpk0zzlKwzX92P3OFLIy0sukft6UUPYCrUSkpYj4AQ8Dmyu4Tkop5VX87H480+kZPh32KX0a9WFpzFJGbhlJ1Nkoy3i7jw+9Hp5F7rM/cKBOf/qceY/LC8M4vOdzj9eton42/AGwG2gjIkki8rQxJh+YDHwBHAHWG2MOVUT9lFLK2zWu1Zil/ZbyVv+3yHPkMe7Lcfz6778mJTPFMj64YVN6PP8Rhwd8yFV7IDXrNvR4nSrtg42eoPdQlFJVQXZ+NisPrmRF3Ap87b5M7DyR8Lbh+Nh8LOOvLXt/oyrDPRSllFI3oIZPDSZ2mcjGoRsJCwljYfRCRm8dTcy5GMv4m0km16MJRSmlqohmgc2I7B/Jkl8uIT03nSc+f4LZ387mQtaFcrm+JhSllKpCRIT+zfuzaegmnu7wNNsTtjP408F8+M8PcTgdZXptTShKKVUF3eJ7C9O6TWPDkA20u7Udc6PmEr49nLjUuDK7piYUpZSqwkKDQll2/zIW3LOA1MxUxm4fyyu7XyEtO83j19KEopRSVZyIMKjlIDYP28yj7R5l448biT7n+V+46s+G9WfDSqlqJvFKIk1qN3G7M+TPcfezYesfKSullKqymga6Xwb/ZuiUl1JKKY/QhKKUUsojNKEopZTyCE0oSimlPEITilJKKY/QhKKUUsojNKEopZTyiGr9YKOIpAKnKroeHhYMnK/oSlQi2l+lo/1VelWxz5obY+oXL6zWCaUqEpFoqydYlTXtr9LR/iq96tRnOuWllFLKIzShKKWU8ghNKFXPOxVdgUpG+6t0tL9Kr9r0md5DUUop5RE6QlFKKeURmlCUUkp5hCYUpZRSHqEJpRoRkbYi8raIfCwiEyq6Pt5OREJFZIWIfFzRdfFW2kelU9U/g5pQKgkRWSkiKSJysFj5QBE5KiLHRWTm9c5hjDlijHkWGA1U6QetPNRf8caYp8u2pt6nNH1XXfuosFL2V5X+DGpCqTxWAQMLF4iIHXgLGAS0Ax4RkXYi0lFEthb7C3EdMwT4FthRvtUvd6vwQH9VU6v4D/uu/KvmlVZRiv6qyp9B3VO+kjDG7BKRFsWKewLHjTHxACLyITDUGDMfeNDNeTYDm0VkG7C27GpcsTzVX9VRafoOOFy+tfM+pe2vqvwZ1BFK5dYYSCz0OslVZklEfikib4rIn4HtZV05L1Ta/qonIm8DYSIyq6wr5+Us+077yC13/VWlP4M6QqncxKLM7ZOqxpivga/LqjKVQGn76wLwbNlVp1Kx7DvtI7fc9dfXVOHPoI5QKrckoGmh102AMxVUl8pA++vGad+VTrXsL00oldteoJWItBQRP+BhYHMF18mbaX/dOO270qmW/aUJpZIQkQ+A3UAbEUkSkaeNMfnAZOAL4Aiw3hhzqCLr6S20v26c9l3paH/9my4OqZRSyiN0hKKUUsojNKEopZTyCE0oSimlPEITilJKKY/QhKKUUsojNKEopZTyCE0oSimlPEITilJKKY/QhKKUckt3ZFSloQlFVUsi4hCR2EJ/LSq6Tp4gIi1EJEtEYj1xPqsdGUUkwNVnuSIS7InrqKpBl69X1VWWMaaL1RsiIhQsS+Qs5zp5ygl3bXNHRDoC84sVP2WMSSkea4zJArqIyMkbr6KqijShKEXBN3vgM+AroA8wTET6AlMAPyAKmGiMcbjiZwOPU7CJUiqwD/gY2GqM6eCKmQ7UMsZEiMijxc9FwfLmn1GwHexdwGkKdpDMch3/ODCdgj1b/gGcBM4bY5a63p8LnDPGvPkz7frcdY3ewAHgXeAVIAQYa4z5wRgTh+5aqW6STnmp6iqg0HTXRldZG+B9Y0wYcAswBviF69u+AxgLICLdKFiOPAwYAfS43oVEpK27cwGtgLeMMe2BNGCk65j2wGygnzGmMzAVWAE84Xrf5qrDmv+grXcAS4FOwJ1AOHA3BcnqpZ+pu+7IqP5jOkJR1VWRKS/XN/lTxpg9rqL+QDdgb8EMGAHAtemfvsBGY0ym69if2+fC3bl2AQnGmGv3O/YBLVz/7gd8bIw5D2CMuQhcFJELIhIGNAD2u3ZM/DkJrhEIInII2GGMMSISV+h6lnRHRlUamlCU+reMQv8W4D1jjLtv5Vb7PuRTdNRf43rnciWxnEJFDgqSzbVjrK6xHHgSaAisdFO34gpfw1notRP9P0B5kE55KWVtB/CQiIQAiMitItLc9d4uYLjr1061gcGu8nNAiGuayJ9/35O43rmud/3RIlLv2jGu8o3AQAqm2b646VYq5UH67UQpC8aYwyIyB/jSdb8iD5hEwbRYjIisA2KBU8A3rmPyROS3FNx0TwD++TPnSr7O9Q+5brr/XUQcwH7gSWNMroh8BaRd+4GAUt5Cd2xU6iaJSARw1RizqByuZQNigFHGmB8t3m9BoV+alXFdTgLdr93nUUqnvJSqJESkHXCcgpvqJZKJiwMI8tSDjW7qEeA6vy8F92GUAnSEopRSykN0hKKUUsojNKEopZTyCE0oSimlPEITilJKKY/QhKKUUsojNKEopZTyCE0oSimlPEITilJKKY/QhKKUUsoj/h8ObQTuzzmvegAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "r0 = np.array([0.1, 0.5])\n", "L0 = [25, 25]\n", "psd_kolm = VonKarmanPsd(0.1, np.inf)\n", "psd_vk = VonKarmanPsd(r0, L0)\n", "spatial_freqs = np.logspace(-4, 4, 1000)\n", "psd_kolm.plot_von_karman_psd_vs_frequency(spatial_freqs)\n", "psd_vk.plot_von_karman_psd_vs_frequency(spatial_freqs, idx=0)\n", "psd_vk.plot_von_karman_psd_vs_frequency(spatial_freqs, idx=1)\n", "plt.plot(spatial_freqs, 10*spatial_freqs**(-11/3), '--')\n", "plt.legend(['kolm r0=0.1m', 'von karman r0=0.1 L0=25', 'von karman r0=0.5 L0=25', 'k^-11/3'])\n" ] } ], "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.6.10" } }, "nbformat": 4, "nbformat_minor": 2 }