{ "cells": [ { "cell_type": "markdown", "id": "9c6e42df-f551-4130-9a07-cfbf05581c42", "metadata": {}, "source": [ "# Orbital Mechanics (Kepler's Laws)" ] }, { "cell_type": "markdown", "id": "eccbe0b6-b9c3-4a67-9e00-4f806200c09b", "metadata": {}, "source": [ "In this lecture, we are going to discuss\n", "\n", "* Conservation law for central forces\n", "* Derivation of Kepler's 2nd Law\n", "* Derivation of Kepler's 1st Law\n", "* Properties of orbits and velocities at perihelion and aphelion" ] }, { "cell_type": "markdown", "id": "af1cfe5e-36b9-4f68-8db8-44c7857641f0", "metadata": {}, "source": [ "
\n", " Key Question \n", "Now that we know how Newton's theory of gravity works for non-point masses, let's use this theory to see if we can derive Kepler's laws.\n", "
" ] }, { "cell_type": "markdown", "id": "d4f796be-a35b-4d16-b14f-d00225023046", "metadata": {}, "source": [ "\n", "## Angular momentum & central forces (recap of PY2101)\n", "So, now that we know that large bodies can be treated as if they were point masses, what can we derive from it? First let's consider two masses, $M$ and $m$, located at the positions below.\n", "\n", "![Setup](Figures/Cartesian_Polar_coords.png)\n", "\n", "While the above has the x and y axes labelled, it's far more convenient to work in polar coordinates. Given that\n", "\n", "\\begin{align}\n", " \\hat{\\textbf{r}} &= \\cos \\theta \\hat{\\textbf{i}} + \\sin \\theta \\hat{\\textbf{j}}\\\\\n", " \\hat{\\pmb{\\theta}} &= -\\sin \\theta \\hat{\\textbf{i}} + \\cos \\theta \\hat{\\textbf{j}}\n", "\\end{align}\n", "\n", "then we know\n", "\n", "\\begin{align}\n", " \\hat{\\textbf{r}}\\cdot\\hat{\\pmb{\\theta}}&=0\\\\\n", " \\hat{\\textbf{r}} \\times \\hat{\\pmb{\\theta}}&=\\hat{\\textbf{k}}\n", "\\end{align}\n", "\n", "From these definitions, we can derive the following time derivatives:\n", "\n", "\\begin{align}\n", " \\frac{ {\\rm d} \\hat{\\textbf{r}} } { {\\rm d} t} &= \\hat{\\pmb{\\theta}} \\frac{ {\\rm d } \\theta } { {\\rm d } t }\\\\\n", " \\frac{ {\\rm d} \\hat{\\pmb{\\theta}} } { {\\rm d} t} &= - \\hat{\\textbf{r}} \\frac{ {\\rm d } \\theta } { {\\rm d } t }\n", "\\end{align}\n", "\n", "The whole point of the above is that the velocity of the planet can then be written as\n", "\n", "\\begin{align}\n", " \\textbf{v} &= \\frac{ {\\rm d} \\textbf{r}} { {\\rm d} t} = \\frac{ {\\rm d} r } { {\\rm d} t} \\hat{\\textbf{r}} + r \\frac{ {\\rm d} \\hat{\\textbf{r}} } { {\\rm d} t}\\\\\n", " &= v_r \\hat{\\textbf{r}} + v_t \\hat{ \\pmb{\\theta}}\n", "\\end{align}\n", "\n", "The angular momentum of the planet is then\n", "$$\n", " \\textbf{L} = \\textbf{r} \\times \\textbf{p}\n", "$$\n", "which, when you take the derivative with respect to time, gives\n", "$$\n", " \\frac{ {\\rm d} \\textbf{L} }{ {\\rm d} t } = \\frac{ {\\rm d} \\textbf{r}}{ {\\rm d} t } \\times \\textbf{p} + \\textbf{r} \\times \\frac{ {\\rm d} \\textbf{p}}{ {\\rm d} t }\n", "$$\n", "using $\\textbf{p}=m\\textbf{v}$ and $\\textbf{F}= m\\frac{{\\rm d} \\textbf{v}}{{\\rm d} t}$, we thus get\n", "$$\n", " \\frac{ {\\rm d} \\textbf{L} }{ {\\rm d} t } = m (\\textbf{v}\\times\\textbf{v})+\\textbf{r}\\times\\textbf{F}\n", "$$\n", "The first term has to be 0 (a vector crossed with itself is 0) while for that last term must be zero as $\\textbf{F}$ is parallel to $\\textbf{r}$. Thus, for any central force\n", "$$\n", " \\frac{ {\\rm d} \\textbf{L}} {{\\rm d} t} = 0\n", "$$\n", "\n", "## Kepler's Second Law\n", "So, if angular momentum is convered under gravity, what does that tell us about the orbit of (for example) a planet around a star? The angular momentum of the planet can be written as\n", "\n", "\\begin{align}\n", " \\textbf{L} &= \\textbf{r} \\times \\textbf{p} \\\\\n", " \\textbf{L} &= m r v_t \\hat{\\textbf{k}}\n", "\\end{align}\n", "\n", "Now imagine a setup like that shown below, where we have a planet moving with a velocity $\\textbf{v} = v_r \\hat{\\textbf{r}} + v_t \\hat{ \\pmb{\\theta}}$. We're interested in how much area is swept out by the orbit of the planet in a time $\\Delta t$.\n", "![Planet_Vel](Figures/Planet_velocity.png)\n", "This is given by\n", "$$\n", " \\Delta A \\approx \\frac{1}{2} r v_t \\Delta t + \\frac{1}{2} r v_r \\Delta t\n", "$$\n", "In the limit where $r>>v_r \\Delta t$ (which is true for planets when considering very short timescales), we end up with\n", "$$\n", " \\Delta A \\approx \\frac{1}{2} r v_t \\Delta t\n", "$$\n", "In the limit $t \\to 0$ we then get\n", "$$\n", " \\lim_{t \\to 0} \\frac{\\Delta A} {\\Delta t} = \\frac{{\\rm d} A}{{\\rm d} t} = \\frac{1}{2} r v_t = \\frac{L}{2 m}\n", "$$\n", "The important thing about this last term is that both $L$ and $m$ are constants, which means\n", "$$\n", " \\frac{{\\rm d} A}{{\\rm d} t} = {\\rm constant}.\n", "$$\n", "This is Kepler's Second Law.\n", "\n", "## Kepler's First Law\n", "The above expressions can be used to state that\n", "$$\n", " L = r^2 \\frac{{\\rm d} \\theta}{{\\rm d} t}\n", "$$\n", "which, following Kepler's Second Law, is constant for any central force. When that force is a gravitational force, we have that\n", "$$\n", " \\textbf{F} = -\\frac{GMm}{r^2}\\hat{\\textbf{r}} = m\\frac{{\\rm d} \\textbf{v}}{{\\rm d}t }\n", "$$\n", "Following some simple math (Section 3.1.2 of Ryden & Peterson), we arrive at the following expression\n", "$$\n", " r = \\frac{L^2}{G M m^2 (1+e \\cos \\theta)}\n", "$$\n", "It can be shown that this expression can be written as\n", "$$\n", " r = \\frac{a(1-e^2)}{1+e \\cos \\theta}\n", "$$\n", "where $a$ is the semimajor axis of the orbit, $e$ is the eccentricity of the orbit, and $\\theta$ is called the \"true anamoly\".\n", "\n", "![Ellipse](Figures/Ellipse_Definitions.png)\n", "\n", "As shown below, this is the equation for a conic section.\n", "\n", "![Conic_Sections](Figures/Conic_Sections.png)" ] }, { "cell_type": "code", "execution_count": 1, "id": "0855ec8b-ef9e-46db-95fc-05f03fcab569", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/tmp/ipykernel_14232/959222672.py:8: RuntimeWarning: invalid value encountered in true_divide\n", " r = (1-e**2)/(1+e*np.cos(theta))\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAF5CAYAAADQ7F7EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAABcSAAAXEgFnn9JSAABMzklEQVR4nO3dd3wUdf7H8ddmN70X0kNLaAlJKAlKR0RFBCkidgRFRUU99e701Dv1fp6n5+mdFbsodgREUU4QRaQTILTQEkoCaQTS+2bn98eGnQSSkECS2fJ5Ph55MN/JzOYDJHnvfOc7369OURSEEEIIrTlpXYAQQggBEkhCCCGshASSEEIIqyCBJIQQwipIIAkhhLAKEkhCCCGsggSSEEIIqyCBJIQQwipIIAkhhLAKEkhCCCGsggSSEEIIqyCBJIQQwipIIAkhhLAKho54UZ1Olwt4AFkd8fpCCCGsVhRQoShKaFtP1HXE8hM6na7E1dXVOzo6ut1fWwhh5xQFio6Bsap+hw78uoKzu6ZlidbJyMigurq6VFEUn7ae2yFXSEBWdHR07N69ezvo5YUQdklRYPGdsOcE4GXed+3rMGimpmWJ1ouLiyMtLe2CesfkHpIQwnqsewX2LFbbl8yVMHIgEkhCCOuw/0dY/X9qu+cYuPIfmpUjOp8EkhBCe/n7YMldQP09bf8eMP0j0HfUXQVhjSSQhBDaqjgNX9wINWXmtos33PQleARoW5fodBJIQgjt1NXC1zOh8Gj9Dh1c9z4E99WyKqERCSQhhHZ+egKO/q62xz0NfcZrV4/QlASSEEIbKR/BlnfVdvz1MPwPmpUjtCeBJITofMc2wo9/VNvhA83PG+l02tUkNCeBJIToXCXZ5vtGJqO57RUCN34uMzEICSQhRCcyVsNXt0F5vrmtd4EbPgOfcG3rElZBAkkI0Xl+/BOcSFHbE16CqGTt6hFWRQJJCNE5Uj6C7R+r7cGzzB9C1JNAEkJ0vKwt5qujMyKT4ep/aVePsEoSSEKIjlWaa75vZKo1tz2DYcZCMLhqW5ewOhJIQoiOY6yBr2+Hslxz28kAMz4BnzBt6xJWSQJJCNFxfvoLZG1S2+NfgG5DtatHWDUJJCFEx9jxKWx9X20PuBWS52hXj7B6EkhCiPZ3Yhssf0Rthw+Ea16WmRhEiySQhBDtq+ykeRBDXbW57REEN3wKzm7a1iWsngSSEKL91NXCollQcsLc1unh+gXgG6llVcJGSCAJIdrP6mfh2Dq1feVz0GOkdvUImyKBJIRoH2nfwYbX1Xb89XDpvdrVI2yOBJIQ4uKdyoBl96vtLv1g0qsyiEG0iQSSEOLi1FSYBzFUl5jbLt5ww0Jw8dS2LmFzJJCEEBdOUeCHRyB/r7pv8hsQ1Eu7moTNkkASQly4bQtg5xdq+9L7IW6KVtUIGyeBJIS4MNk7YMWf1XbUpXDFs9rVI2yeBJIQou0qTpuXIa+rMbc9u5ifN9I7a1qWsG0SSEKItjGZYOk9UJRpbuucYPqHMoO3uGgSSEKItln3MhxaqbbH/hV6jNKuHmE3JJCEEK13eA38+rza7n01DP+DVtUIOyOBJIRondI8WHwXKCZz2787TJ0PTvJrRLQP+U4SQpyfqQ6WzIHyfHNb72pe+dXdX9u6hF2RQBJCnN/af8ORtWp7/PMQlqhdPcIuSSAJIVp25Hf47QW1HTsZku7Urh5htySQhBDNKzsJi+c0vm907esyaaroEBJIQoimmUyw9G4oyzW3nZxh+kfg5qttXcJuSSAJIZq2/j+Q8YvavvI5iBikXT3C7kkgCSHOdWwj/PIPtd13Ilxyj3b1CIcggSSEaKz8FHxzByh15rZvV/OSEnLfSHQwCSQhhEpRYNl9UJptbjsZ4PqP5Hkj0SkkkIQQqi3vwcH/qe1xz0BkkmblCMcigSSEMMtLg5VPqe2YK2DoPO3qEQ5HAkkIAbWVsPhOqKs2tz27wJT5ct9IdCoJJCEErPob5Kep7Slvg1cX7eoRDkkCSQhHd/An2PKu2r70Pug1Trt6hMOSQBLCkZXmwrf3qu2QePNABiE0IIEkhKMymcxhVHHK3Da4wXXvg8FV27qEw5JAEsJRbZ7feGqgq56H4L7a1SMcngSSEI4oZyeselpt97kGku7Qrh4hkEASwvHUVJiXlDDVmtteobKkhLAKEkhCOJqfn4GCg/UNHUx7BzwDtaxICEACSQjHkvErbHlHbQ+bBz3HaFaOEA1JIAnhKCqLYNn9ajs4Fi57qtnDhehsEkhCOIr/PQ4lJ8zbTs4w9R1wdtO2JiEakEASwhGkfQc7v1DbYx6HsATt6hGiCRJIQti7snxY/ge1HZkMw//Q3NFCaEYCSQh7pijw3YPqbAzOHuauOr1B27qEaIIEkhD2LPUzOLhCbV/xdwiM1q4eIVoggSSEvSo8BiseV9s9L4PkOdrVI8R5SCAJYY8UBb57AGpKzW03X5j8pszGIKyaBJIQ9mj7J3DkN7U94d/gG6FdPUK0ggSSEPam+ASsbPDAa58JEH+9dvUI0UoSSELYE0UxD/GuLjG3XX3hmlekq07YBAkkIezJrq/h0Eq1Pf558AnTrh4h2kACSQh7UZoHK/6stqPHwoBbtKtHiDaSQBLCXvz4R6gqMm+7eMGkV6WrTtgUCSQh7MHeb2Hfd2r7imfBr6tm5QhxISSQhLB1lYXmq6Mzuo2AwbIcubA9EkhC2Lqfn4Hyk+Ztgztc+xo4yY+2sD3yXSuELcvcDNsWqO0xj8lcdcJmSSAJYavqahsvKxEcC0PnaVaOEBdLAkkIW7XxDchPU9sT/wt6Z83KEeJiSSAJYYsKj8KaF9X24NnQ9RLNyhGiPUggCWFrFAV++CMYK81tzy4w7mltaxKiHUggCWFr0r6F9FVqe/wL4O6vWTlCtBcJJCFsSVUxrHhMbfe8DPpfp109QrQjCSQhbMmaF6Esz7xtcIOJMpO3sB8SSELYivz9sOUdtT3yUQjoqV09QrQzg9YFCNEe6kwKJ0uryS+toqiiluLKsz4qaimvMVJtNFFtNFFjrDNv15qorTOh04EOnflPnQ4d5skO3Ax63F30uDnr8XDR4+5sbvu6O+Pv4YKfhzN+Hi74ezgT4OlCsLcbLoYOeJ+nKOaZvE1Gc9u/Owx7sP2/jhAakkASNqGqto5jpyo4UlDO0VPl5BRVkltSRW5xFbklVZwsrcakaF2lWZCXCyE+boT6uBHi60akvzs9Aj3pFuhJ9yAPPFwu4Mdu3/eNlyS/6p/g7NZ+RQthBSSQhFUpqaplX3YJaTklpOeXcfRUOUcLKsgurkSxksA5n4KyGgrKatibXdLk54O9Xeke5EmPQE96h3rTL8yb2DAf/Dxcmn7B2kr46Um1HX059Lm6AyoXQlsSSEIzeSVV7D5ezN7sEtJyiknLKSHrdOVFvaa3qwF/T3NXmq+7Mz7uzvi5m7c9XQ24Gpxwddab/zQ44WrQ46w3DwowKaAoCgrmP+tMUG2so7K2jsqaOqpq66ioMX8UVdRQWFFr+bOwoobSKmOraswvrSa/tJotR0432h/m60a/MB/6hnoTH+HLoG7+hPi4wfpXoTjTfJCTwTzMWwYyCDskgSQ6hbHOxP7cUrZnFpJytJBtxwo5UdS28HHW64gK8KBHoCdRAR7mbjFfV0v3WKiv24V1h7WTamMd+SXVlq7EvJIqsouqyDxdwdFT5WSeqqCmztTs+TnFVeQUV/HL/nzLvkG+pXxZ8zJnrp3qhsxF36V3B/9NhNCGBJLoEIqicCCvlPXpp9iQXsDmI6cpq27dFYS3m4HY+iuFHkGedA/ypGeQF+F+bhj01jsw1NWgJyrAg6gAjyY/X2dSyCmutNwLyzhZxv6cUvblllBUUdvkOXdWfICLvgaAk4ovV68fTGz2FkbEBDIsOojYMB+cnORqSdgHCSTRbgrLa1hzMJ9f959kQ0YBBWU15z0nzNeNuHBfYsN9iA3zIS7ch0h/d3R22CWld9IR6e9BpL8Hw2OCLPsVRSG3pIp9OSXsyyklLbuE1Kwiwop3cI1+i+W4F2pvosDkytqDJ1l70Lz+kb+HM0OjAxnduwtj+4bQxdu10/9eQrQXCSRxwRRFIeNkGT/vy2f1vjy2HStscaSbwUlHXIQvg7v6M7ibP4O6+RHm6955BVspnU5HmK87Yb7ujO0bYt6pKNS883+Qa24eMvRmuXEUnNXjV1hRy4+7c/lxdy463W4SI/0Y1y+Yy/uF0DfU2y6DXdgvCSTRZofySvl+Vw7Ld2Vz+GR5s8fpdBAf4cuw6CCGxwSS1C0Adxd9J1Zqw/YuxSV3u6XZ69b/sjPiUrYfK2R9RgHr00+x63hRozcAigKpWUWkZhXx75UHifR355qEMCYlhBMX7iPhJKyeTumAsbQ6nW5vbGxs7N69e9v9tYU2jp0q57vUbJbvyuFAXmmzx3XxduXyvsGM6dOFS3sGNj+UWTTPWA1vDjEvMQHQ5xq46fNzDiupqmXz4dOsPXiS1fvyyC6uavYlewR5MjEhjEmJ4fQO8e6gwoWAuLg40tLS0hRFiWvruRJIolnl1UZ+3J3DopTjbDl6utnj4sJ9uLxfCOP6BdM/3Fdusl+sjW/BT38xb+v0cP9mCOrV4imKorAvp5TV+/L4eV8eO48XN3tsQqQv1ydFcW1iOL7usqCfaF8XE0jSZScaURSF7ZmFfL31OMt3ZVNeU9fkcf0jfJiYEM418WHNjioTF6CyEH5rsPBe0uzzhhGY70PFhvsQG+7DA5f3Iq+kih9357B8Vw7bjhU2OnbX8WJ2HS/mueVpjO8fyoykKIZFB0qXntCcBJIAzFPzfJeazYINR0nLaXqGgd4hXlybGM7EhHC6B3l2coUO4veXoarIvO3iDaMfv6CXCfFxY/bwHswe3oPjhRX8sCuH73dls+eE+n9bbTSxLDWbZanZxAR7cfvQbkwbFImnq/xaENqQLjsHl11UycJNx/hySyaFTTwL4+1q4NoB4cxIiiIh0lfeRXekkmx4dQDUVZvbY/8Ko/7Yrl9if24Ji1KOs3THCU6Xnzss39vVwPVJUcwc2k3edIgLIveQRJsdyC3lrTXpLN+VQ10TY7Uv7RnAjclduSouVEbGdZblD0PKh+Zt7zB4YDu4dEx3aI3RxC/78/hyaxZrDpw85/M6HVzRL4T7L4shMcqvQ2oQ9knuIYlWS80q4s1f01mVlnfO59yd9UwbFMHtw7rLSKzOVngUtn+itkf9scPCCMDF4MT4/mGM7x/G0YJyPtl4jEUpWZTWz6ahKLAyLY+VaXmM7BXEfWNiuLRngFwhiw4lgeQgthw5zWurD7EuveCcz0UFuDPz0u7MSIrC10NGXWlizYvqWke+XWHgzE770t2DPPnbpFgeubI3S7cfZ8GGo2Q0eL7s90MF/H6ogMHd/HlgbAyje3eRYBIdQgLJzqVll/DST/v5tYlumb6h3tx3WQwT+oda9Rxxdu/kQdj1pdoe8xgYOv/5LS9XA7cN7c4tl3Tjl/35vPFrOqlZRZbPbztWyKyPtjKkRwCPje/D4G4BnV6jsG8SSHbq2KlyXll1kO92Zp+zjtDArn7MuyyGsX2D5Z2uNVjzT1Dq5wQKjIGEGzUtx8lJx7jYEC7vF8zGjFO88Ws6GzJOWT6/5chprpu/kXH9QvjTVX3oEyrdu6J9SCDZmeLKWv7780EWbjyG8azBCoO7+fPoFb0ZKs+cWI+8vbB3idoe8xfQW8ePpU6nY1hMEMNigtieWch/Vh3k90Nql+/P+/JYvT+PGYOj+NP4PgR5ycSu4uJYx3e+uGgmk8LXKVn866cD5wzn7RPizZ+u6sPl/eSKyOqs/be6HRwLcdO0q6UFg7r6s/DOS9iQXsCLPx1gZ31XnqLAVylZ/Lg7h4fG9eL2Yd1xlu5fcYEkkOxAalYRTy/bc850MRF+7jx6ZW8mD4hAL9P5WJ+CdNi7VG2P+hM4Wfcv82ExQXwbHchPe/P498oDpOeXAVBabeS5H/bxxZZMnrk2jpG9umhcqbBFEkg2rKSqln/+uJ8vtmQ22u/urGfe2BjmjOyBq0GeIbJa6/8D1HerBvaC2MmaltNaOp2O8f1DGdcvmE83HeOVVQcpqV++PeNkObd9sIVpAyN4amIsAZ4yua5oPQkkG/XL/jyeWLKH3JLGMzxPSgzniQl9ZZ0ha1eUBTsbjKwb8TA42dabB4PeiVnDezApMZyXVx3kiy2ZlgE0S3acYM3Bk/xtYiyTB4RLV7FoFevuHxDnKKqo4ZGvUrljQUqjMOoT4s0Xd13K6zcNlDCyBRtea/zcUcIMbeu5CIFerjw/NZ7v541oNKvD6fIa/vBVKrMXbCWvpPmlMYQ4QwLJhvyyP49xr6xlyY4Tln0ueif+eGVvlj84gqHRgRpWJ1qtLL/xrAzDHwS97T+Q3D/ClyX3DuOvE2Nxd1av9tYcOMlV/13L//bkaFidsAUSSDag2ljHs9/v5Y4FKRSUVVv2J0b5sfzBEcwb20tGNtmSTW+Bsf6KwTMYBt6mbT3tSO+k484RPVj58ChG91YHNhRV1DL30+38+ZudlNVPTyTE2eS3mJU7fLKMaW9t4KP1Ry37XA1OPDGhL0vuHSZzztmamnJI+UhtD5sHzm7a1dNBogI8WDA7mRemxTe6Wvo65TgTXv3dMmxciIYkkKzYku3Hmfj6OvZmq2vY9A31ZvkDI7h7VLQM5bZFO79osN6RFwyepWU1HUqn03HjkK78+NDIRveWMk9XcP3bG/l00zE6YrUBYbskkKxQbZ2Jp5ft4ZGvd1LRYMXW2y7txrf3D6eXXBXZJpMJNs1X2wNvAzdf7erpJD2CPPlm7lAevLwXZ95D1dSZeOrbPTy6aCeVzaxKLByPDPu2MqfLa7jvs21sOnzass/X3ZkXr0tgfP9QDSsTFy19FZxKr2/o4JK7NS2nMznrnXjkit4Mjw5k3hc7OFlqvhe6ZPsJ0rJLePe2JLoGdtxyG8I2yBWSFUnLLuHaN9Y1CqOESF9+fGikhJE92PSWut33GgjoqV0tGrmkZyA/PDCC5O7+ln37c0uZ8tZ6Uo6ebuFM4QgkkKzEyr25XDd/A8cLKy37pgwI5+t7hhLhJ88V2by8vXB4jdq+9F7NStFasI8bn991KXeO6GHZd7q8hpvf38yy1BMtnCnsnQSSFfhySyZzP91GZa25L91JB09M6Mt/bhiAm7NtPb0vmrH1fXU7NAG6DdeuFivgrHfirxNj+e8NA3Cpf2ShxmjioS9TeX31IRns4KDkHpKGFEXhrTUZvPTTAcs+bzcDr980kDF9gjWsTLSr6jLYtUhtX3IPyFQ6AEwZGEG4nzv3LEyhsKIWgJdXHSSvtIq/X9sfJxlJ6lDkCkkjJpPCs9+nNQqjYG9XFs0dKmFkb/YshppS87arr9UuMaGVIT0CWHLfcHoEeVr2fbopk0cX7cRYZ9KwMtHZJJA0UGdSeHTRThZsOGrZ1yPIk8X3DqNvqI92hYmOsa3Bg7CJN4CLjCY7W48gT5bcO4wBDZ5XWrrjBPd+tp2qWhkW7igkkDpZnUnhka9TWdpgPrr4CF8WzR1KVID8orI72amQvUNt2/GDsBfL39OFT+dcwrAGczKuSsvj7oXbqDZKKDkCCaROVGdSePTrVJalZlv2DYsO5Iu7L5Xln+3VtgXqduQQCInTrBRb4OVq4MNZyYzrp3Zbrz14knmf76BWuu/sngRSJ1EUhSeX7ubbs8Log9uT8XKVsSV2qbbSfP/ojKTZ2tViQ9yc9cy/dTDXJIRZ9q1Ky+MPX6XKPSU7J4HUSf710wG+3JplaZ8JI3cXGdZttw78CNX18xC6eEHsFE3LsSXOeif+e8MArogNsez7YVcOjy/ZLUPC7ZgEUid4//fDzF+TYWkndfOXMHIEu75Wt/tdK4MZ2shZ78QbNw9stIzFN9uO859VBzWsSnQkCaQOtmJ3Ds/9sM/S7hvqzQezJIzsXnkBpP+sthNv0K4WG+Zq0PPObYMZ0iPAsu+1X9L5ckumhlWJjiKB1IF2HS/i4a9TLe2oAHc+uWMIvu62vzqoOI89i9Ulyr3DoftIbeuxYW7Oet67LYlewV6WfU9+u4ffDp7UsCrRESSQOkhOcSVzPk6hqtZ8E9bHzcBHs4YQ7GN/i7GJJuz6St2Onw5OckV8MXw9nPlodjJdvM2jUetMCg98vp1jp8o1rky0JwmkDlBjNHHfZ9vJr59i3+CkY/6tg4lp8A5P2LGiTDixTW0nSHdde4j09+CjWcmWFWhLqozcs3AbFTWyJLq9kEDqAM//uI8dmUWW9t8n92d4TJB2BYnOtW+5uh3YS549akf9I3x5cXqCpb0/t5THF8vIO3shgdTOvtuZ3WhKoBuTo7j5kq7aFSQ6377v1e1+k2Qi1XZ2bWI4d41Ul674bmc2n2w8pmFFor1IILWjrNMVPLlkt6UdF+7DM9fKu2OHUpoHmRvVdr9J2tVixx4b37fRFEPP/7iPQ3mlGlYk2oMEUjsxTwu0k9Jqc3+2t6uB+bcMlvWMHM2BH4D67iPfKAgfqGk59sqgd+K1mwZaptyqNpp48MtUmfPOxkkgtZO3f8tgS4MlmJ+b2p+ugfIgpMOR7rpOE+TlykvXq/eT9uWU8PJKeWjWlkkgtYP0/FJe/fmQpX1tYjiTB0RoWJHQRE05HF2ntvtO1K4WB3FZn2BmDu1mab/3+2G2HSvUsCJxMSSQLpLJpPCXJbupqZ/0MdTHjf+b0l/jqoQmjq6DuhrztqsvRF2ibT0O4okJ/YjuYl7cT1HgiSW7ZWZwGyWBdJE+35LJ1qPqO7L/m9JfZmJwVA2nCuo5GvQyi3tncHPW88J1atfdgbxS3v/9iIYViQslgXQRCstr+Nf/9lvaE+JDG81OLBxM+mp1O+Zy7epwQMndA7hpiPp4xaurD5J1ukLDisSFkEC6CP/9+SAlVeqoumcmyRBvh3X6CJxWZ3QnWgKpsz0+vq9l1F1VrYkXG7xZFLZBAukCHcwr5dPN6ozDD1weI/PUObKMX9TtoD7gF6VdLQ7K18OZJyb0tbSX78phe6YMcLAlEkgX6F//20+dyfy8SbdAD24f1l3bgoS2jm1Qt3uO0awMRzdlQARx4T6W9vM/7JNphWyIBNIF2JlVxM/78i3tv1zdD1eDPADrsBSl8ewM3YZpV4uDc3LS8eSEfpZ2yrFCVjf4WRXWTQLpArzSYMXKhEhfroqTgQwOrTgLSk6o7a5DtatFMCwmiMv6qKvMvv5rulwl2QgJpDZKzSpqtDDYw1f0RidP4zu2Yw2ujgJ6gre8QdHaH8b1tmzvzCpiXXqBhtWI1pJAaqP31h62bA/s6seY3l1aOFo4hIbddV2lu84aJEb5MbKXuuTL67+ka1iNaC0JpDbIOl3Bij05lvbc0dFydSTgeIq63VVmZ7AWD4ztZdnecuQ0u48Xa1iNaA0JpDZYsOEo9QPr6Bbowbh+0jXj8Gqr4OQ+tR0xWLtaRCNDegQwsKufpb1w01HNahGtI4HUSlW1dXyz7bilfcfwHuid5OrI4eWngal+CW2Dm/kZJGE1Gk68uiw1m+KKWg2rEecjgdRKK9PyKK40fzO7O+uZNkhm8xZAzk51OyRO5q+zMhPiwwj0dAHMayYt2palcUWiJRJIrfT1VvUb+ZqEMLzdZAJVQeNACkvUrg7RJFeDnhnJ6qwZ36aeaOFooTUJpFbIK6lifYY6bHRGkkwLI+rlqkvWE5rQ/HFCM9c16M3Yc6KE9PwyDasRLZFAaoWf9uZy5rm6SH93krv7a1uQsA6KAgXqwowEx2pXi2hWTLA3sWHqdELf7czWsBrREgmkVlixO9eyPSE+TIZ6C7OyfKhuMJQ4qFfzxwpNTR4Qbtn+cXdOC0cKLUkgncfp8ho2HzllaY/vH6phNcKqnGpwdeQRCB4B2tUiWjQhPsyynZ5fJmslWSkJpPNYn15gefYo2NuVAZF+mtYjrEiBOqchgXJ1ZM2iAjwsy5wDrDkgE65aIwmk81h3SB3MMLJXF5zk2SNxxqkGC/IFxWhXh2iVsX2DLdu/HjjZwpFCKxJILVAUpdGkjCN6BWpYjbA6ReoCjfh316wM0Tqje6uBtOXIact6ZsJ6SCC14ERRJSeKKi3t4dFBLRwtHE5Jg9FaPpHa1SFaZVA3P8vsKmXVRg7klmpckTibBFILdjWYjDEqwF2WKBeNNVwDyVdm7rB2Hi6GRsO/tx07rWE1oikSSC3YmVVk2U6UwQyiobpaKFUfB8BHAskWDO6mPkO4PbNIu0JEkySQWtDwCkkCSTRSmgs0uAfhE97socJ6JET6Wraly876SCC14FCDKUZiw31aOFI4nAr12TRcfcHZXbtaRKv1DvG2bKefLJOBDVZGAqkZJVW1FJRVW9o9GzzDIASVheq2u59mZYi2iQn24sxEKzVGE8dOlWtbkGhEAqkZR06q36geLnpCZUCDaKhRIMnchrbCzVlPlL+HpX3slMzYYE0kkJqRVah+o3YL9JT560RjEkg2K9xPfXOZXVzZwpGis8lqYs3IL1G760J8XDWsRFglB+yyK6gsYMmhJaTkplBuLMfT4ElyaDJTe00lyN12ntEL91Xv9+UUVWlYiTibBFIzTja4fxTsLYEkzlLb4J21i5d2dXSCKmMVL2x5gWUZyzCeWa693sacjby18y2mxEzh8SGP46q3/p+VUF/1CimnWALJmkggNaOgVA2kIC/r/yETnczY4BeZwX7vL1YZq7j353tJyUtp9hijycg3B7/haPFR5o+bj5uV/3v4e7hYtkurajWsRJxN7iE1o7xGfSfo4y7LlYuzGNU3LBjs9w3LC1teaDGMGkrJS+HFrS92cEUXz9tNfR9eWmVs4UjR2SSQmlFZU2fZdnfWa1iJsEqNrpDsM5AKKgtYlrGsTed8m/4tBZUF5z9QQ14NAqmsWgLJmkggNaPh9PQNpxASAoC6GnXbBu6bXIglh5acc8/ofIwmI0sPLe2gitpHwzeYVbV1LRwpOpsEUiuUyGW9OIf9PwaQktu6rrqzbc3d2s6VdBx5msO6SCC1gruLdNmJs+ga/Ogo9vkuu9x4YbMYXOh5QkggNWNEjPpcxaU9AzSsRFglp4aBZNKujg7kabiw6bIu9LzOItPXWS8JpGa4GtR/mhqjff7CERdB1+Cq2WSfV0hJoUkXdF5yaHI7V9K+yhsMZPBwkSdfrIkEUjPcGnTTVdTY5y8ccRGcGvwia+ONf1sxrdc0DE5t+4VtcDIwtdfUDqqofTR89qjhEHChPQmkZgQ0eHjuVFlNC0cKh+SiTtBJjX3eMwlyD2Jy9OQ2nTMlZorVTyPUcJCSj5s8Y2hNJJCaEejVIJDKq1s4UjgkF3VdHWrKmj/Oxj0+5HGSQlrXdZcUksTjQx7v4IouXsM3mL4eEkjWRAKpGYENpgs6WSqBJM7i2mD+umr7XXnUzeDG/HHzmd57erPddwYnA9N7T+ftK962ibnschrM8B3hJwsrWhPpQG1GRIMp6hsuRSEEAK4NrpDsOJDAHEpPD32a+wfcz9JDS9mau9WmZ/vObjChasOlKIT2JJCa0S1QHbp6orCSGqMJF4NcUIp6DQOpqli7OjpRkHsQdyXcxV0Jd2ldykU50eANZpivXCFZE/kN24wofw+c6p/iNilylSTO4tlF3S637rnbhOpUWTUFDe4h9exi3c9MORoJpGa4GJyIClBHUh3Mte9uGdFGnsHqdnk+mORZNVtwIE/9OQ7wdKGLLC1jVSSQWhAX7mPZ3n3CMbplRCt5NQgkkxGqijQrRbRewzeWvUO80MlkdlZFAqkF/SN8LdsSSKIRN1/Qq48GUJanXS2i1VIbzNzfN9Sn+QOFJiSQWpAQ4WfZ3nW8GJNMgiXO0OnAK0Rtl2RrV4totZRjhZbtpO7+GlYimiKB1ILEKF/09SMbiitr2ZdbonFFwqr4dVO3i45pV4doldziKo4Xqs8gJXWTSZOtjQRSC7zdnIlv0G23MeOUhtUIq+PfXd0uPKpVFaKVNh9Rf34j/NwJ9ZVnkKyNBNJ5DIsOtGyvS5fhvaKBhoF0+ohmZYjW+WV/vmV7eExgC0cKrUggncfwBusibcg41WjqeuHg5ArJZtSZFH47eNLSHts3uIWjhVYkkM5jSI8AfOqnqK8xmhp9UwsHF9BD3T59GBQZ9GKtdmQWUlRhXnbCWa9jRK8u5zlDaEEC6Tyc9U5c3k8dTfXT3lwNqxFWJaiXul1TBkWZ2tUiWvT9TnUU5CU9AvFylVnTrJEEUitcFacG0qq0POm2E2ZuvuAbpbbz92lXi2iWsc7E8l05lvbEhDANqxEtkUBqhTF9gi0rS1bU1PG/PXKVJOoFx6rb+Xu1q0M0a116AafKzfPXueiduLq/BJK1kkBqBTdnPZMSwy3tRduyNKxGWJWQBoGUl6ZdHaJZi7Ydt2yP6dNFFuWzYhJIrTR9cKRle9Ph0xw+ab+rhIo2aHiFlLdHuzpEk/JKqvipQY/GdQ1+joX1kUBqpYFRfvQOUVcJ/XjDUe2KEdYjLFHdPnkAqmQ2D2vy+eZMjPVTfoX7unG5DPe2ahJIraTT6Zg1TB3mu2jbcYorazWsSFiFwF7gemY2DwWyt2tajlBVG+v4Yos68vGWS7th0MuvPGsm/zttMHVgBH71/c8VNY2/2YWDcnKCyMFq+3iKdrWIRhZvO0F+aTVgHsxwQ3LUec4QWpNAagN3Fz03D+lqab+39jAVNTIE3OFFJKnbEkhWobbOxFtr0i3t6UmRBMlifFZPAqmN7hjRA3dnPQCnymv4dJPM8uzwIpPV7eNbZcYGK7AsNdsys7feSce9o6M1rki0hgRSGwV5uXLbUHXZgXd+OywPyjq6yAZXSBUFUHBQu1oE1cY6Xlt9yNKeOjCCqAAPDSsSrSWBdAHuHtUTN2fzP92p8hre+S1D44qEpjwCICRebR9Zq10tgoUbj5F5ugIwXx3dN0aujmyFBNIFCPJyZc6Inpb2u78fJruosoUzhN3rMUrdPvKbdnU4uKKKGl7/Rb13dPOQrvTs4tXCGcKaSCBdoLljoi03SatqTfzrf/s1rkhoqudodfvI72AyaVeLA3ttdbrlcQwvVwN/GNfrPGcIayKBdIG8XA388crelva3qdlsyJAF/BxW16GgMw92oaoIcndpWo4j2nOimAUb1IUS77ssmkAZWWdTJJAuwvVJUfSP8LG0n1q6h6raOg0rEppx84GIBs8jHf5Vu1ockLHOxF+W7KZ+Uga6BXpwx/AeLZ8krI4E0kXQO+n459QEnHTm9uGCct76Nb3lk4T9ih6rbh/4n3Z1OKCPNx5j94liS/v5qfG41T+eIWyHBNJFio/05fZh3S3tN9dksDOrSLN6hIb6jFe3j2+B8lPa1eJADp8s498/HbC0pw2KYHhMkIYViQslgdQOHr2yDxF+7gDUmRQe/jqVyhrpunM4YQPAu36tHcUEh1ZqWo4jqK0z8YevUqms7yr393DmqWtiz3OWsFYSSO3Ay9XAKzMS0Z3pujtZzvM/yuqhDkeng95Xqe2DK7SrxUH89+eD7DqudtX9c1oCAZ4uGlYkLoYEUju5pGcgd49Un01auOkYPzRYNlk4iN5Xq9vpv4CxRrta7Nymw6d4a436UPpNQ6IY3z9Uw4rExZJAakePXNmb2DB11N1ji3fJQn6OpudoMJi7b6kphcNrNC3HXuUUVzLv8+2WaQN7BHny14nSVWfrJJDakatBz1u3DMLb1QBAWbWR+z7bLnPdORJnd+g1Tm3vXaJdLXaq2ljH3E+3U1Bmvvp00Tvx6o0D8HAxaFyZuFgSSO2se5AnL12fYGnvzy3lka9TMZlkBmiH0f86dXv/D1BbpV0tdkZRFP727d5GI1n/b0ocCZF+mtUk2o8EUgcY3z+Mu0aqD+X9tDePV1bJDNAOo9dV4Oxp3q4ugfSfta3Hjny0/ihfpWRZ2jdf0pUbkru2cIawJRJIHeTxq/sxpk8XS/uNX9P5ZttxDSsSncbFA/o0GNwg3XbtYsXuHP7vhzRLe1BXP56eJPeN7IkEUgfRO+l47aaB9ApWZxp+bPEuftmfp2FVotP0n6ZuH1gB1TK45WKkHD3NQ1+lWgYxRPi58/atg3E1yGwM9kQCqQP5uDnzwe3JBHmZn4uoMync99l2th0r1Lgy0eFixoGrr3m7tgLSlmlbjw3bn1vCnE9SqDGaZ1D3cTOwYHYywT5uGlcm2psEUgfrGujBgtlD8KofeVdVa2L2R1vY02DeLWGHDK4Q32BwQ+pn2tViw9Lzy7j1/c0UVZiXlHDRO/HezCR6hXhrXJnoCBJInaB/hC/v3jYYF735n7ukysitH2wmLbtE48pEhxpwq7p9bD2ckpWF2+JoQTk3v7fJMrzbSQf/uWEAl/QM1Lgy0VEkkDrJsJggXrtpIPr6qcGLKmq55f1N7MuRULJbEYOgS1+1nfq5drXYmKMF5dzy/mbyS6sB86xM/74+kWsSwjSuTHQkCaRONL5/KK/dqIZSYUUtN7yzkZSjpzWuTHQInQ4G3KK2d34BJpl093z255Zw/TsbOVFUadn3z6nxTBsUqWFVojNIIHWyaxLC+O8NAyxrKJ3pvvt1f762hYmOkXijupJsyQnI+EXbeqzc9sxCbnhnEyfrr4wA/m9yHDcOkWeNHIEEkgYmJYYz/9bBuBjM//xVtSbu+iSFb3ec0Lgy0e68ghs/k7TlPe1qsXK/HzrJre9vprjSPIDBSQcvTU/gtqHdtS1MdBoJJI1cFRfKgtnJltF3RpPCH75K5f3fD6MoMs2QXUm+U90+tBIKj2pWirX6Yksmsz7aSkX9OmLOeh1v3TKI65OiNK5MdCYJJA0Niw7iy7svJbDB+i3P/bCPP3+zi2qj3GuwGz3GQGCv+oYCWz/QsBjrYjIpPP/jPv6yZDd19fM9ujvr+eD2ZMb3lwEMjkYCSWP9I3xZNHeoZcVZgEXbjnPze5sb9aMLG+bkBMlz1PaOhVBb2fzxDqK82si9n23j3bWHLfuCvV35+p6hjOrdpYUzhb2SQLICPbt4sWzecJK7+1v2bTtWyOQ31skDtPZiwE3qhKuVhbBnsbb1aCw9v5TJb67np73qVFr9wnxYNm848ZG+GlYmtCSBZCWCvFz5dM4lzEhSh7ZmF1dx3fwNfLElU+4r2To3X0i8QW1vfgcc9P/0+53ZXPvGetLz1fn9xvYNZtHcoYT5urdwprB3EkhWxNWg58XrEvjbxFjLsPBqo4m/LNnNA1/soLSqVtsCxcVJvkvdzt0FR37TrhYNVBvrePb7vTzwxQ7L4AWA+y+L5r2ZSZYBPsJxSSBZGZ1Oxx0jerBg9hACGgx2WL4rh2teW8eu40XaFScuTkgsRF+utte/pl0tnexgXilT3tzAR+uPWvZ5uxl4f2YSf7qqr+VhceHYJJCs1KjeXVjx0Egu7Rlg2Zd5uoLr5m9g/poMjHUmDasTF2z4g+p2xmrI3aNdLZ3AZFL4cN0RJr6+rtE0WbFhPvzwwEjGxYZoWJ2wNhJIVizEx43P5lzKw+N6W7rwausUXvzffq57eyOH8kq1LVC0XY/REJaotjfY71VSTnElMz/cwt+Xp1mWjgC47dJuLLlvGF0DPTSsTlgjCSQrp3fS8dC4Xnx+16WENlj/ZWdWEde8to43f02XqyVbotPBsAZXSXsWQ7F9rSRsMiks3HSMK19Zy7r0Asv+IC9XPpqVzP9N6Y+bsyysJ84lgWQjLu0ZyE9/GMV1DSaYrKkz8dJPB5g2f4MMD7clsVPAr35uNpMRNs3XtJz2dCivlBnvbOSv3+6htNpo2X9lbAg//WEkl/UN1rA6Ye0kkGyIr4czL89I5MNZSYT4uFr27zpezKQ31vHUt7spqqjRsELRKnoDXHq/2k75CMpPaVdPO6ioMfLyygNMeO13UhqsiOzjZuBf0xN457bBBHq5tvAKQkgg2aSxfUNY+fBorh+sXi0pCny6KZPL/r2GL7dkYjI55jMuNmPQbeBRv9BcbTlsfEPbei6QoigsSz3B5S//xuu/pFNbp37fXZMQxs+PjmZGUhQ6nYyiE+cngWSjfN2deen6RD698xKiu3ha9hdW1PL4kt1MfWs9mw7b9rtuu+biCUPnqe0t70KFba2LtedEMTPe2chDX6aSU1xl2R/m68YHtyfx5s2DCPZ2a+EVhGhMAsnGjegVxIqHRvGXq/vi4aLeKN55vJgb393E7I+2yKq01mrIXeBeP11UTRlsfFPbelrpaEE5D36xg4mvr2PrUbV7zkXvxL1joln1yGgu7yfDuUXbSSDZAReDE/eMjuaXR8cwKTG80ed+PXCSCa/9ziNfpZJ1ukKjCkWTXL1haIN7SZvfseqrpNziKv6yZDeXv/Ib3+3MbvS5K2JDWPXIKB4b31dmXBAXTNcRc6TpdLq9sbGxsXv37m331xbnl3L0NP9csZ9tDW4ug3mNmemDI7l3dIw8A2Itqkrgv/FQVWRuj/ozjH1S05LOlldSxXtrD7Nw0zGqjY0fMegT4s1TE/sxspfMzi3M4uLiSEtLS1MUJa6t50og2SlFUfh5Xz7/+t9+DjWYxBLMzzZNTgznvsuiiQn21qhCYfHbv+DXf5i3XX3goZ3gEdDyOZ3g2Kly3v7tMIu3HafmrGfdogLceeSK3lybGCHT/ohGJJBEs+pMCou3H+e/qw6S3eDGM5if0RwfF8qckT0Z1NVPRkJppaq4/iqp/lmyYQ/Alc9pVk5adglv/5bB8l3ZnD1Ys4u3Kw9e3osbkqJwMUiPvziXBJI4rxqjiaU7jvPWmgyOnTr3XlJipC+zh/dgQnyY/KLRwu+vwOpnzdsGN3hgO/hGdNqXN9aZWJWWx4INR9l85Nz7WKE+btw1qic3D+mKu4vMsiCaJ4EkWs1YZ+KH3Tm88Uv6OV15YH4HfOsl3bhxSBQhPjJkt9PUVMBrA6Es19wedDtc2/Hz3J0qq+bLrVl8uulYo6HbZ3QP9GDu6GimDorA1SBBJM5PAkm0mcmksGpfHh+uO9LkO2InHVzWJ5gZyVGM7RuMs16umjrc1g/gh0fM2zo93L8Zgnq1+5epMyn8fugk32w7zsq0vEYTn57RP8KHe0ZFMyE+TO4RiTaRQBIXZW92MQvWH2XZzuwmfzkFebly3aAIrhscSe8QGQTRYepq4Y1kKDxibsdOgRkft9vLp+eXsXj7cZZsP05eSfU5nzc46ZgQH8btw7rLPUVxwSSQRLsoKKvmi82ZfLk1ixNFlU0e0yfEm2sHhDMxIYxugZ5NHiMuwu5vYPGdavvuNRA+8IJfLruokhV7cvl+ZzapWUVNHtPF25Wbh3Tllku6EizdtOIiSSCJdmUyKazPKOCrrVms3Jt3zpDfMxIjfZkQH8YVsSH07OLVyVXaKZMJ3hkFebvN7e4j4fbvzUMiW+lEUSUrdufww+4cdmQWNXmMwUnH2L7BTB8cyZg+wTKQRbQbCSTRYQrLa/g29QRLtp9gdwtLXPTs4skV/UIYFxvCoK7+ct/hYhz6GT67Tm3f+Dn0vabZw+tMCqlZhfy6/yS/7M8nrYWpovqGenN9UhRTBoTL7NuiQ0ggiU5xpKCc5Tuz+W5ndpMj9M7w83BmeHQQw2OCGBETJLNCtJWiwKfTIOMXczugJ9y3GQwulkNyi6vYeLiA3w6c5LeDJymsqG325boFejAhPoyJCWHEhft2dPXCwUkgiU53ILeU5buyWZWWx/7clpdS7xrgwfCYIIb08CepWwCR/u5yw/x88tLg7eGgmLtLS0c/yy/+17Pp8Ck2ZpziaBPPkjXUI8iTCfGhTIgPIzbMR/69RaeRQBKayjpdwep9eazen8/GjFMYz7MWU7C3K0nd/RncLYABUb70C/PBw0Um5DyjqraOPSeK8Vz1J/qd+AaAYsWD0dX/oYimRzkanHQkdw9gbN9gLusbTHQXTwkhoQkJJGE1Sqtq2XT4NOvTC1iXXkB6C117ZzjpILqLF/ERvvSP8KVvmDcxXbzo4u1q179UFUWhoKyGg3ml7M8t5WBuKXuyi9mfW0qdSSGAEta4PoyPzjzi8SPjVTxrvN1yftcADy7tGcBlfYIZ0SsIbzdnrf4qQlhIIAmrlVtcxbr0AjYfPsW2Y4UcLihv9bnebgZigr2I6eJFdP2fXQM9iPBzx9OGljgorzaSVVhB5qkKMk+bPw7mlXIwr4zT5S0vOX+3/nuecP4CACNOvBzzMTH9BnFJzwAi/eXenLA+FxNItvNTLWxSqK8b0wdHMr1+ufWCsmq2HStk27FCth8rZG92CZW1dU2eW1plZEdmUZNDl33dnYnwcyfcz51If3fCfN0I8HQh0MsFfw8XAj1d8fd0xsvV0CFXWYqiUF5TR0llLUUVtZwsqya/pIqTZdWcLK0mv7Sa7KJKsk5XUFDWcug0xd1ZT3ykL/qIeynftw7P8iwMmHiMT2DQpDYNAxfCVkggiU4V5OXKVXGhXBUXCpiHLB8pKGPPiRL2nChmT3Yx6fll5/0lXlxZS3FlbYtDnMG8iqm3mwF3Fz0eLnrcnfX12wbcnJ3Q6XQ46XToMHcdOtX/oq+pM1FjNFFtNP9ZU2ei2lhHaZWRkspaSqqM1J3nXlnr/01c6B3iTZ9Qb/qEeJMY5UevYC8MZ6Zr6vkP+HqmeTt9FRz4scVh4ELYKgkkoSm9k46YYG9igr2ZMlCd3bqooob0/DL142QZh0+Wk1NcSW1d64Ogps7EqfIaaH1PYYdw0TsRGeBO1wAPovw9iO7iSe9Qb3qHeBN0vueB+l0LPUbDkd/M7RWPQ8/LwEW67IR9kUASVsnPw4Wk7gEkdW+8UJ3JpHCyrJoTRZWcKKwku6iSE0WV5BZXUVhRw+ly80dRZS0dcHu0Se7Oerp4u9LF25XgBn+G+LjRNcCDroEehHi74XShDwvrdDDh3zB/GJhqoTgT1r0CY59q37+IEBqTQBI2xclJR4iPGyE+bgzq6t/scXUmhaKKGgoraiitMlJZU0dFTR0VtXVU1hipqKmjqtaEgoKimO8JKQqYFFBQcNY74Wowf7jUf7ga9Hi7GfBxc8bH3RkfNwPebs6dM+1Ol94w9H5Y/19ze/2rkHgTBEZ3/NcWopNIIAm7pHfSEejlal/T44z6E+xeBCUnoK4GVvwZbvlGBjgIuyEzKgphK1y94Krn1Xb6z7D/B+3qEaKdSSAJYUtiJ5sHNJzxv8fNq80KYQckkISwJTodTHgJnOpnZSjOgrUvaVuTEO1EAkkIWxPUC4bNU9sbXoPcPdrVI0Q7kUASwhaN+hP4dTVvm4zw/UNganrGCyFshQSSELbIxRMm/kdtn0iBre9rV48Q7UACSQhbFTMO4meo7dV/h+Lj2tUjxEWSQBLCll31PLjXPyBcUwY//JFOm6JCiHYmgSSELfPqAlf+Q20fXAFpy7SrR4iLIIEkhK0bcDP0GKW2V/wZKos0K0eICyWBJISt0+lg4n/B4GZul+XByic1LUmICyGBJIQ9CIyG0X9W2zs+NU8tJIQNkUASwl4MexDCEtX2dw9CVbF29QjRRhJIQtgLvTNMma9OK1RyAlb+VduahGgDCSQh7ElIHIx+TG1v/xjSV2tXjxBtIIEkhL0Z8QcITVDb3z0IVSWalSNEa0kgCWFvzum6Ow6rpOtOWD8JJCHsUWj/xqPuti2AjF81K0eI1pBAEsJejXgYQuPV9ncPSNedsGoSSELYK0vXncHcLs4yrzArhJWSQBLCnoXGNx51l/oZpH2nXT1CtEACSQh7N+IRiEhS298/BKV52tUjRDMkkISwd3oDTHsXnD3M7crTsOx+WaZCWB0JJCEcQWA0XNVgmYr0VZDyoXb1CNEEgxZfVFEUFHl35vB0Oh06nU7rMhzH4NlwYAUcWmlur3wKeoyGoBht6xKiXqcFkqIolJaWkp+fT21tbWd9WWHFdDodXl5ehIeH4+QkF+sdTqeDa9+Aty41d9vVVsDSu+GOleZuPSE01mm/BXJzczlx4oSEkbA48yYlOztb61Ich3cITHpVbZ/YBr+/rF09QjTQKW+LFEWhpMT8QF5AQACBgYHyjlhYwqisrAxFUaT7rrPEXgsDbjEPAQf47UWIGQeRg7WtSzi8Tgskk8kEQJcuXSSMBADe3t6Aek9RAqkTjX8Bjv4ORZmg1MGSOXDPWnD11roy4cAkGYRwRG4+MPUdoP5NwOnD8OOfWzxFiI7msIG0Zs0a4uLiiImJYc6cOdTV1XXK112xYgUDBgzAYDDw6aefXvTrGQyGJreTkpKaOlwIVbdhMPJRtb3zc9j9jXb1CIfnkIFkMpmYM2cOixYtIj09nZKSkvOGQ1lZGTU1NRf09WpqaigrKwOgV69efPrpp9x8880X9FqtlZKS0qGvL+zEmMchcojaXv4wFB7VrBzh2DQLJEVRKK6s7bCPlp5z2rp1K+Hh4cTGxgJw5513snjx4iZrXLNmDbNmzaJ///4UFhYCcOzYMSZOnEhSUhJJSUn89ttvTX6dlJQU5s2bR+/evTl06BAAMTEx9O/fv0330UwmE08++SRDhgwhISGBJ5544rznnH3l9MQTTxAfH8+gQYPYvXs3AOvWrWPw4MEMGDCA/v37s3nzZgB27drF2LFjGTx4MCNGjLAcL+yQ3hmuex9cfczt6hJYPAfqZDSs6HyaPXxQUmUk8dmVHfb6O5++El935yY/d/z4caKioiztrl27kpWVZWkfOXKEjz/+mC+//JJ+/foxc+ZM3n33XVxcXAC44447eO2114iLiyMzM5MxY8aQkZGBTqcjNzeXTz/9lIULFxIcHMzMmTN58cUX8fT0vOC/y4IFCwDYsmULJpOJKVOmsGLFCq6++upWnV9XV0dkZCS7d+9m2bJlzJ49m5SUFP71r3/xxhtvMHToUIxGI1VVVdTW1nL33XezePFiIiIi2Lp1K3PmzLGElbBD/t1g4n9g8Z3m9vGtsOafcPnftK1LOByHfBqupaunxYsXM2PGDB566CE2bNhAQEBAo8+XlZWxbt06brnlFsu+mpoa8vPzOXbsGMOHD2fGjBn8+OOPREREtEu9P/74Izt37uSHH34AoLy8nEOHDrU6kABuv/12ACZPnszs2bMpLy9n9OjRPPjgg8yYMYNJkybRt29f9uzZw969e7nmmmss554+fbpd/h7CisVPh4xf1KHgv78CPcdAj1GaliUci0MGUlRUVKMroszMTCIjIwG44oorePPNN/nkk0+YNGkSt912GzfccAP+/v6AufvMw8OD1NTUc17X19eXjz76iI8//phJkyZxyy23cPPNNxMWFnZR9SqKwksvvcSUKVMu6nXO9uijjzJx4kRWrlzJtGnTePrpp4mNjSU6OrrJv5+wc1f/CzI3wekMQIEl98C968Ej4LynCtEeNAskHzcDO5++skNfvzlJSUkcP36ctLQ0YmNj+eCDD5g2bZr5PB8f5s6dy9y5czlw4AAff/wxycnJJCYm8uGHH+Lr60tcXBwffvghd9xxBwDbt29n0KBBuLm5ceutt3LrrbeSlZXFwoULGTduHJGRkbz77rt069atxZpnzpzJ1KlTmTp1aqP9V199NfPnz2f8+PG4ubmRnZ2Nk5MToaGhrf73WLhwIXPnzmX58uX07NkTT09PDh06RJ8+fejTpw+lpaWkpKQwbdo0SktLWb16NZdffjmKopCamsrAgQNb/bWEjXL1gukfwvvjwFQLpdmwbB7c+Jl52iEhOtqZhxLb8wPYGxsbq5xRV1enpKWlKWlpaUpdXZ1iDVavXq3069dP6dmzpzJ79myltra22WPr6uqUlStXKkVFRYqiKMrRo0eViRMnKgkJCUq/fv2U2bNnt/i1NmzYoGRmZiqKoig///yzEhERoXh4eCj+/v5KRESEUlVVpSiKosTFxSlbt24953yTyaQ8/fTTSlxcnNK/f3/lkksuUfbt26coiqLo9XrLcS1tP/HEE0p8fLwycOBAZefOnYqiKMp9992nxMbGKgMGDFAuu+wyS427du1SxowZY/n7Pfnkky3+/S6UNX5fCEVR1r+mKE/7qB8b52tdkbAhsbGxCrBXuYDs0CkdMOu2TqfbGxsbG7t3717A3M114MABAPr06SMzNTShqKiIGTNmsHJl+w/0MBgMGI3Gdn/diyXfF1bKZILPpkPGanPbyRnu/AkiZGohcX5xcXGkpaWlKYoS19Zz5TeAlfDz8+uQMBKizZyczLM4eNV3CZtqYdFsqCzStCxh/ySQHIA1Xh0JK+fVBaZ/ALr6XxFFx+C7ebLKrOhQEkhCiKZ1HwGXNXgIe9/3sOVd7eoRdk8CSQjRvBGPQs/L1PZPT8KJ7drVI+yaBJIQonlOTjDtvbPuJ82S+0miQ0ggCSFa5tXFPN9do/tJD8j9JNHuHDaQbGn5iQULFjBnzpwOruxca9asYdy4cR1+jrABPUbCmIb3k76DLe9pV4+wSw4ZSI6w/MTFkpF54hwjH2l8P2ml3E8S7Uu7QFIUcz90R33Y0fITACdPnmTixIn07t2bmTNnAublI8aOHWs5Jisri5iYGBRFYdasWcydO5chQ4bQq1cvXnnlFctxa9asYcSIEQwePJgrr7zSMq/frFmzuOeeexg6dKhlMtaKigqmTp1KbGwsEydOpKioCDDPmD5+/HgSEhJITk5m06ZN59SclZXF6NGjGTRoEPHx8e2yIKHQkJMepr0LXiHmdl0NfH07VMjku6J9aDe5alUxvNjy3G4X5bFj4O7X5KdsbfkJMIfbrl278PX1JSkpifXr1zNixAhOnjxJRkYG0dHRfPjhh8yaNQtd/bxje/bsYd26dVRUVDB48GCuuOIKIiIieOKJJ1ixYgW+vr4sWrSIP/7xj3z11VcApKens3btWpydnVmzZg1btmxh+/btJCQk8PDDD/Pcc8/x73//mwceeIAJEybw4IMPsnXrVq6//npL6J4RFBTEihUr8PDwoKSkhMGDBzNx4kT8/Jr+fxE2wCsYrvsAPrkWFBMUZ8KSu+Hmr80DIIS4CA4523dL0yVZ4/ITAGPHjiUwMBCAgQMHcuTIEYYPH85dd93F+++/zz/+8Q8++eQT1q5daznnpptuwsXFBRcXF6699lrWrl1Lt27dOHDgAKNHjwbM3Zfe3t6Wc2bMmIGzs7qO1MCBA0lISADMV1Bn7mX99ttvLFy4EIDk5GT8/PzOCSSj0chDDz3E1q1bcXJyIicnh/T0dFle3db1GAmXPw0/P21up6+CtS/BmMe0rUvYPIcMJFtbfgLA1dXVsq3X6y33eGbOnGlZ2TU2NrZRCOqamKFZURRGjhzJt99+2+TXOftKrqnXaOpzTR33yiuv4O7uTmpqKnq9nsGDB1NVVdXs6wkbMvwh80J++5eb22v+CZGDIUYGtIgLp10gufmau9U68vWbYWvLT7TEz8+PUaNGcffdd/PWW281+tyXX37JXXfdRUVFBd9//z1Lly4lLCyMe+65hz179tC/f39qa2vZv38/8fHxTb7+jh07LMcuXLjQcmU1ZswYPvroI+bNm8e2bdsoKiqiV69e5OfnW84tLi4mNDQUvV7P5s2b2blzZ6v/XsLK6XQw5S14Nw1OHwYU89Ln96wFv65aVydslHadvjqd+R5PR3208M5er9fz/vvvM336dKKjo/Hy8uK2224757g+ffrw/PPPc/DgQebOnWvZ/9lnn7F06VISExOJjY3ljTfeOOfcqKgonnjiCfbu3cszzzxjGcSwevVqIiMjWbRoEQ8++CCRkZFUV1cD5mBreG+rtWbOnImiKI1WeQWIjY1l+PDhJCUlcd999xEfH09QUBBffPEFc+bMITExkQEDBjTq5jtbcnIyf/vb34iNjWX//v089dRTALz22mssX76chIQE5s6dy+eff96oqw9g3rx5LF26lISEBF599VWSk5Pb/HcTVszNF2YsBIO7uV1ZaB7kYKzWti5hs2T5CStxMctPPPvss9TV1fH3v//dsm/WrFmMGzeOW2+9tT3LbFfyfWEnUr+Ab9U3bCTdARP/o109QlMXs/yEQ95DskYXuvzE8OHDqaysZPXq1R1QlRCtMOAmOL4FUj40t1M+hMgh5v1CtIEEko1bv359k/sXLFjQuYUIxzb+BcjeYf4AWP4whMZDaH9t6xI2RfpIhBAXz+AKMz4Bd/NoVIyV8PVt5ucNhWglCSQhRPvw6wrT3gfqBxSdPgxL7zUviS5EK1h1l11BZQFLDi0hJTeFcmM5ngZPkkOTmdprKkHuQVqXJ4Q4W69xMOZx83NJAAd+gN9fhtF/0rYuYROsMpCqjFW8sOUFlmUsw2hqPMnnxpyNvLXzLabETOHxIY/jqndt5lWEEJoY9Wc4nmKewQHg139AWAL0vkrbuoTVs7ouuypjFff+fC+LDy0+J4zOMJqMfHPwG+aumkuV8cKe/NfpdAwYMMDy8euvv15M2UKIM5yc4Lr3wL9H/Q4FFt8FpzI0LUtYP6sLpBe2vEBKXkqrjk3JS+HFrS9e0NfR6/WkpqZaPi677LLznySEaB13f7jxc3Cun4qquhi+vBmqS7WtS1g1qwqkgsoClmUsa9M536Z/S0FlQbt8/eLiYm644Qbi4+NJTEzk+++/t3zu66+/ZsCAASQmJjJy5EgAnnnmGZ577jnLMePGjWPNmjUA/O1vfyMuLo6EhASuuOKKdqlPCJsSEmueXuiMk/th6VwZ5CCaZVX3kJYcWtJsN11zjCYjSw8t5a6Eu9p0Xl1dHQMGDADAYDCQkpLCM888Q3h4OF999RVHjx5l6NCh7Nq1i4KCAh577DE2bNhAWFgYp06davG1T58+zTfffMOePXtwcnKyrKMkhMOJmwI5j8C6+vW49i+HdS/DKBnkIM5lVVdIKbmt66o729bcrW0+p2GXXUqK+euuWbPGsrxC9+7dueSSS9i6dSurV69m2rRpllm7zywD0RxfX188PT2ZNWsWn332GXq9vs31CWE3xj7VeBbwX/4BB9s+K4mwf1YVSOXG8k4972xnL6Fwpt3cfH8GgwFTg+6HM0sr6PV6NmzYwO23387OnTtJTEykuFgeEBQOykkP171/1iCHOTLIQZzDqgLJ03Bhq6pe6HlnO7OkApjXSNqyZQtDhgxh3LhxLFmyhJycHABLl12PHj3Yvn07ABkZGezYYZ42pbS0lFOnTnH55Zfzwgsv4ObmxvHjx9ulRiFskgxyEK1gVYGUFHphK4kmh7bPsgZPP/00WVlZxMfHM2nSJN5++22CgoLo168fL7zwAuPHjycxMdGydtJ1111HTU0NsbGxPPnkk5Z7UsXFxUyePJmEhAQSEhKYPHkycXFtnvhWCPsSEgtT3lTbJ/fDt/dBB6w4IGyTVS0/UVBZwBXfXNGmgQ0GJwOrpq+SmRtskCw/4aB+fgbWNVieYuxfYdQfNStHtK+LWX7Cqn4DBLkHMTl6cpvOmRIzRcJICFsy9q8Qfbna/uU5OLRKu3qE1bCqQAJ4fMjjJIW0rusuKSSJx4c83sEVCSHalWWQQ/f6HQp8cycUpGtZlbACVhdIbgY35o+bz/Te0zE4Nf2YlMHJwPTe03n7irdlLjshbJFHwLmDHL64ESqLNC1LaKtTHozV6XTodDoURcFoNOLi4tLi8W4GN54e+jT3D7ifpYeWsjV3q8z2bYcaDpk/e8i9cAAhcTB1Pnw909w+dQgW3wk3f22+ihIOp9MCycXFherqak6cOEFERAQGw/m/dIBrAHf2v5M7+995zudMMv2IzTszfN7Z2VkCyVHFTobRj8NvL5jb6T/Dz0/Dlc+1fJ6wS502dVB4eDiZmZlUVVWRkSEPxAlVcHCw1iUILY1+DPL3wr76uSM3vA7BcTDgJm3rEp2u0+4hubm50bVrV1xdXeXdsADAyckJPz8/fHx8tC5FaMnJCaa8DSH91X3fPwhZbZ8STNi2Tp1c1c3NjZ49e6IoSrPT8QjHcebeohC4esFNX8C7l0FFAdTVwFe3wF2/gm+E1tWJTqLJbN/yi0gIcQ6/rnDDQvh4EpiMUJZnnl7ojv+Bs7vW1YlOYHXDvoUQDqzbMLjmZbWdkwrL5sn0Qg6io6YOKnF1dfWOjo5u99cWQjiAsjyobLCOmGcX8Gh52RdhHTIyMqiuri5VFKXNN4c7KpByAQ8gq91fXAghhDWLAioURQlt64kdEkhCCCFEW8k9JCGEEFZBAkkIIYRVkEASQghhFSSQhBBCWAUJJCGEEFZBAkkIIYRVkEASQghhFSSQhBBCWAUJJCGEEFZBAkkIIYRVkEASQghhFSSQhBBCWAUJJCGEEFZBAkkIIYRVkEASQghhFSSQhBBCWAUJJCGEEFZBAkkIIYRVkEASQghhFSSQhBBCWIX/Bx0prQHncM+/AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "e = np.array([0.0,0.6,1.0,1.5])\n", "# Polar coordinates\n", "t = np.arange(-np.pi,np.pi,0.01)\n", "theta = np.stack((t,) * 4, axis=-1)\n", "r = (1-e**2)/(1+e*np.cos(theta))\n", "\n", "#Cartesian coordinates\n", "x = r.T*np.cos(theta.T)\n", "y = r.T*np.sin(theta.T)\n", "\n", "fig,ax = plt.subplots(nrows=1,ncols=1,figsize=[3.3,3], dpi=150)\n", "#ax.plot(x[0],y[0], label='e=0, circular')\n", "ax.plot(x[1],y[1], label='0-np.pi/1.4) & (t-np.pi/1.4) & (t-np.pi/1.4) & (t-np.pi/1.4) & (t1$).\n", "\n", "As such, we've arrived at a generalisation of Kepler's First Law.\n", "\n", "We can use the above expressions to determine what the velocity of the body throughout its orbit is by realising that\n", "$$\n", " \\frac{L^2}{m^2} = GMa(1-e^2)\n", "$$\n", "Knowing that $L=m v_t r$, we get\n", "$$\n", " r^2 v_t^2 = GMa(1-e^2)\n", "$$\n", "Now, typically the velocity of the body is given by $\\textbf{v} = v_r \\hat{\\textbf{r}} + v_t \\hat{ \\pmb{\\theta}}$ - that is, there's both a radial and a tangential component. However, at perihelion, the velocity is entirely tangential, and we know that the distance between the planet and the Sun is $a(1-e)$. As such, at perihelion, the velocity is given by\n", "$$\n", " v_{pe} = \\left[\\frac{GM}{a}\\frac{1+e}{1-e}\\right]^{1/2}\n", "$$\n", "A similar analysis at aphelion provides us with\n", "$$\n", " v_{ap} = \\left[\\frac{GM}{a}\\frac{1-e}{1+e}\\right]^{1/2}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "5d892669-c109-4c35-ae78-204be565ddd2", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.10" } }, "nbformat": 4, "nbformat_minor": 5 }