From 428d0c1670317dd373ee345c76b015c6175f08fa Mon Sep 17 00:00:00 2001 From: Grant Curell Date: Mon, 29 Apr 2024 11:47:13 -0400 Subject: [PATCH] Updated the gradient_based_optimization tutorial with my understanding of what is going on with an explanation at a level that an engineer with a bachelors in comp sci understands - Added/updated comments for the gradient_based_optimization tutorial - Shamelessly added my name to the bottom of that tutorial so that I have something to point to for the time spend with leadership. Feel free to delete it if my explanation is terrible :-p - Added comments to initl.py Signed-off-by: Grant Curell grant_curell@dell.com --- .gitignore | 2 + .../gradient_based_optimization.ipynb | 138 +++++++++++--- doc/notebooks/tle_object.ipynb | 173 +++++++++++------- dsgp4/initl.py | 157 ++++++++++------ 4 files changed, 320 insertions(+), 150 deletions(-) diff --git a/.gitignore b/.gitignore index f192351..8dfcabd 100644 --- a/.gitignore +++ b/.gitignore @@ -139,3 +139,5 @@ dmypy.json cython_debug/ .vscode + +.idea/ diff --git a/doc/notebooks/gradient_based_optimization.ipynb b/doc/notebooks/gradient_based_optimization.ipynb index 207c070..08b874a 100644 --- a/doc/notebooks/gradient_based_optimization.ipynb +++ b/doc/notebooks/gradient_based_optimization.ipynb @@ -1,24 +1,47 @@ { "cells": [ { + "attachments": { + "1f7a6234-cf18-4c3d-a0e4-cc66a5d5b53a.png": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAIAAAC6s0uzAABxf0lEQVR4XuydB3gUVReGN9KLIEiXjnTBH0Gk2EFABEEREVBREWlKEcHQpfdeBJEiJfQeIBAgtFBDCZ2EkhBaKiG9bHb/sxky3Nwt2TLt3j3vkwd250w2u9/M+b65s1N0RgRBEARBFEdHT0AQBEEQRH4wgBEEQRBEBTCAEQRBEEQFMIARBEEQRAUwgBEEQRBEBTCAEQRBEEQFMIARBEEQRAUwgBEEQRBEBTCAEQRBEEQFMIARBEEQRAUwgBEEQRBEBTCAEQRBEEQFMIARBEEQRAUwgBEEQRBEBTCAEQRBEEQFMIARBEEQRAUwgBEEQRBEBTCAEQRBEEQFMIARBEEQRAUwgBEEQRBEBTCAEQRBEEQFMIARBEEQRAUwgBEEQRBEBTCAEQRBEEQFMIARBEEQRAUwgBEEQRBEBTCAEQRBEEQFMIARBEEQRAUwgBEEQRBEBTCAEQRBEEQFMIARBEEQRAUwgBEEQRBEBTCAEQRBEEQFMIARBEEQRAUwgBHJ6NGjR6VKlcSnOp1u7Nix4lNEMfz8/EB8+JcuZCHMsHnzZrqAKIXki4DqPoQJMIDdgrt37/bv37969eoFMqldu3a/fv0CAwPp+VxDvgBet27dnDlz6KnZgT+ty8TDw6No0aJvvPFGr169Tp8+Tc+nMR4+fAgqXbx4kS5ksXHjRvhQ27ZtIyfWr18fJh4+fJicWKFChaZNmxrNAthcPUfdf+XKlTB/vnz5Hjx4QE7/4IMP6tatS05xncTERBDExtaDE8C64cSqKKgErFmzhio1a9YMptv/2V1fBDmCAcwiGMD8s3v37oIFCxYpUqRv375Lliz5559/fv/998qVK0NQhYSE0HO7AGUBycnJ6enpL8ou8Nlnn+VoLjDD//73vzWZLF68+LfffitTpgx43ODBg+lZtcS5c+fgTULC0YUsIKFhBlhk4pRnz5699NJLuXPnnjBhgjjx/v37MNvQoUPhcUZGBogP/wolc/UcdX8hgIFff/2VnC5HAEdGRuqk23QTcCWA8+fP/+mnn5LT7927J0y3/7O7vghyBAOYRTCAOef27duFChWCIe+jR4/I6RCN8+bNA9cmJ4okJCTQk+xAPgsw9y9zYAaYjZySlJTUsWNHsDnIY3K6psgxgIEqVao0btxYfOrj4wMbT127dm3durU40cvLC15n586d4hQRc/UcdX8hgGH7BgbBsEEgTneHAP7yyy9hWwfelTh90qRJpUuXfvfdd+3/7K4vghyRr/sQ+cAA5pxffvkF+jzHPbHQvZDTkNawsV+4cOEOHTrAxGPHjn311VcVKlTImzdv+fLlBw0aBJFG/tb27dvBg8CU4d9t27bZ3gX94MGDH3/8sVSpUvBqderUWb58uVgSzGjjxo0TJ0587bXX4AU//vjj4OBgoQouryOw5jLmAQzEx8cXL14cXtNgMAhTYFw4Z84ceAPwV+DNgD4xMTHi/BCHrVq1evXVV2F8U7lyZXjDYgl+ce7cuW+88Qb8YokSJSD8YGaxCsPut956C36rWLFiXbp0IbdshJS6du3ahx9+WKBAgXLlyk2bNk0oCR+cxGISf/fdd3ny5BHFHz16NLyN1atXFy1aVBzm9u/fH1I5KirKmH0XtEX1bAtujhDAmzZtgij67bffxOnmAWxNB9jag1H706dPhaczZ87UETsn9Ho9rHXDhg0TBpck4ip06NAhyLyCBQvCp/7888+vX78uTAdgHpgT3j+sgVAtUqTIDz/8kJiYKM5ABnBaWtpff/31+uuvw6eGdaN58+YHDhwQ5yQRVPrvv/+gNchtOPjIIIL9n93pRQCCCy8IK2T37t2p/f+2uw9hAgxgzgG7B6+hp5oB3QudXK1aNXiwZMkSMHeYCC7Ttm3byZMnL126tGfPnrly5YI8Fn9l//79YKmQBLNnzx45ciQYHxiBtQB+8uQJRDhk+fjx4//++28wUKiK34oJZtSgQYOGDRvCRPBH8FlxzAf+CGMvyDxh9zL4jjCdwmIAA/DO4cWvXr0qPP35558hRXr16gUf888//wRvffvtt8GUoRQeHg7WWaNGjRkzZixbtgw+VO3atcXXAU+H14ENFIhhyA/YRlmwYIFQAg+F8APPBZseN24cvFUIbzFswH9hKcBnHzhwIMwAPguvs3fvXmOmLCAIPIXtAOHT3blzR/yLIqC/jvhOF14B5oetJZgofpEPEonvlgxgi+rZFtwcIYBhg+Onn36CPBAHwVQI2dDhwoUL8Aq7d+8W5gT1YOVp1KiR8FTYDeDt7Z2QkACrBzz+4osvhDcsfEBfX19YarBopk+fLrwyLClIa+HXhQCGjwOjVfjTsIjhKcS5UDVmD+ARI0bAm4QVABbxrFmzunbtOnXqVHFOEkElGKR269btvffeEyZeunQJJp46dcr+z+7cIhA0h5UTqp6enrDpRq5UOXYfwgQYwDzz7Nkz6OGOHTuSE6GHI7MQB1WQuzAn9Dk5JzXenTJlClhMaGio8BQ8pWzZsrGxscJTcBld9uGpjghgSEGYWRifCXzzzTfgGsKfEMwI8iM1NVWowoAJply5ckV4ar4HzxxrAQz+pcvaN3v8+HF4vG7dOrHq4+MjTgFn1GXGjFgVOXz4MJQGDBhAThRG1SEhIbBpMmnSJHE6vG1IC3GKMAAStmkA+IxlypTp1KmT8NSeXdAweoZ5hG9809PTYaMBhmXwuHTp0osWLYIHcXFx8B4gVIT5yQA2WlIvR8EpxACG7QP4aKIOZAjZ1gFG6jAwFUIRdIMhXefOnWH++Ph4mAIpIo6PLe6ChpWtVKlS0dHRwlNIZZj/+++/F54KAQwbB+L8kN/wJ8SnJG+++abF9cQcMYBhywDWfGFEO3To0KpVqxod+exGxxcBbBHC54V8TU5OFqrwHqA6ZswY4WmO3YcwAQYwz4SFhUFbfvvtt+REMCBdFjDUEyYKASyGKwWMS8AWjx49CvPs2LEDpjx69EhnFth16tSxGMBguK+88goM2sTgBwRPP3HihDHLjGBwI/6uMGASv9E09y9zrAUwDHTgpdauXQuPITkg9SMiIsh3UrhwYRgzGbPeBrxnYUBMIuzgFQOABMIDSsHBweRrgrG2bNlSmAGcGv6EuA8c+Pzzz2HoIzy2J4CFxBK+8Q0ICNBl7m41ZsYMDM6MmeMhXebOUmF+OwPYhuAUYgDD4x9//BEGwcIhBWQI5ahDmzZtmjRpYszanjh//jyEqLD7Fz5I/fr1hdkizQJYWNnIES0AasCYUngsBPDZs2fFKrwZmAIboOIUEXjPMJQMCgqiC2aIAQzrA+gPcsGCqFChAow4jQ5+dkcXwcmTJ3Vmxy7UqlULhstG+7oPYQIMYJ6BDWSd2Qj49OnTvr6+EEi67AEMG+ziF4oCkMcwvVixYjoCweVPnToFj8nvcY2ZNmoxgMPDw8lXIBHOrhHMaMOGDeLvCt8Frlq1Snhq7l/mWAtgcgT86aefUm9AABLRmJlzMDCFpzBWgykrVqxISUkRXgTC47XXXiNfVqRv3770y2UiJgo4NVgn+SugKmSA8NieAAbat28vfOML4yQYGwkTYfEJssDACF5E3H1tZwDbEJyCDGByEEyGUI46wHBQ+CZ7yZIlMHqDKbAVMnr0aHhQsmRJ2MQRZjMPYIsr26BBg2CicLSgEMBPnjwRq8IbtniQP2xHwuYgVGF8+ccff9g4GU8MYHjcp08fGHQeOXIEpsAGhNHBz+7oIli/fj08PnTo0ItfMBqhkYVtDouCUN2HMAEGMOeA01n8DljodjKACxUqRM6g1+tr1KgBDT916lQY9UJmgzXosqLCHgvQZdno48ePdZkDcV8zIJuN2Z1OQHh7YiyZ+5c51gJY+A5YME0YNkF60W/C1/fSpUvi/PDRRowYAUMNXeaJnsI+UhsB3Lt3bxj9+Pj4UK8JryPMQDq1QA/ieBk7A3jatGkwG7zPr7/+Wtyi8vf3h4kPHjz4+OOPy5UrJ85sZwDbEJyCDGBj5tfhwiCY/Gg56iDsQYFQ6d69O3wKY+YOiQ8//PDGjRs6IoqcDuBI4kBl4Q2LXxJTREdHw9bVN998A0mcK1euZcuW0XNkQqp07NgxeAzv9s033xSqDn12RxcBBrCbgAHMOcIBKWfOnKGm5xjAFy9e1BF7NY1Z3zMJBmHPTjDRRiHLX3755a5du4olCttmBLRr1y5Hc7EYwMJR0BUqVBD2APfr1w8Ml/pu2xrr1q2D9yC4s41d0NOnT4fZbt26RReysB3Awi5la8knImTtwoULYTtAXGowQM+XL9/atWth2XXu3FmcmQpgc/VyFJyCCuDbt2/DIHjgwIHkR8tRh+Tk5Lx58/71118VK1YUjl/bsmVLgQIF4EPBL4qnyUVFRYlrjoDFXdCwSUTtgrY/gEVg9YBRuLVNK1IlWH/gbcNT8Qh2hz67o4vA4i7o2rVr4y5ozsAA5pygoKCCBQuCU5A76IyZ18bS2Qzgy5cv64h9kmBAEG+iQRjtOAyEtFEYM4H5Usf4RERECA9smxHQpUsXGKyIVYuYB7B4HvCSJUuEKcIuxOHDh5OzpaenC4f/xMTEkN/UCl9VQjwYbR6EBWkEod6tWzfyd+GxeMSZ7QAWxn85XucrNTUVBp1NmzaFmSGMxelNM4GJ8+bNEydSAWyuXo6CU1ABbMwaBNesWVP8aDnqADRv3hx+RZc5lDdmfTdRo0aNatWqifPAUoOJkO7iFGPmyla6dGnxGGBYkcwPwrIzgMn3A8CGixjkFJRKO3bsgD8kbiiQizXHz+7oIhAOwqpfv774JcjevXt1eBAWd2AA8w8YB4wzihYtCuO/pUuXCqffwKAQLGz9+vXCPOYBDBYAtgjeNGnSJBivCDvfSI/et2+feCLEqFGjzE+E0GU/DQlKsCkAxgrvYcqUKWB8xYoVE6q2zciYNcIYPHiwl5fXrl27xNlIKhFXwvr777/FK2ENGTKEnK137966zLOJIPMgXOH9lCtXTvjTMKV69eow0oJ3OHPmTIiKIkWKwJaK8Ivfffed8IsQdTDnl19+KZ6GBB8HSs2aNYP3CX8aXgFeR9y4sR3AoDNYM/ytf//9FxaH+OfMee+993SZ14MUTRmAT6fL5Pz58+JEKoDN1ctRcArzAA4ODobI0WW/HKNtHQAYtMEM5OnLQh5DnIvzGDMHc7DsFi1aBIIIG22+mach1apVC15t/PjxJUuWhJVH1MqhAIZg+/rrr2Egu2zZMmHXMXlmM4m5SiTUYrX92Z1YBMJHeOedd+bOnQubjNA75GlIOXYfwgQYwG4BbKH37dv39ddfh1ELhDEYWZ8+fcgvPs0DGLh+/XrLli0LFy4MMdyrV6/AwEDKo7du3Vq7dm2IBHBM80sB6LLvSIThTv/+/SH48+TJA/baokWLf/75RyjlaEYJCQkwvBCOnbHmMpWIa0FDcIIfwXs23/cOwN9t2LAh6PDyyy/Xq1cPvFIY1ly4cKFr164VK1YUrtHRrl27gIAA8bf0ej34KUgHQ3kIAEhiMvNAinfffbdQJjAPfFJxh6TtAAZ27twJAkLAUPJSgAsLFk9OBNlhInwQeHviRCqAzdXLUXAK8wA2Zh05T300GzoAe/bs0WVuxIhThK9IqK8zT548CQsIdNYRq9DBgwdhAA1LDRZu+/btzS/EYWcAT5w4sXHjxqCG0AiwfWl+0LuAuUok5ovVxmd3bhFs3LixQYMGwgVDzC/EYbv7ECbAAEYQBEEQFcAARhAEQRAVwABGEARBEBXAAEYQBEEQFcAARhAEQRAVwABGEARBEBXAAEYQBEEQFcAARhAEQRAVcK8AzsjICAsLi42NfYYgCIKwDDg5+Dl1Dze2cK8AFu6PiyAIgvABuDpt9OzgXgEs3B8XFhi9KcUgUVFRXl5e8C9dcEtQDRJUgwTVIOFJDWFAJd6RgkXcK4BhmcECg3/pAoOkpaXt2LHD2mVs3Q1UgwTVIEE1SHhSgwM/xwBmFZ4ayXVQDRJUgwTVIOFJDQ78HAOYVXhqJNdBNUhQDRJUg4QnNTjwcwxgVuGpkVwH1SBBNUhQDRKe1ODAzzGAWYWnRnIdczUMBgM8TXZL4uLivL294V+64B7Acoelb2PdcGd4UoMDP8cAZhWeGsl1KDVSU1NDQkKuuyvXrl0LCAiAf+mC2wBLH9YBi+uGm8OTGhz4OQYwq/DUSK5DqpGRkXHz5s3g4ODY2NikpCR6fOQGJCYmhoeHw790wQ2AJQ7LHZY+rAPCJRqwU0h4UoMDP8cAZhWeGsl1SDXAhWEMBPFDz+Q2QPA8ffqU6SsEuQgsfVgHYE0wYqdkhyc1OPBzDGBW4amRXMc8gAXzdU8wgMl1ADuFhCc1OPBzDGBW4amRXAcDmAQDGAPYGjypwYGfMxDAkydPbtSoUeHChUuWLNmhQ4ebN2+KpQ8++EBH0Lt3b+L3LMDBAhPhqZFcBwOYBAMYA9gaPKnBgZ8zEMCtW7deuXLl1atXL1261LZt24oVKyYkJAglCOBevXo9ziLHJcHBAhPhqZFcBwOYBAMYA9gaPKnBgZ8zEMAkERERoPjRo0eFpxDAAwcOzD6LLThYYCI8NZLrYACT2A7g+/fvQ+PUrl27Xr16mzZtostcgAFsDZ7U4MDPGQvg4OBgUPzKlSvCU/CREiVKvPrqq3Xr1vX09LR44GtKSsqzLIS7Z0RFRaU5y/4rD3uuPLP7YhhdUBz4sNBI8C9dcEtINeLi4q5duwaPM9wVvV4PAQz/0oVMHjx4cP78eXjw8OHDcuXKgVz0HOwDSx/WAfho1LqBaESNXRfDwEjBTumCI4CTYwArB/TVZ5991rx5c3HK0qVLfXx8Ll++vHbt2tdee+2LL74gZn/O2LFjs74jfo6Xl9cOZ+k5f1elP73bT9tFFxDN4O3tHRAQEB4e/lRLTJs2rXz58rly5erWrRtsNQYGBtJzZAdW5gkTJtBTpQa2XGFzlp7KPrD0YR2ANYFeORBt0G6qyUh/nu+SkYKTYwArR58+fSpVqmTt9suHDh2ChXH79m1qurQj4EuhUbDeVB+5Nzouia4pi0a2ZDWCxkfAFy5cyJ079+7du2H0OWjQoJ49e9JzmAEJXaxYsZiYGLpgB7ZHwCJnz56FAKancgGOgK2hBTXAPMFCwUgDQ6PpmiPgCFg5+vfvDwOIu3fv0oUsEhISYGHAgJguEDxz+TsDg8Hw0Qw/WHW2XbC8HaAYaRx9l+M6pBoa/A544sSJzZo1M2ZeIKJIkSKnTp2i57BEo0aNFi5cSE+1gwyb3wELREdH16lTx9/fny5wAX4HbA0tqLH1fBhY6Ecz/chLdjuB636uOgwEMCwkSN9y5coFBQXRNYITJ07AwoBxA10gkGSBzdp/E9aenqvO0gVl0UIjaQctB3C1atXErz8KFChQsmRJsurl5ZU/f/5Hjx4JT3/44Yd69erFxsbC43Hjxr377rvkzHYiBPDatWutvXJKSsp77723evXqbL/GERjA1tCCGj+tPAsWOuvALbrgIJL4ubowEMB9+/YtWrTokSNHxNONkpKSYPrt27fHjx8fEBBw7969nTt3Vq1a9f3336d/OTuSLLBbT+Jg7Xl9xJ7YJDVXYi00knawEcCwAZeYmi7Hj53b7+Hh4bByzpgxA1bdn376qU2bNmQVXqR+/fq//vorPB4zZkz58uUfPHgglPbt25c3b14IS3L+SZMmFbJCaGioMI8QwHq93uIrw1/85ptvxo4dS7wqb2AAW0N1NWIT08A8wULBSOmag0ji5+rCQACLoweRlStXGjPPpoDELV68eL58+V5//fWhQ4fmuCSkWmAtZx2BFWjTuft0QUFUbyRNYSOAISlhYcnxA6+c7U1YITEx8aWXXhJ2O3fo0AEymJph9+7dsA5PnDixWLFiV69eFacHBgbC6hoSEkLMa9p1HGyF9PTn70fcBW3xlY8fP+7h4fFmFpcvX37x6ryAAWwN1dXYeO4+9M4ns4/QBceRys9VhIEAlhCpFthc3yBYh75ffoYuKIjqjaQptBzAEL0QwMLVY1q1atWvXz96DqOxQYMGMNg9ciSbKwUFBcHqCp+FnGgP5HfAFl+ZezCAraG6Gt8tPwO9M+9gEF1wHKn8XEUwgJ3hdkQ8rEPVhu+JSXh+z1HlUb2RNIWNAFZ9F/Tff/9do0YN4XG3bt26du2avW7a1VygQIFcuXLduHGDnH769GlYXSMjI8mJ9u+Chn+tvTL3YABbQ101wDCrDjftf74TEU/XHEcqP1cRDGAnaTP3GKxG6888tzzlUbeRtIaNAFad3r17d+7cWXg8Y8aMN998k6yeP3/+5Zdf9vLygsHxV199RZb+/fff8uXLk1OMjuyCPnfunLVX5h4MYGuoq4bXmVCwzU/nHqMLTiGhn6sFBrCTLDwcDGtS92Wn6YJSqNtIWkPLAfzOO+/AsFV4fPny5dy5c8fExAhP7927V6ZMmSlTphgzx7seHh6Qx+Iv9ujRw/wLY3vIyDyN2MYrcw8GsDXUVaPbslNgm2CedMEpJPRztcAAdpKQqARYk6p4ekfGZztIVTHUbSStodkAhiwsWLCgt7e3OKVx48ZLliwxZo5la9asSd7Cq23btq1btxYew/svWrSonWcMU0RGRlavXv2XX34Rp5Cv7A5gAFtDRTUi4lLAMME2Q6MsXDPYCST0c7XAAHaedvOPw8q05lS2g1QVQ8VG0iCaDWBzIIxr165t+yoZwOLFiz/55BN6qn3YcyEOvsEAtoaKaqw+FQKG2X7BcbrgLNL6uSpgADvPkiO3YX3qsvQkXVAEFRtJgzAUwMCcOXPu38/hHLZly5aRt752CAxgDGBrqKjG10tOgmEuPUpfLdhppPVzVcAAdp770YmwPlX29A5/poLXq9hIGoStAJYbDGAMYGuopcaTZ8mVM/c/h8VIs//ZKLWfqwIGsEt0XHQCVqkVJ6xeoVo+1GokbYIBTIIBjAFsDbXUWH78LljlF4tO0AUXkNzPlQcD2CUgemGt6ijpWmUnajWSNsEAJsEAxgC2hlpqdFhoGquslHSsIrmfKw8GsEuEP0sWjuu7Hy3ZfhU7UauRtAkGMAkGMAawNVRRQ/i2DqwyPE7KrpTcz5UHA9hVuv5jOrNtkZ80Z7bZjyqNpFkwgEkwgDGAraGKGsJVE7otc+acOhvI4ecKgwHsKsK1XdpIdG0X+1GlkTQLBjAJBjAGsDVUUaP1nKOVZLhuoBx+rjAYwK7yNDFVuLtWcLird9dyCFUaSbNgAJNgAGMAW0N5NYLE+7cmSvxH5fBzhcEAlgCp7i/tEMo3kpbBACbBAMYAtobyaszafxPsseeqs3TBZWTycyXBAJaA7RcewBr24Qw/O2+PIwnKN5KWwQAmwQDGALaGwmqAJX4w/TDY446LD+iay8jk50qCASwBCSnpNUfthZXsyoNYuiYbCjeSxsEAJsEAxgC2hsJqXA6LBWMEewSTpGsuI5OfKwkGsDT0W3se1rNJexy+d7rTKNxIGscNA7hHjx4dOnQQHn/wwQcDBw4USxjAGMDWUFiNid7XwBj7rZPlTlzy+bliYABLw74rj2E9azL5YEaGQnuhFW4kjaPlAH78+PGAAQOqVauWL1++UqVKNWvWbPHixYmJrp44TgZwdHR0XNyLYwAdDWDypcyBqi6T3Llzw/tv2bLl8uXL7X9x+Rg7dix1c2URDGBrKKkGmCFYIhijz9XHdE0K5PNzxcAAlobkNP0bY3xgVTt7L5quyYOSjaR9NBvAd+7cKVOmTK1atTZu3AjvCp7C+2zbtu3OnTvpWTM/BT3JOjZSU/IAbtOmDWxGPHjw4Pz585MmTSpcuPCnn36ani79TkWHwAB2AiXVOHM3GizxjbE+Kel6uiYF8vm5YmAAS8bvGy/B2jZq+xW6IA9KNpL20WwAt27dunz58gkJCdR08Xg9WCFhQNy+ffuCBQtCouj1+p9++qly5cr58+evUaPG3LlzxV+B0uDBg4sWLVq8ePGhQ4d+//33FndBp6Sk/P7772XLloUXbNy4sZ+fnzB95cqV8Ls+Pj6wNVCoUCF4Y48ePTJmxpgwwBUQ5xcxj+dDhw7BnMuWLROeQtj37NmzRIkSL7/88kcffXTp0iVhOjz48MMPIa1h+ltvvXXu3Dlh+okTJ+ANFyhQ4JVXXmnVqlVMTIwxc6Nh8uTJwgevX7/+5s2bhZnh/cDfOnjwYMOGDeFXmjZtKtwkCj4O+bbhqTC/AAawNZRUY+T2y2CJQzY9Xx8kR1Y/VwYMYMk4cisC1rYG4w+k6e0debiCko2kfXIOYIhA8x9yHvMq/CQlWZ3BDqKiojw8PKZMmUIXCGCFLFWq1IoVK2BwHBoaCh9hzJgxkFV3795du3YthCgMnYU5p02bVqxYsa1bt8Kng8CDVLMYwD///HOzZs327NkTFBQ0Y8aMfPnywQNjZmLlyZOnZcuW8OIwkK1du3a3bt1genx8/Ndffy2McYHU1FThdUTMAxiAoScMgoXH8JqwAQEvC39oyJAhr776anS0aT9Q3bp1v/322xs3bsD0TZs2CcF88eJFeEt9+/aFp1evXl2wYEFkZCRMnzhxImwZwPYB6ABvFeY5cuSIMSuA33nnHXh67dq19957Dz4dTE9KSoK/BX9CeNtJ5JLCALaOYmqADYIZgiUevRVB1yRCVj9XBgxgyUjXZzScYFrhDt8Mp2syoFgjMUHOAazTWfhp2/bFDAUL0lX4+eCDFzOUKJGtZAenT5+G9W3btm3iFAinQpkMGzZMmAIzDBo0SJyBon///p06dRIew6B2+vTpwuP09HQYWJsHMER4rly5wsLCxF3QLVq0GD58uDFryHj79vO7sS5atKh06dLCY4sRK2Kx2qVLF4hweHD8+PEiRYrAsFssVatWbenSpfAANhFWrVolThfo2rVr8+bNqYnw67CpcfLki1trwxYGzGkkRsDCdNiwgKfCwsVd0E6gmBqHb4SDGYIlpss2IJHVz5XBLh/hBrkX2NidV2GdG7j+Al2QAcUaiQlYCWAY3gUHBzdu3FgcsMIMMNIVZwAWLlz41ltvlShRAnIaxqxvv/02TIyNjYU5jx49Ks7WsWNH8wD29vaG2YSMF8idOzcMcI2ZAQwhJ/46vCsYnQuPLUasiMUqvGadOnWMme/2pZdeIv8iPBU2LyAg4a/DFsCUKVPE4IfYhiE++VIADIWptw0fHFQyZgVwRMTzUdSFCxfgKWxnGDGAnUIxNQasvwBmCJZIF6RDbj9XALt8hBvkXmDnQ2Ngnas9el9iquzHpyjWSEyQcwCb715OUG0XNLnHGFbI7du3i6X169fnz58fhqeQNBDVv/zyi5Axdgbwhg0bYAQMH//8+fO3bt0KzuTxY9MxqMJ3wOKvwx/VZW1GWIxYEYvVevXqffbZZ/Bg6tSpr732mvCHRIS9ygC8h9mzZ3/yySd58+YVNkRg28I8gIUtlSNHjpAvcv/+fWNWAMOAXpjz4sWL8PTevXtGDGCnUEYNMMBao/aBGV4INX3BLxNy+7kCYABLicFgeG+aXJd9oVCmkVgh5wBWiVatWkE+UQdh2QjgX3/99eOPPxafwvBRzBhqF3SFChXMAxgCT0gy86OgbQRwr1692rVrJ5YozANYOAhrxYoV8PjAgQMQ+UIi2uCbb75p3749PPjhhx/Md0HHxcXly5dv9erV1HSjzQCeNGnSG2+8Qc4sggFsDWXUAAMEG3x/+mFZLw4ot58rAAawxMzMvPDpTyulv/AphTKNxAqaDeDbt2+XLl26Vq1aMDaFd3Xz5s01a9bAlN9//12YgQrgefPmFSlSxMfHB6J01KhR8FgMYBhrFi9eHGa+ceMGRKa1g7C6d+9euXJlCDP402fOnJk8ebK3t7fRZgBDklWsWBHeG4xczdeoHpZOQ4LA1utN55aAw7777rvwJvfv3w+56O/vP2LEiHPnziUlJfXv3x/iMyQk5MSJE9WqVRP2S8PngtFw3759AwMD4YMsXrxYGC6PHDny1VdfXbVqFbxt+Cvz588Xvj+2EcDr1q0rVKgQTIFXIL+ENmIAW0cZNX4ULo+/33S8unwo4OdygwEsMcHhplt/VBu+JzqBPppUWpRpJFbQbAADjx49gnFtlSpV8uTJA9HVuHHjGTNmiBfioAIYggTGiJCUr7zyCqSUp6enGMAw6oWUhUiGEuS3tdOQQITRo0dDoMKfg0HzF198cfnyZaPNAI6IiPjkk0/gvemsnIakyyR37twlS5Zs2bIljH3J4TWMX3/77bdy5crBX4RxOWwB3L9/PzU1FUa98BTiFkqggLhEYIDerFkzGPLCB2ndurUQrhDkc+fOrVmzJrwI/BWYLuxvtxHAoFWnTp3gRXR4GpLdKKAGWB8YYCXTDeLi6ZqkKODncoMBLD1t5x2DlW/NqRC6ICkKNBJDaDmAlcfRC3HwBwawNRRQY/WpEDDAz+bLfot0ZfxcVjCApWfp0duw/n31tz9dkBQFGokhMIBJMIAxgK2hgBqdFvuDAf5z9A5dkBpl/FxWMICl51FsUmVPb1gFw2Jcvd6vDRRoJIbAACbBAMYAtobcatyPTgTrAwN8HCt7Ayrj57KCASwLXZaehLVwkV8wXZAOuRuJLTCASTCAMYCtIbcaCw8Hg/V9s/QUXZABxfxcPjCAZcHrTCisha3nvDhrU3LkbiS2wAAmwQDGALaG3Gq0mn0UrG/9GdOVUuRGMT+XDwxgWYhNTHt9hOk4wOuP5PpbcjcSW5gHMHVxYLcCAxiWPgawRWRVA+wOTK/6iL1ggHRNBhTzc/nAAJaLXv+dg3Vx8t7rdEEiZG0k5iDV0Ov1YL5RUVH0TG4DBjAsfVgHhDOVsVNIZFVj8p7rYHpgfXRBHpT0c5nAAJaLvZcfwbr4zqSD+gxZrgUjayMxB6XGo0ePhAyGkVCy+5GYmBgeHg7/0gU3AJa4kL7CzRbN1w03Rz41wOgaT/IF09t35bnycqOkn8sEBrBcJKfp6431gdXRP/j5dXGlRb5GYhFKDYPBIGSwe3Lt2rWAgAD4ly64DbD0xYsgYqeQyKfGieBIsDswvZR0044HBVDSz2UCA1hGPLcGwhr5hzz3o5avkVjEohp6vZ4eH7kHcXFx3t7e8C9dcA+EPc8iFtcNt0U+NYZsugR257nVdNk1ZVDYz+UAA1hGztyNhjWy7hifpFTpNwnlayQWQTVIUA0SVINEJjXA4uqMNt3+6Oy9aLomGwr7uRxgAMtIRoah2ZRDsFLuuvSQrrmMTI3EKKgGCapBgmqQyKTGzksPweiaTz2UIc8hLxZR2M/lAANYXmb4yHVzJJkaiVFQDRJUgwTVIJFJDeH2RzNlvv0RhfJ+LjkYwPIi3hwpKj7b7dJcR6ZGYhRUgwTVIEE1SORQIzI+paoitz+iUN7PJQcDWHbazT8Oq+Yq/3t0wTXkaCR2QTVIUA0SVINEDjVWnrgLFtd+wXG6IDOq+Lm0YADLzvLjprXz84Un6IJryNFI7IJqkKAaJKgGiRxqfL7ANMZYceIuXZAZVfxcWjCAZSci7vn+mTsRUu6fkaOR2AXVIEE1SFANEsnVuB0RD+YGFhcp9bdsOaKKn0sLBrAS9FhxBtbRWQdu0QUXkLyRmAbVIEE1SFANEsnVmLXfdJzpDyvO0AX5UcvPJQQDWAl2XHwA6+i706Q8Rl/yRmIaVIME1SBBNUikVQMMrflU05mWYHF0TX7U8nMJwQBWgqRUfd0xpstSnrkr2Vnq0jYS66AaJKgGCapBIq0ap+9EVZLtWkM5opafSwgGsEIM3Wy6TtufWwLpgrNI20isg2qQoBokqAaJtGoM22y62i78SxcUQUU/lwoMYIU4lbmp+MYYn+Q0aTYVpW0k1kE1SFANElSDREI1xB17MA6ma4qgop9LBQawQohfluyU6LKUEjYSB6AaJKgGCapBIqEachza4hAq+rlUYAArh3C44PfLpTlcUMJG4gBUgwTVIEE1SCRU47vl0p/c4RDq+rkkYAArx73IBFhfq3h6hz9LpmuOI2EjcQCqQYJqkKAaJFKp8eRZMlgZGFpIVAJdUwp1/VwSMIAV5cvF/rDKLj16my44jlSNxAeoBgmqQYJqkEighl5v9PPzHT2nS9fJnRceo6sKorqfuw4GsKKsOx0KAdxq9lGDwdVvTSRoJI5ANUhQDRJUg8RVNbZuNZYvb9TphJ+EUmVNU1RCdT93HQxgRYlNSqs+ci9k8JUHsXTNQVxtJL5ANUhQDRJUg8QlNSBrPTzE9IUfAzyFH5UyWHU/dx0MYKXpt+48BPC4XdfogoO41EjcgWqQoBokqAaJ82ro9eTY98UPBHCFCqaq4mjBz12EgQCePHlyo0aNChcuXLJkyQ4dOty8+eKez8nJyf369StevHihQoW+/PLLJ0+eEL9nAS0ssEM3nkAAN5xwIE2fQdccwflG4hFUgwTVIEE1SJxXw8+Pjl7yB6qKowU/dxEGArh169YrV668evXqpUuX2rZtW7FixYSE58fd9enTp0KFCocOHQoICGjSpEmzZs2y/yqNFhZYuj4D0hcy2PdaDpsLtnG+kXgE1SBBNUhQDRLn1fDyokOX/IGq4mjBz12EgQAmiYiIAMWPHj0Kj2NjY/PkybN582ahdOPGDSidOnUq2y9kRyMLbMLuaxDAvVcH0AVHcL6ReATVIEE1SFANEufVwBGwDDAWwMHBwaD4lStX4DEMfOHx06dPxSoMjmfPnv1ibjM0ssBuPH4GAfz6iD3RCal0zW6cbyQeQTVIUA0SVIPEeTUyvwM2mEcvfgfsAiwFcEZGxmeffda8eXPh6bp16/LmzUvO8Pbbbw8bNoycAqSkpDzLIiwsDBZYVFRUmtp8Nu8YZPCyo8F0wW4SExOhkeBfuuCWoBokqAYJqkHiihrpGzdSAWzw8IAfmE7Pqgjg5BjAytGnT59KlSpBiApP7QzgsWPH6rLj5eW1Q21+/3snBHCzcbvpAoIgiCbZv2yZPvtpSIklSpz58096PqUAJ8cAVoj+/fuXL1/+7t274hQ7d0FrcwQcEZv4+og9kMGXQp18M65syfIHqkGCapCgGiSuqBHe5zcIXf+K9Z/u3pe+enW6r29acjI9k4LgCFgJDAYDpG+5cuWCgoLI6cJBWFu2bBGe3rx5U8fIQVgC/daaTgj+a9dVumAfaU5/l8MjqAYJqkGCapA4r0ZKSlLhIhDAi4fMoUsqoSk/dw4GArhv375FixY9cuTI4yySkpKEUp8+fWDUe/jw4YCAgKaZZP9VGk0tsMM3wiGAG4w/kJruzAnBzjcSj6AaJKgGCapB4rQaYFNf9l+6qMlXh68/pmsqoSk/dw4GAlj87lZk5cqVQkm4EEexYsUKFiz4xRdfQDZn+00zNLXA0vUZb0/0hQzed+URXbMDpxuJS1ANElSDBNUgcVqNvZcfgVk1nuSrV+nuv+Zoys+dg4EAlhCtLbDJe6/DOt1z1Vm6YAdONxKXoBokqAYJqkHipBoZGT+tPAtmNWXvDbqkHlrzcyfAAFaT4PA4WKerDt8THufwHYKdbCROQTVIUA0SVIPEOTVSWrXZVvej5r3/DQ6Pp2vqoTU/dwIMYJXpsPBEJafuEOxcI/EKqkGCapCgGiTOqHHlilGn03u81Gv888sOagQN+rmjYACrzNrTIRDALWcdcfQOwc40Er+gGiSoBgmqQeKEGoaffoIA3lOj2brToXRNVTTo546CAawyz5LTamTeIfhCaAxds4kTjcQxqAYJqkGCapA4rEZEREa+fBDA3/SYAWZFV1VFg37uKBjA6jNow0UIYM+tgXTBJg43EtegGiSoBgmqQeKwGhMmQPpeKlN98PoLdElttOnnDoEBrD4nb0dBANcd45OYmk7XrONwI3ENqkGCapCgGiSOqZGamlGmDATwgHZDTt2Joqtqo00/dwgMYPUxGAzvTz8MGbzp3H26Zh3HGol3UA0SVIME1SBxTI01ayB9Hxcu3mLyfkcPUlEAbfq5Q2AAa4IFh4IggL/6258uWMexRuIdVIME1SBBNUgcUyM6ek3HviNa9Vt4OJguaQDN+rn9YABrgsexyVU8vSGD70TYe5qdY43EO6gGCapBgmqQOKTG7Yh4MCWwJjAouqYBNOvn9oMBrBV+WHHGoQvNONRI3INqkKAaJKgGiUNqCJfq+3GlM5fqUwAt+7mdYABrhX1XTJdabTTRN11v170ZHGok7kE1SFANElSDxF417t7NaNZ8WJdRlUwXq8/hGvtqoWU/txMMYK2Qmp7x1vgDsLofuPaErlnC3kZyD1ANElSDBNUgsVeNwYONOt2RKm81nHAgzb4hgfJo2c/tBANYQ0z0vlbJdG+Gc3TBEvY2knuAapCgGiSoBoldaoBDvvwyBPD3ncdN2nOdrmoGjfu5PWAAa4gX92Z4lvMhD3Y1ktuAapCgGiSoBoldasydC+kb/Gr5ysN2aeruCxQa93N7wADWFl8sMt2bYbFfzvdmsKuR3AZUgwTVIEE1SHJWQ683Vq0KATyiVb8vFztwYqTyaN/PcwQDWFtsPHsfAviD6YdzPO0950ZyJ1ANElSDBNUgyVmNHTsgfZ8VeLnW4C1gR3RVS2jfz3MEA1hbJKam1x3jAxnsfzuSrmUn50ZyJ1ANElSDBNUgyVmNDz6AAF7U5Ks3HLw4rvJo389zBANYcwzfdhkC+FevHC59nnMjuROoBgmqQYJqkOSghsFg3Lw5qNZb7/RdNWLbZbqqMZjwc9tgAGuOy2GxEMDVR+yNTkilawQ5NJKbgWqQoBokqAZJjmqA7YD5gAVdeRBL1zQGE35uGwxgLfLZ/GPQAMuO3aELBDk2kluBapCgGiSoBkmOaoDtgPm0m3+cLmgPVvzcBhjAWmTNqRDogRazjtg4FCvHRnIrUA0SVIME1SCxpcaCBYbJkzuO2wHms/Z0CF3VHqz4uQ0wgLVIXHJarVH7oA3O3Yuma1nYaiT3A9UgQTVIUA0Sq2okJxtLlTLqdP0+/xPMByyInkF7sOLnNsAA1ih/bLoEATx440W6kIXVRnJLUA0SVIME1SCxqsaKFZC+Ma+WqTp059DNl+iqJmHIz62BAaxRzofGQADXHLU3NsmsVTKx2khuCapBgmqQoBokltUwGIz160MAT//4R7CdC6Ex2apahSE/twYGsEYxGAytZh+FZvjv5D26lonlRnJXUA0SVIME1SCxrMbhw5C+afkL1Bu4ofWcozYOPdEUDPm5NTCAtcuKE3chgNvMPWaxHyw3kruCapCgGiSoBollNdq3hwDe1bwjGM7KE3ezlTQMW35uEQxg7RKbmFZjpOmEPIt7hCw3kruCapCgGiSoBokFNYKCjB4eEMAf/bwEDAds50VJ27Dl5xbBANY0gzdehAD+Y5OFYyIsNJIbg2qQoBokqAaJBTXu3DF+//3ltz8Cq/l9owWr0SzM+bk5GMCaJiAk61Ass81SC43kxqAaJKgGCapBYlENsJeaI7zBas5b2tmmWZjzc3MwgDWNwWBoPcd0KNYKsy9mLDaS24JqkKAaJKgGiUU1lh83HW7C0OFXAsz5uTkYwFpn9cl70Bstza6KZbGR3BZUgwTVIEE1SLKpkZ5uHDLEcPlyi1lHwGRWn2Lg6lckLPo5BQaw1olLTqs92nRVrDN3s10VC22FBNUgQTVIUA2SbGps3mw6+6j4q6//sR1MhomrX5Gw6OcUGMAM8OeWQAjgAeuz3aAQbYUE1SBBNUhQDZJsajRvDgG8r+PPYC+eWwPpWTUPo35OggHMABZvUIi2QoJqkKAaJKgGyQs1zp6F9DXkydP8tzVgL2Ay9Kyah1E/J8EAZoP2C45Dkyw9elucgrZCgmqQoBokqAbJCzW6d4cAvvVJBzCWzxcwcPNBc9j1cxEMYDZYfyYU+uSD6YczMp4fioW2QoJqkKAaJKgGyXM1QkKMuXNDAP88YAkYy4azofR8LMCun4tgALNBQkr6G2N8oFVOBEcKU9BWSFANElSDBNUgEdTQe3pC+sY2agKWAsaSmJpOz8cC7Pq5CAYwM4zecQW6pffqAOEp2goJqkGCapCgGiTPA3jmTGOpUn//Nh0sBYyFnokRmPZzAQxgZrj5OA66perwPY9jk41oK9lBNUhQDRJUg0RU49GTmNf/3AWWcutJHD0TIzDt5wIYwCzReclJaJhZB24Z0Vayg2qQoBokqAaJqMas/TfBTL5ecpKegx1Y93MjBjBb7A58CD3TaKJvmj4DbYUE1SBBNUhQDZJ0f//TI0YkxCeCjYCZeAc+oudgB9b93IgBzBap6Rli26CtkKAaJKgGCapBktGmjVGnu9FroLgpT8/BDqz7uREDmDnEHUdoKySoBgmqQYJqvODGDdPFNzw8+o/fKH6ZxS4c+DkGMGM8ik2qOnwPNM+1sBi0FRE0WRJUgwTVeEHfvhDAdxs2Jg/nZBcO/BwDmD16rw6A/hm5LRBtRQRNlgTVIEE1nhMdbSxYEAJ4Yr9JYCB91jw/oZFdOPBzDGD2OBEcCf1TZ8y+DVvRVp6DJkuCapCgGs+ZOhXSN71evRqeprOP/LMu6cMuHPg5BjB7GAyGj2b6QQv9sWQn2ooAmiwJqkGCapiAj1++PATwsRGmi298NMOPur84i3Dg5xjATLLixF3ooiZ/7U5NfXF/JHcGTZYE1SBBNUyEhhrfestQqtSn0/aDdfx77MVtXdiFAz/HAGaS2KS0WqP2QiMdv/WErrklaLIkqAYJqvEcgyHghOlyttWH746OS6SrDMKBn2MAs8qfWy5BL/3y31m64JagyZKgGiSohohw/Ga32bv4UIMDP8cAZpWrYdHCuQQPnybRNfcDTZYE1SBBNYze3uB9D54mVfH0BtP424sTNTjwcwxgVoEWajnJdDTjtH036Jr7gSZLgmqQuLsa9+8bc+UyvvzyfK/jYBffLD3JjRoc+DkGMKtAC41bvhM6qsH4A8lperrsZri7yWYH1SBxdzWGDTPqdPoPP/zfONPhV3suPeBGDQ78HAOYVaCFtm7f0XTyQWiqTefu02U3w91NNjuoBolbq5GQYHzlFdPZR7OWg1E0m3IoOSWVGzU48HMMYFYRbGXBwVvQV5/NP8bBWX2u4NYmawaqQeLWaixebLr4c7Vqn805Yvr298htntTgwM8xgFlFaKTw2IQaI03nIwWERNNzuBM82YrroBok7qtGRoaxRg0I4NCxU8AiwChiElJ5UoMDP2cggI8ePdquXbuyZcuC1tu3bxen9+jRQ0fQunVr4pcsw8ECExEbaehm0/lIv3pdoOdwJ3iyFddBNUjcV409eyB9jUWK/P7vMbCIYZsDjXypwYGfMxDAe/fuHTly5LZt28wDuE2bNo+ziImJIX7JMhwsMBGxka4+jIXuqjZ8z5NnbN/bxBV4shXXQTVI3FeNyZONL72U8OuAasL90x6afI8nNTjwcwYCWMQ8gDt06EDUc4aDBSZCNtJXf/tXYv/unq7Ak624DqpB4tZqhIT8veEEmEPnv08KE3hSgwM/ZzuAixYtWrJkyRo1avTp0ycqKoqY1zIcLDARspF2Bz6EHms4wX3PR+LJVlwH1SBxZzXAEN4afwDMwTvwkTCFJzU48HOGA3j9+vU7d+68fPkyTKxdu/bbb7+t11uIn5SUlGdZhIWFwYtAVKexT2JiIjQS/AuPk5JTmmSej7T+9D16PveAVANBNUjcUY0nT9KuXoX/vU7fA1sAcwCLECo8qQFOjgGsHFQAk9y5cweqBw8epAtG49ixY3XZ8fLy2sEdvy0yXZSj2bjd27fTJQRB3Irr3bsbPDyCOnYEQwBbAHOg5+ACcHIMYOWwEcBAiRIllixZQk91jxEwEBGbKN4fKfuMbgFP2/Wug2qQuJ0aCQmGsmXBLm/OWASGALYQ+ezFZ+dJDRwBK4qNAIZk9fDw2LlzJ13IzjP2vzMQSTP7Lmfk9svQb73+O0fM5S6Yq+HOoBokbqfGmjWms4/Klu29/CQYwqjtV8giT2pw4OcMBHB8fPzFTEDr2bNnw4PQ0FCY+Mcff5w6derevXsHDx586623qlevDoNd+pezw8ECEzFvpODweOi3yp7eoVE83OzTIczVcGdQDRL3UsNgML71FgRwzIgxlTPvfXQ7Ip6s86QGB37OQAD7+fllfXv7nB49eiQlJbVq1apkyZJ58uSpVKlSr169njzJ+db0HCwwEYuN9P3yM9By43ZdIye6AxbVcFtQDRL3UuP4cdPwN3/+aWtN9z7qseIMVedJDQ78nIEAlhAOFpiIxUY6cisCuq7uGJ+4ZB4azH4squG2oBok7qVGp04QwKk/9QQTACs4eiuCqvOkBgd+jgHMKhYbyWAwtJhluur68uN3yencY1ENtwXVIHEjNeLjjSVLQgBvWW1KX7AC83u08KQGB36OAcwq1hpp7ekQ6L33px/WZ9C9xzHW1HBPUA0S91IjIUG/bft70w6DCaw7HUpX+VKDAz/HAGYVa42UlKp/M/PO2z5XH1MljrGmhnuCapC4mxr7rjyG9gcTACuga3ypwYGfYwCzio1GmrbvBnTgV3/70wV+saGGG4JqkLiLGmFhpkOgjcZOi01Xhp/uc4OeIROe1ODAzzGAWcVGIz15lvz6CNP9Ty7ef0rXOMWGGm4IqkHiFmro9cZq1Yx16lz3NZ37W33E3nAr90bjSQ0O/BwDmFVsN9LvG003Ce637jxd4BTbargbqAaJW6ixc6fp7KNixQb9a7r30ZBNl+gZsuBJDQ78HAOYVWw30vVHz6APq3h63492i4ty2FbD3UA1SNxCjQ8/hAB+NuD3KpkX37jx2KrF8aQGB36OAcwqOTZS92WnoRXH73aLi3LkqIZbgWqQ8K/GxYum4W+uXHOWm26J9u2/p+kZCHhSgwM/xwBmlRwbye9meKXMi3I8c4OLcuSohluBapDwr8YPP5guvtH56zqj90HLHzG7+AYJT2pw4OcYwKySYyMZDIZPZpsuyrH06G26xh05quFWoBoknKvx5Ikxb14I4G2Lt0Czt5p91PziGyQ8qcGBn2MAs4o9jbTx7P1KmbfjTtNn0DW+sEcN9wHVIOFcjXnzIH0z3mkCbQ7NvvHcfXqG7PCkBgd+jgHMKvY0Ukq6vuEEX2jLHRcf0DW+sEcN9wHVIOFcDRjv+voeW2Ya/jaa6AstT8+QHZ7U4MDPMYBZxc5Gmn8wCDqz7bxjtndMsY6dargJqAYJ92pAa0ODQ5svOBRE18zgSQ0O/BwDmFXsbKSYhNRao0yHZhwPiqRrHGGnGm4CqkHCrRqwSZ1oOsnwWJDpHmjQ5tDs9Dxm8KQGB36OAcwq9jfS2J1Xczw5gXXsV8MdQDVIuFXj8GFjsWLGceOEEw7/2nWVnsESPKnBgZ9jALOK/Y10Pzqx6nDTlSmvPIila7xgvxruAKpBwq0an39u1Omivu8JrQ0NHhZj1yV3eFKDAz/HAGYVhxppwPoL0KW/eV2gC7zgkBrcg2qQ8KlGcLDRwwMCeMKMrdDaA9fb29o8qcGBn2MAs4pDjXT1YaywmczrlSkdUoN7UA0SPtUYMADSN6llK+Hak9ce2utpPKnBgZ8rF8DXr18fM2bMRx99VLVq1TJlytSrV+/7779ft25dSkoKPatscLDARBxtpG//NX1RNHanXV8UMYejavANqkHCoRqxscbChSGAV/21FJr6u+Vn6Bmsw5MaHPi5EgF8/vz5Fi1a5MuX7+OPP/b09Jw/f/6yZcumTZv23XffVa5cuXjx4lOnTlUmhjlYYCKONtKJ4Ejo1Zqj9kbbcagkcziqBt+gGiQcqjFrFqRveu06NUeaju3wD3bgBAee1ODAz5UIYEjZRYsWPX1q+d60J0+e7NKly6RJk+iCDHCwwEQcbSSDwfDZfNPJgnN9g+ga+ziqBt+gGiS8qWEwGKtWhQA+MHgitHO7+ccdOsWfJzU48HMlAtiehW3PPK7DwQITcaKRdl16CB3bYPyBpNQcLpfDHE6owTGoBgmHagQFpQ3+/Z1Ru6Cddwc+pKs24UkNDvxciQC2TWLmueTKwMECE3GikdL1Ge9OOwRNu8r/Hl1jHCfU4BhUg4RLNVaeuAuN/N60w+kOXuadJzU48HNFA/jjjz9+8CDbRYnPnDlTvXp1coqscLDARJxrpNWnQqBvm005xNntGZxTg1dQDRKu1MgwtW1qekbTzFsvQDvTM+QET2pw4OeKBnDbtm2LFy++YcMGo2lFyhg7dmyePHkGDhxIzycbHCwwEecaKTnt+e0ZNgeE0TWWcU4NXkE1SLhSo0cP41df7dtoSt9GE32hnekZcoInNTjwc0UDGFi4cGHBggW7du3atGnTcuXK7d+/n55DTjhYYCJON9LfR25D93480y8jw4FjNzSO02pwCapBwo8aDx8ac+c26nS/DFoCLbzkiDP3+eZHDS78XOkABjw9PT08PGDs6+/vT9dkhoMFJuJ0I8Ulp9Ub6wMNvO/KI7rGLE6rwSWoBgk/aowcCekb3aAxNC+0cHxKOj2DHfCjBhd+rmgAx8TEfPnll0WLFv3nn3+6d+9eqFChRYsW0TPJCQcLTMSVRprhcxN6uP0Cx05g0DKuqMEfqAYJJ2okJRlffRUCePJPE6B5Z+6/Sc9gH5yokQkHfq5oAJcrV6558+Z3794Vnm7YsKF48eJt27bNPpeMcLDARFxppKj4lJqj9kIbHwuKoGts4ooa/IFqkHCixrJlkL7Jr1WsMnRnrVH7nL6cDidqZMKBnysawOPHj8/IPIpPJCwsrGXLluQUWeFggYm42EjCPQq/WXqKLrCJi2pwBqpBwoMaBoOxbl0I4P+++q2S3XcetAgPamTBgZ8rGsCqw8ECE3GxkR4+TaqWeY/C86ExdI1BXFSDM1ANEh7UOHAA0ldfqPAbgza+PmIPNC89g93woEYWHPi5EgEcGhpKTzKDOj9YJjhYYCKuN9Ifmy5BAPdcdZYuMIjravAEqkHCgxpgWXPnbvmqHzTs0M2X6Koj8KBGFhz4uRIBXKpUqV9++eXsWQtGHxsb+88//9StW3fevHl0TQY4WGAirjfS7Yj4yg7ey0yzuK4GT6AaJHyoIdxRFBoW2pauOQIfaghw4OdKBHB0dPTgwYOLFi1aunTptm3b/vzzz7/++mv37t0bNGiQN2/eJk2a7Nmzh/4deeBggYlI0kj9152Hru67NoAusIYkanADqkHChxp91gRAq/7qdYEuOAgfaghw4OdKBHBgYKBer09KStqyZcvAgQM7duzYunVrCOCZM2deuXKFnltOOFhgIpI00o3Hz4TN6qAncXSNKSRRgxtQDRK21YiONjZv/mj+ksrDTLdeuPnY1T5lW43scODnSgTwSy+9FBFhOt2lSpUqUVFRdFlBOFhgIlI10i+rz0FjD1jv6pa1ukilBh+gGiRsqzFtmlGnC6tUo9Kw3b1XS7Cnim01ssOBnysRwMWLFz99+jQ88PDwEJJYLThYYCJSNdKVB6bvlqp4et+NTKBr7CCVGnyAapAwrAa85/LlIYD/aDsImhRalZ7BcRhWwwwO/FyJAO7Vq1e+fPkqV64MQ+GKFStWMYP+BdngYIGJSNhIP648C+09ZJNLR1eqi4RqcACqQcKwGhs2QPrGFS1eY8i2n1ZaOIjVCRhWwwwO/FyJAAb27du3YMECGAFPmDBhrhn03LLBwQITkbCRLoTGQABXHb7nfrRy92aWFgnV4ABUg4RhNZo0gQCe8253aM+L95/SVadgWA0zOPBzhQJY4IcffoiLc/UgAlfgYIGJSNtI3/57Gprcc+tlusAI0qrBOqgGCatqnDoF6ZueO0/DX9d8t/wMXXUWVtWwBAd+rmgAqw4HC0xE2kY6ey8aAtjFi+yoiLRqsA6qQcKqGl26QABvrv8JNOa5e9F01VlYVcMSHPg5BjCrSN5IXZaehFYfvUPRE8OkQnI1mAbVIGFVjT17btd7p82P86W9YDuraliCAz/HAGYVyRvp5O0oCODqI/ayOAiWXA2mQTVIGFUD2hCaEVry1B0pT91kVA2LcODnGMCsIkcjCYPgkdvZ+yZYDjXYBdUgYVSNEdsuV5LhfmWMqmERDvwcA5hV5Ggk2NYWvgl+wNogWA412AXVIGFPjTVrng0b8c6AtdCMpyUd/hpZVMM6HPg5BjCryNRIsMUNbT98G2ODYJnUYBRUg4QxNTIyjDVrGnW6v1r06vqPxMNfI3Nq2IQDP8cAZhWZGunM3eeHQ4fFsHROsExqMAqqQcKYGnv3mi6+ka9gnUGbzkp38LMIY2rYhAM/xwBmFfkaqdsy0yDYc2sgXdAw8qnBIqgGCWNqtGoFAfzP2x27LzNdvldyGFPDJhz4OQYwq8jXSMI5wdWYujCWfGqwCKpBwpIaV69C+uo9Xmre+18Jz/0lYUmNnODAzzGAWUXWRoKtb8jgP7cwMwiWVQ3mQDVIWFKjVy8I4L01mn37ryzDXyNbauQEB36OAcwqsjZSQMjzQXBIFBu3SJJVDeZANUiYUSMyMiN/fgjgTt2nBYTE0FWJYEYNO+DAzzGAWUXuRvp++RnI4EEbLtIFTSK3GmyBapAwo8ajR/4ff3GmfJ0ey+Ua/hoZUsMOOPBzDGBWkbuRLoeZ7hNc2dM76Ima98+wE7nVYAtUg4QVNW49iYN2qzxslyT3/bUGK2rYAwd+jgHMKgo0Uu/VAZDBfdYE0AXtoYAaDIFqkLCihjLtxooa9sCBn2MAs4oCjSRskoMpyLpJLgkKqMEQqAYJA2oYDOG//fFZj7kK7HBiQA274cDPMYBZRZlGGrj+AgRwjxWS3Y5UJpRRgxVQDRIG1Dh+3KjTJefO67n8GF2SGgbUsBsO/BwDmFWUaaR7kQlVh++pJOkdSeVAGTVYAdUg0b4a0W3aQwCv/19rBU460L4a9sOBnzMQwEePHm3Xrl3ZsmVB6+3bt4vTDQbD6NGjy5Qpkz9//hYtWgQFBRG/ZBkOFpiIYo3kuTUQAvjrJSdBcLqmGRRTgwlQDRLtqqHXG/38DHPnZug8IIBnz9lGzyAD2lXDcTjwcwYCeO/evSNHjty2bRsVwFOnTi1atCisTIGBgZ9//nmVKlWSk5OJ37MABwtMRLFGEu9Lejwokq5pBsXUYAJUg0SjamzdaixfHnJX+EnJlSdmtRc9jwxoVA2n4MDPGQhgETKAYTQGY98ZM2YIT2NjY/Ply7d+/foXc1uCgwUmomQj/bXrKgRw+wXHNTsIVlIN7YNqkGhRDUhfD9OoV/wxwL8wBabLjBbVcBYO/JzVAL5z5w48vXjxxWUi3n///QEDBohPRVJSUp5lERYWBr8VFRWVxj6JiYnQSPAvXZCBxzHxdUbvgwzeeeE+XdMGSqqhfVANEs2pkZxseO01U+Jm/zF4eBjKl4cqPb+kaE4NFwAnxwBWDjKA/f394emjR4/EaufOnb/++mvxqcjYsWN12fHy8tqBOEjvBTshgBuN3b11G11CEMR+jk+YQEUv+QNV+hcQK4CTYwArh3MBjCNgSYiJT2owfj9k8KoTd+iaBlBYDY2DapBoTY301avNc1f8gSr9C5KiNTVcAUfAikIGsP27oEmesf+dgUia4t/lrDxx1zQInuibmJpO19RGeTW0DKpBojk1/PzMc/fFD1TlRHNquAAHfs5qAAsHYc2cOVN4CssAD8KSm9T0jHenHYIMXng4mK6pjfJqaBlUg0Rzauj1xvLlzb8DNh2EVaGCqSonmlPDBTjwcwYCOD4+/mImoPXs2bPhQWhoqDHzNKRXXnll586dly9f7tChA56GpADbLzyAAH5jjE9MQipdUxVV1NAsqAaJBtVI6tvfKBz5TKYvHgXtIBz4OQMB7Ofnp8tOjx49jFkX4ihdujSMfVu0aHHr1i36N83gYIGJqNJIGRmGNnOPQQZP2H2NrqmKKmpoFlSDRHNqREUlFCkGoZuQv9CLAIaxr/zpa9SgGi7AgZ8zEMASwsECE1GrkQ7fDIcArj5ib1hMIl1TD7XU0CaoBonW1Ijv/j0k7s0SFY9cDDF94+vlZfpX5j3PIlpTwxU48HMMYFZRq5EMBkOXpSchgweuv0DX1EMtNbQJqkGiLTWOHBGGvKOHLVXlsjbaUsM1OPBzDGBWUbGRLofFQgDDT2DYU7qmEiqqoUFQDRINqZGSklytOqTv2v+1UesWnxpSw2U48HMMYFZRt5EGbbgIAQxDYVW24s1RVw2tgWqQaEcNw44dkL4RhV4ZLv9tB62hHTVchwM/xwBmFXUbKSwmsfpI0x0afK89oWtqoK4aWgPVINGOGgeuPencbWrPLuMePE2ia0qhHTVchwM/xwBmFdUbacreGxDAH830S9Nn0DXFUV0NTYFqkGhEDWgTaBZoman7btA1BdGIGpLAgZ9jALOK6o30LDntf+NMF6dcfSqErimO6mpoClSDRBNqHDiwbYvp/L0G4w9A49BVBdGEGhLBgZ9jALOKFhpplf898JSGEw7EqeopRm2ooR1QDRL11YiMzHj11eQ8+T7/btZ/J+/RVWVRXw3p4MDPMYBZRQuNlKbP+HCGaa/adB8196oZtaGGdkA1SNRXo0cPo053vWTlllN9Vf++Rn01pIMDP8cAZhWNNJLP1cem63KMVPm6HBpRQyOgGiQqq3HoEKRvhs6j47czoVnoquKorIakcODnGMCsopFGEq/L0W/debqmIBpRQyOgGiRqqpGcbKxuOvH3vwaffbP0lBbO2VNTDanhwM8xgFlFO4107eGzyp6m63KcvRdN15RCO2poAVSDRE01xoyB9H1SuHj9QRuvP9KE7aiphtRw4OcYwKyiqUby3BoIAdxu/vGMDHW28TWlhuqgGiSqqXHzpiFPHgjgPh08PbdepqsqoZoaMsCBn2MAs4qmGikyPuWNMT6QwRvP3adriqApNVQH1SBRTY3k5Ks9B+6r0fSN0fugQeiqSqimhgxw4OcYwKyitUZaevQ2BHCjib7xKel0TX60poa6oBokaqkBjdBwgm+lYbv/OXqHrqmHWmrIAQd+jgHMKlprpJR0/fvTD0MGT1PjQj9aU0NdUA0SFdSIizOmp0/dZ7pU3AfTD6emq3zqEYkKasgGB36OAcwqGmyk/VmnJIVGKX1KkgbVUBFUg0QFNb79NrX+/9r3XAjtcEAbF0sXUUEN2eDAzzGAWUWDjWQwGLotOwWm8/N/5+iazGhQDRVBNUiUVsPX13Tir4fH59/NgnbQwqlHJEqrIScc+DkGMKtos5GCnsRVG74HMvjwzXC6JifaVEMtUA0SRdVISjK+/joE8IqG7aERoB3oGdRGUTVkhgM/xwBmFc020oTd1yCAP5zhl5Kup2uyoVk1VAHVIFFUjZEjTXf8LVqi7qBNE72v0VUNoKgaMsOBn2MAs4pmGykuOc108Oef3ov8gumabGhWDVVANUiUU+PqVWPmib+/fDGi0URf1e9QYhHl1JAfDvwcA5hVtNxIW8+HQQDXGrXvUaxCNx7XshrKg2qQKKRGRoaxeXNI34M1mlQatnvbhTB6Bm2gkBqKwIGfYwCzipYbyWAwdFrsDxn8q9cFuiYPWlZDeVANEoXUiIoyNmuWkq9A074rvvrbX2vHXokopIYicODnGMCsovFGuvIgVrhA9MnbUXRNBjSuhsKgGiSKqeEfFN7qp4VVPL2vPoyla5pBMTUUgAM/xwBmFe030sjtlyGAW846osCFCLSvhpKgGiTKqAErOazqsMKP2n6FrmkJZdRQBg78HAOYVbTfSE8TUxuMPwCWtNjvNl2TGu2roSSoBonsavj4GIcMWbrXdD8SWOFhtadn0BKyq6EgHPg5BjCrMNFImwOeH40VFiPvtbGYUEMxUA0SedVITDRWqWLU6Ra+1w1W9S0BGj32SkReNZSFAz/HAGYVJhrJYDB0XnISjKnnKnmvjcWEGoqBapDIq4anJ6RvVPHStQdv/nrJSc0eeyUirxrKwoGfYwCzCiuNdCvr2lj7rz6ma9LBihrKgGqQyKjG5cvG3LkhgHt+OVqb170yR0Y1FIcDP8cAZhWGGkm4LUyzKYcSU+W6UyFDaigAqkEilxoZGcYmTSB9D9d9t5JKNwFzArnUUAMO/BwDmFUYaqSkVD2kL5jU5L3X6ZpEMKSGAqAaJHKpsXgxpG9KgUKN+61qPvUQrOT0DJpELjXUgAM/xwBmFbYa6eD1JxDA1YbvufFYFvHZUkNuUA0SWdSAV6tQAQL4r096w4p96Ia27jloA1nUUAkO/BwDmFWYa6RfVp8Dq+qw8IQ+Q/oDVZhTQ1ZQDRKZ1NDfD9vU6rsqQ3fCik3XNIxMaqgCB36OAcwqzDXS49jkN8b4QAavOHGXrrkMc2rICqpBIpMay4/fhZUZVmlYsemahpFJDVXgwM8xgFmFxUZafSoEPKv26H0Pnkp8kwYW1ZAPVINEYjUSEoy+vmExibAaw8q85lQIPYO2kVgNVeHAzzGAWYXFRsrIeH6Thh9WnJH2jEkW1ZAPVINEYjX++MOo0+1tbbrsxld/+2fI8H2KrEishqpw4OcYwKzCaCMFh8dVH7EXzGvnpYd0zQUYVUMmUA0SKdW4eNGYKxcE8A9fjYXVODg8np5B80iphtpw4OcYwKzCbiPN9Q2CAH5r/IGYBMmumsuuGnKAapBIpoZeb2zcGNL3QN33YQWedzCInoEFJFNDA3Dg5xjArMJuI6WmZ3wy23TfmN83XqJrzsKuGnKAapBIpsbChZC+SQUKv93vP1iBFbjHlxxIpoYG4MDPMYBZhelGCgiJEe4WfPhmOF1zCqbVkBxUg0QaNR4+NL78MgTwyE/6wqp7PjSGnoERpFFDG3Dg5xjArMJ6I43ffQ0C+J1JB58lS/ARWFdDWlANEmnU6NQJ0vdKhVpVhu6csPsaXWUHadTQBhz4OQYwq7DeSEmp+g+mH4YMHrY5kK45DutqSAuqQSKNGlu3Rpcs2/rHBR/O8GPlqpMWkUYNbcCBn2MAswoHjXTmbrSwI/rIrQi65iAcqCEhqAaJJGr43QyvOnQnrK5n70XTNaaQRA2NwIGfYwCzCh+N9NeuqxDATSa7uiOaDzWkAtUgcVWN5GRYOWEVhRV13C6Gdz4LuKqGluDAzzGAWYWPRkpMTX8/c0e051aXdkTzoYZUoBokLqlx4YKxVKlN/cZVGrYbVlSmdz4LuKSGxuDAzzGAWYWbRjp9JwoC2MUd0dyoIQmoBonzauj1xkaNjDrdrlrvVfb0PnOX7Z3PAs6roT048HMMYFbhqZHG7jTtiG48yfdpopOX5uBJDddBNUicV2PePEjfuPyFGvVfA6soXWUT59XQHhz4OQYwq/DUSEmp+o9m+kEG/+p1ga7ZB09quA6qQeKkGvfvGwsXhgAe3ro/rJwc7HwWcFINTcKBn2MAswpPjQRcuv+06vA9Tl8jmjM1XATVIHFSjY4dIX3PvVa7muduWDnpKrM4qYYm4cDPMYBZhadGEph94BYEcL2xztxglT81XAHVIHFGje3bIX3TXsr1yU8L5/jeoqss44waWoUDP8cAZhWeGkkgTZ/RfsFxyOBv/z3t6M0K+VPDFVANEifUMIwZk+HhsbBJ588XHIfVki6zjBNqaBYO/BwDmFV4aiSR4PD4GiNNNyv87+Q9umYTLtVwGlSDxAk1Vvnf6/DdrPp/br8dwd4NB23jhBqahQM/xwBmFZ4aiWTlibsQwBDDQU/i6Jp1eFXDOVANEkfVuPUkTtgKhBima+zjqBpahgM/xwBmFZ4aiSQjw/Dd8jNgf63nHE1Os/fQU17VcA5Ug8QBNdLT0/v2+3HEOlj9YCWEVZGegX0cUEPzcODnGMCswlMjUYTHJTeccABM0P6TLzlWwwlQDRIH1Jg1y6jTPSlc/J2xeyLiUugqFzighubhwM8xgFmFp0Yy5/DNcOHyWIduPKFrluBbDUdBNUjsVSMkRF+gAATwsDa/+Ul0m2oNYq8aLMCBn2MAswpPjWSRcbtMNwxuMP5A+LOcz0riXg2HQDVI7FLDYEhp0xbS93T5uuN3XqGrHGGXGozAgZ+zHcBjx47VEdSsWZOeIzscLDARnhrJIinp+jZzj0EGd192Osdv47hXwyFQDRJ71MjYtBnSN/Wl3L94roYVjy5zhD1qsAIHfs58ANetW/dxFpGRkfQc2eFggYnw1EjWCA6PqznKdDzqYr/bdC077qCG/aAaJFbV0OuNfn5GLy/j7t0JxUtCAC96t2twOG/nHVFYVYNBOPBz5gP4zTffpKdah4MFJsJTI9lg/ZlQCOCqw/fYvhG6m6hhJ6gGiWU1tm41li8PoSv+PC5UbOMxri56ZRHLarAJB37OfAAXLFiwbNmyVapU6datW2hoKD2H0ZiSkvIsi7CwMFhgUVFRaeyTmJgIjQT/0gW+SE1N/W3decjgxhN9Hz9NoMtZuIkadoJqkJirkb5xo8HDw0CkryHzJ23jRuL3+MRcDXYBJ8cAVpO9e/du2rQpMDDQx8enadOmFStWjIujr95AfU8MeHl57UDYYePWHY3/2m06M3jyrm3b6SqCOMbWrUmvvkqmr5jBiSVKQJWeH9Eq4OQYwFrh6dOnRYoU+ffff6npOALmgKth0cKXwXMP3KRrmbiVGjmCapBQaqT7+lLRS/5ANftv8wZP6waOgLVFo0aNPD096akEz9j/zkAkjaPvcuxh47n7EMBVPL1P3o6ia+6nhm1QDRJaDS8v89x98QNVrqHVYBkO/JyfAI6Pjy9WrNi8efPoAgEHC0yEp0ayk983XoIMbjTR1/zMYDdUwwaoBgmthp8fHbrkD1S5hlaDZTjwc7YDeMiQIUeOHLl3756/v3/Lli1LlCgRERFBz0TAwQIT4amR7CQxNb3V7KOQwZ0W+6emZ7tJnBuqYQNUg4RWQ683lCpl/h2w0cPDWKGC6dwkrqHVYBkO/JztAO7SpUvZsmXz5s372muvwePbt3M4W5SDBSbCUyPZz93IhDfG+EAGj9mR7XJF7qmGNVANElqN4OCkgoWNmUddZUtf+Nm6Ndtv8gitBstw4OdsB7CjcLDARHhqJIc4eP2JcJnoLQFh4kS3VcMiqAZJNjWePIl/rSIkbkjR0smlyrwIYBj7ukH6GvlaNzjwcwxgVuGpkRxl9oFblTLvGXzlQawwxZ3VMAfVICHVCPY5Gl2wSMgrZZZu8n9xJSz4l/c9zyI8rRsc+DkGMKvw1EiOkpFh+GGF6Z7BzaYciklINbq3GuagGiSiGtEJqbDCfNhrqeeULTleXZxXeFo3OPBzDGBW4amRnCA2Me396Ychg7v+cypNn+HmalCgGiQmNbZvT7p2DVYVWGE+mH44Nsl9leFp3eDAzzGAWYWnRnKOm4/j6ozeB5Y6avsVVIME1SABHW516pSSv8B3ncfBCgOrDT2HO8HTusGBn2MAswpPjeQ0B649qexpOiBr5fHbqIYIrhsk+jlzhCOt/mg7EFYYuuxm8LRucODnGMCswlMjucIiv2DhdknT/tuJagjguvGCTZsMHh6QvtPf+y7Hm1q6AzytGxz4OQYwq/DUSK5gMBgGb7gIGVx7xO6gx0/psluC68Zzjhwx5M0L6ftfg88GeQXAqkLP4H7wtG5w4OcYwKzCUyO5SHKavsPC45DBH804HJuIguC6kcnly4aiRSF999Zo9tG4HfGJ9OVL3ROe1g0O/BwDmFV4aiTXeRgd/+Zo0y0Lv15yMiXdXc7ptAauG4D+x58gfc+Ur/PBuD2rN7m7GiI8rRsc+DkGMKvw1EiuAzosWrejzhjTQdG/eV1w27M8BXDdgBVgwOozc5t1bTZ085X70W6uBglP6wYHfo4BzCo8NZLrCGocvv642vA9kMHT9t2g53An3HrdSE01GgxT992A1QBWhuNBkW6thhk8qcGBn2MAswpPjeQ6ohqbMm8bDD9rT4fQM7kN7rtu6PXGjh2DOnStOnQnrAObM68W7r5qWIInNTjwcwxgVuGpkVyHVGOOr+lK0VU8vQ/dcNOTPt103TAYjH36GHW6lFy52/w4f65vkDDZTdWwAk9qcODnGMCswlMjuQ6phsFgGLLpEmRwzVF7A0Ki6VndADddNyZMgPTN0Hn06eD5x6ZL4klHbqqGFXhSgwM/xwBmFZ4ayXUoNdL0GT0y79ZQb6zPjcc8LG6HcMd1Y/ly4XJXY1r2hkUPK4BYcUc1rMOTGhz4OQYwq/DUSK5jrkZSqv7Lxf6QwW9P9A2NSiTm5R9zNTjH29uQKxek7+J3vuq02B8WPVl0OzVswpMaHPg5BjCr8NRIrmNRjdjEtNZzjkIGvz/9cHicG12HwaIa3BITk1GkCKTvlroftZ59xPxKLO6lRk7wpAYHfo4BzCo8NZLrWFMj/Fnyu9MOQQZDEptbM69YU4NLYLGO/nnqvhpNP5x8ABY3XXYzNXKEJzU48HMMYFbhqZFcx4YaIVEJDSf4QgZ/vuD4s2QLM/CHDTU4AxYoLFZYuI2sf9HgPmrYA09qcODnGMCswlMjuY5tNW48fva/cfvBpr9c7B+fkk6XucO2GpwQF5fW6as+f22ExQoL18ahdm6hht3wpAYHfo4BzCo8NZLr5KjGlQex9cb6gFl3XnIyMZXzDM5RDeZJTdW3aGnU6W6UqPTmmH1XH8bSMxDwr4Yj8KQGB36OAcwqPDWS69ijxqX7T98YY8rgrv+cSk7j+YYN9qjBMBkZ6d26Q/om5snXpee8wLAc7kHJuRoOwpMaHPg5BjCr8NRIrmOnGgEhMXVGm27Y8O2/p6mTVXjCTjUYJW3IH5C+6R4v/dJ1/PnQGLpsBt9qOApPanDg5xjArMJTI7mO/WqcuRtda5Qpg7ssPcnr98H2q8EGer3Rz8/o5QX/pkyfIVxww7P977Ao6TktwZsarsGTGhz4OQYwq/DUSK7jkBpn70XXzdwX/cWiE1weF+2QGlpn61Zj+fJC6MKPIfPfuR/9cO6eXelr5EwNl+FJDQ78HAOYVXhqJNdxVI1L95/W/8t0XHS7+cdjElLpMuM4qoZ2gfT18BDTVwhg+Lm35D96Tuvwo4YU8KQGB36OAcwqPDWS6zihxrWHzxqMPyBco4Oz62Q5oYYW0evJse+LDIZIrlDBVLUPTtSQCJ7U4MDPMYBZhadGch3n1Ah6Evf2RNM1Ot6bdvheZAJdZhbn1NAcfn7m6fviB6r2wYkaEsGTGhz4OQYwq/DUSK7jtBqQu5C+kMENJxy4HGbrdFKGcFoNbeHlRYcu+QNV++BEDYngSQ0O/BwDmFV4aiTXcUWN8LjkT+cegwyuM3rfsaAIuswgrqihIQ4epEOX/MERsFPwpAYHfo4BzCo8NZLruKhGXHJa139OQQa/PmLPjosP6DJruKiGJoiMDH/nPTp0hR/8DtgFeFKDAz/HAGYVnhrJdVxXIyVd33/dechg+JnrG2QwGOg52MF1NdTFcPbss1LlIGtTcuU2HfZMHggNj+Fn61b6d6zDuhrSwpMaHPg5BjCr8NRIriOJGhkZhgm7rwkZPGD9BXYvVymJGmoBsi8fOB2y9m6xsksXbs/YvCXbsdAw9nUkfY2MqyE5PKnBgZ9jALMKT43kOhKq4XUmtNrwPcJlOiLjU+gyC0iohsJExKWA7CD+kHa/bz505flU4kpY9u95FmFXDTngSQ0O/BwDmFV4aiTXkVYN/+BI4dZJzaYcuv6IvbVFWjUUIiQkrmXrz4ebbi8I4sMioGdwFibVkA2e1ODAzzGAWYWnRnIdydW4HRH///buPSiq644DOIYIziTVzVhJBijgA2M7PproMOnojG2TqQ9S/KPJ1DqKxlSjwvhIwmiikaZGMOZhLFZHk4AxFonG+kiE4GN8pCJIMioIPkCJoMkuogICu1fZpQeurIeDGte7u/ee3/1+5o6z7D2m9rv3nO9e7j5GLm99e9KTi7Kle1mW19Pwud27HZbHWgICdveL+f17+89VXxcHaCBfGr5EKQ0C6zkKWFaUJpJ2vkjjWqMy6dMC9ZJw8o6Tyk2nOMKofJGGrzidN5e8o77M6vgT0XPf3c5iF8doI1MavkcpDQLrOQpYVpQmknY+SqPZ6fog97TawX9ZfdhaJ8cnVvooDe+7ds0+9nn11VWZg//00VcnWODiGM2kScMvKKVBYD1HAcuK0kTSzqdp7CmxDmy7JDx0yZ6DZyT4pA6fpuE15eWNkb3b3mvUdfGf5+wttYoDvESONPyFUhoE1nMUsKwoTSTtfJ1GxeWGUSsOqqfC73xdYvBfR/s6De0cN5uXbT56+pcRVd17Jc5b69MP4jZ+Gv5EKQ0C6zkKWFaUJpJ2fkjDfqP5re3FagfH/utQuVdfKORdfkjDM/z7iOz2clv92JWtn/054pVPlm045Ov3WxsuDV1RSoPAeo4ClhWliaSd39LYXWL97dutXyQ8YFHOxvwfjPmBWX5L475s3cp/kkZz167bB/6RBchiZGGKg33AWGnojVIaBNZzFLCsKE0k7fyZxk+19vFrWz84mm0TPj5SeaVRHKE3f6bxM1j78h8k2ba5AgL+/dJiFqM42DcMlIYBUEqDwHqOApYVpYmknZ/TaHa6Pj507slF2ayDf/1WzoYjPzh98PLdB+bnNO6qubnDp0hyBezy5NsUNDJKGsZAKQ0C6zkKWFaUJpJ2uqRx/nLDi2vy1FPhv67NK7MZ5aqwLmncwf79ndv39nbf3yeokVHSMAZKaRBYz1HAsqI0kbTTKw124pvxv/MDFuVEtn2VYWr2qQbHTXGQ3+mVhsD+2edi6fJbZqb4F3zDIGkYBKU0CKznKGBZUZpI2umbRuWVxqkZR9VT4WdS9u4q+lHfF2fpmwbjstlYCDP//r5YuvyGM2A9UEqDwHqOApYVpYmknRHS2FNiHfHuPrWGx689cqLqmjjCX/RMo7bWmviqo2vwpBff7p20w9ajV4dv81U3dg+uAeuEUhoE1nMUsKwoTSTtDJKG/Ubzh7vPRC9sfXEW22b95/sfanz4ERN34480On9FoKLUpLxX/wuL2rKbnhrDolC+2NJat3wHqz96+J2+WvgjDXlQSoPAeo4ClhWliaSdodKouto4L+tY1ILWDu77xq7F24v99pYblc/T6PjWXnb7euKcmid+pf5Y1jP8k6QVVVfan3kIg9m5rx/bt8UPaUiFUhoE1nMUsKwoTSTtDJhGyaW6+PYvU4p+M/uN/xb57R3Dvk2j01t7Xe03qh+xpMcvKLlwRfwrnU+X/ci3aciGUhoE1nMUsKwoTSTtDJvG4bLL7rcq9XljFzszPmOtFwd5mw/TuPtbexu7PXLk+3JxvAH4MA0JUUqDwHqOApYVpYmkncHTKDh/xf3Vwmz727ojOcU/3Wz21Tc6eCeNO562GuOtvR7xThpUUEqDwHqOApYVpYmknRRpHK+89sqG73q3XRtm2+9S9qbtO/tjbZM4TjMvpNHpKm/Lli01O3MuDnxaLF1+89dbez3ihTQIoZQGgfUcBSwrShNJO4nSuHitafk3p5765261hqMWtJ4Qby6srLd77R/vcRrCye6drvK6L/Tea8MZsOFRSoPAeo4ClhWliaSddGnYbzRv/b7KfXmYbU8uyp658Tt259UGRRztoRsNDUVTpjjj4lqee65l/PiWBQtaPvigZcOGlhUrWjZubK1JRbnduJs3iye7PXuKzdq2OQO67B4eZ3+sp+5v7fWIdMeGT1FKg8B6Ln0Br1q1KjIyMjg4OCYmpqCgQNzdEYEHzI3SRNJO3jQqrzSm7Tv7h/f3u5u494KvX1hzeM2B8hNV1x7kOnFS0h0KUtgCA8V72refOdNlha2eH+v61l6PyHts+AKlNAis53IXcFZWVlBQUHp6eklJybRp0ywWi81mEwdxCDxgbpQmknayp+FyuYqqat/PPT36o9Zvqndvv3krZ+In+ayhD5ddvq8z46Skny/RjptHg29d5dX7rb0ekf3Y8C5KaRBYz+UuYHbWm5CQoN52Op2hoaGpqakdh3RA4AFzozSRtKOURtXVxs/yKl7KODow+Ru+jNkWs3TP5PSClOzSzIILB85Ul9muNyncb30V5R6ntt7Z3Fd57/gaaUOidGxoRykNAuu5xAWsKEpgYOC2bdvc98THx8fFxXFDWjkcjrp2VVVV7AGrqam5Ib/GxkY2kdif4g5TIpmG3aGcuHDl00PlMzYUjlh261OmO2+Dkr8Zvmzf6BUHM16YLfal9zZXly6u8PAbdrv4rzQ8ksfGA6OUBlvJUcC6uXTpEks/Ly/PfU9SUhI7J+aGtEpOTg7oKDMzczuAbLK2bl/x+faktTsmfrhzTOrOZ/7xVf83v+KbOOPp5zsX5wNv/K+m1VdBF8yfL/6bAPTDVnIUsG7us4BxBmwG5kxDUZTLdY0lF68WlFfvLfmx6NXFnXv0AbbWk92ePV2hobfvCQ+/+cUX4v+8JMx5bNwNpTRwBqyn+/wVNK9O/msGbjcIXcvRDmm08so1YPdLmuW5yntvODZ4lNIgsJ5LXMAtbS/CSkxMVG87nc6wsDC8CMuckMYtnr8KWtyM/ZLmB4Bjg0cpDQLrudwFnJWVFRwcvH79+tLS0unTp1ssFqvVKg7iEHjA3ChNJO2Qxm2evg+YNe6WLTROdu8IxwaPUhoE1nO5C5hJS0uLiIgICgpiZ8P5+fni7o4IPGBulCaSdkiD59knYZFrXAGODR6lNAis59IXsEcIPGBulCaSdkiDhzR4SINHKQ0C6zkKWFaUJpJ2SIOHNHhIg0cpDQLrOQpYVpQmknZIg4c0eEiDRykNAus5ClhWlCaSdkiDhzR4SINHKQ0C6zkKWFaUJpJ2SIOHNHhIg0cpDQLrOQpYVpQmknZIg4c0eEiDRykNAus5ClhWlCaSdkiDhzR4SINHKQ0C6zkKWFaUJpJ2SIOHNHhIg0cpDQLrOQpYVpQmknZIg4c0eEiDRykNAus5ClhWlCaSdkiDhzR4SINHKQ0C67m5Cri2tpY9YFVVVXXyq6mpyczMZH+KO0wJafCQBg9p8CiloX69LFvVxYVeHuYqYPUBAwAAGtiqLi708jBXATudTvZosWdM4lMpCalPJmiczWuHNHhIg4c0eJTSYCs5+z/CVnVxoZeHuQqYkjr5r394EdLgIQ0e0uAhDUNBAcsKE4mHNHhIg4c0eEjDUFDAssJE4iENHtLgIQ0e0jAUFLCsHA5HcnIy+1PcYUpIg4c0eEiDhzQMBQUMAACgAxQwAACADlDAAAAAOkABAwAA6AAFDAAAoAMUMBEOh2PIkCEBAQHHjh0T95lGRUXF1KlTo6KiunXr1qdPn8WLFyuKIg6ibtWqVZGRkcHBwTExMQUFBeJuM0lJSRk2bNijjz7aq1evcePGnT59WhxhVqmpqWytmDNnjrgD/AsFTMTs2bPHjBlj8gLOycmZMmVKbm7uuXPnduzYERIS8tprr4mDSMvKygoKCkpPTy8pKZk2bZrFYrHZbOIg0xg1alRGRsbJkyePHz8+duzYiIiIhoYGcZD5HD16lD1JHTx4MApYdyhgCrKzswcMGMDWXJMXsGD58uW9e/cW7yWNnfUmJCSot51OZ2hoKDvX6TjEpKqrq9nsOHjwoLjDZK5fvx4dHb1nz56RI0eigHWHApae1WoNCwsrLCysqKhAAfMWLlw4dOhQ8V66FEUJDAzctm2b+574+Pi4uDhuiHmVlZWx2VFcXCzuMBl2SMydO5fdQAEbAQpYbi6Xa/To0UuWLGlpuwKKAnZjC2737t3XrVsn7qDr0qVL7ADIy8tz35OUlMTOibkhJuV0OmNjY4cPHy7uMJlNmzYNHDjQbre3oICNAQVsdPPnzw+4i1OnTq1cuZItK83NzS3mKOB7p+EedvHixb59+7788svcX6UPBXw3M2bMiIyMlPqLY7WrrKwMCQk5ceKE+iMK2AhQwEZXXV196i4URRk3btxDDz0U2I6tv+zP+Ph48b9Cxb3TUMewHoqOjp40aZLUXxT6APAr6DtKSEgIDw8/f/68uMNk2IGhrg/utaJLly7shvr0HXSBApbbhQsXitvl5uaySfXll1+a+Zk+O/dl7Tt+/HhzLivsfDcxMVG9zZ5/hIWFmflFWC6Xi7VvaGjo2bNnxX3mU19f714rmGHDhk2cOLHY9BfF9YUCpsMMv4K+N9a+/fr1e/bZZ9mNn9qJg0jLysoKDg5ev359aWnp9OnTLRaL1WoVB5nGzJkze/ToceDAAffB0NTUJA4yK/wK2ghQwHSggDMyMriLwreIg6hLS0uLiIgICgpiZ8P5+fnibjMRD4WAAHaEiIPMCgVsBKZbngAAAIwABQwAAKADFDAAAIAOUMAAAAA6QAEDAADoAAUMAACgAxQwAACADlDAAAAAOkABAwAA6AAFDAAAoAMUMAAAgA5QwABEVFdXP/7440uXLlV/PHz4cNeuXffu3dtxFAAYBQoYgI5du3ax0i0sLKyvr+/Tp8+8efPEEQBgGChgAFJmzZrVv3//CRMmDBo0yOFwiLsBwDBQwACkNDU1sXNfdh5cVFQk7gMAI0EBA5BSXFzcrVu3wMDAnTt3ivsAwEhQwAB0KIoyZMiQyZMnp6SkhISE2Gw2cQQAGAYKGICO119/PSoqqq6uzul0jhgxIjY2VhwBAIaBAgYgYv/+/Q8//PC3336r/lhRUdG9e/fVq1d3HAUARoECBgAA0AEKGAAAQAcoYAAAAB2ggAEAAHSAAgYAANABChgAAEAHKGAAAAAdoIABAAB0gAIGAADQAQoYAABAByhgAAAAHaCAAQAAdIACBgAA0MH/ASLBv2WdbDF9AAAAAElFTkSuQmCC" + } + }, "cell_type": "markdown", "id": "c40d45ab", "metadata": {}, "source": [ "# Gradient Based Optimization\n", "\n", - "## Problem description:\n", - "We have a TLE at a given time, which we call TLE$_{0}$, and we look for a TLE at a future observation time ($t_{obs}$): TLE$_{t}$. \n", + "## How it Works:\n", + "We start with a Two Line Element (TLE) at a given time, which we call TLE$_{0}$, and we want to use it to predict the actual position of a satellite at some future observation time ($t_{obs}$).\n", "\n", - "We can propagate the state from $t_0 \\rightarrow t_{obs}$, and obtain the state at $t_{obs}$. In general, we define the state (i.e., position and velocity), as:\n", + "We can use the SGP4 algorithm to propagate the state from $t_0 \\rightarrow t_{obs}$, and obtain the state at $t_{obs}$. We define the state (i.e., position and velocity) of an object in 3D space, as (this is described in more detail below):\n", "\n", "\\begin{equation}\n", "\\vec{x}(t)=[x(t), y(t), z(t), \\dot{x}(t), \\dot{y}(t), \\dot{z}(t)]^T\n", "\\end{equation}\n", "\n", - "We then have: TLE$_0$, $\\vec{x}(t_0)$, and $\\vec{x}(t_{obs})$, but we want to find TLE$_{obs}$. That is, the TLE at the observation time, that when propagated with SGP4 at its time, it corresponds to that $\\vec{x}(t_{obs})$. In general, this means that we are able to invert from the state to the TLE, at any given time. \n", + "where:\n", + "\n", + "- **$\\vec{x}(t)$:** This is the state vector of the object at time $t$. The vector encapsulates all the necessary information to describe the object's state in space at some specific moment. In our case, the object is a TLE.\n", + "- **$[ \\cdot ]$:** Indicates this is a vector\n", + "- **$x(t), y(t), z(t)$:** These are the positional coordinates of the object at time $t$. They represent the location of the object in a three-dimensional space\n", + "- **$\\dot{x}(t), \\dot{y}(t), \\dot{z}(t)$:** These terms represent the velocities of the object in the direction of each corresponding axis (X, Y, and Z). The dot above each symbol signifies these are first derivatives of the position coordinates. $\\dot{x}(t)$ is the velocity in the X direction, $\\dot{y}(t)$ is the velocity in the Y direction, and $\\dot{z}(t)$ is the velocity in the Z direction.\n", + "\n", + "We then have: TLE$_0$, $\\vec{x}(t_0)$, and $\\vec{x}(t_{obs})$, but we want to find TLE$_{obs}$. Said with an example: Imagine we have some satellite with TLE$_0$. For TLE$_0$ we have all the orbital parameters required to define the satellite's orbit at $t_0$; inclination, right ascension of the ascending node, eccentricity, etc. Now we use TLE$_0$ as input to an orbital propagation model (SGP4) to compute the satellite's state vector at time 0 ($t_0$). This calculation gives us $\\vec{x}(t_0)$ which is the satellite's position and velocity in space. Now fast-forward a bit into the future to a new time ($t_{obs}$) where we again observe the satellite's actual position and velocity ($\\vec{x}(t_{obs})$). The problem is that our initial estimate $\\vec{x}(t_0) \\neq \\vec{x}(t_{obs})$. Our estimate was wrong. What we now want to do is to calculate TLE${_{obs}}$. TLE${_{obs}}$ is the TLE set that when applied at our new observation time **would** give us $\\vec{x}(t_{obs})$. Said another way, we want to be able to invert from state to an accurate TLE that correctly produces TLE$_{obs}$.\n", + "\n", + "This is where gradient descent comes into play. If you have no prior experience with gradient descent this is very likely going to be confusing. You can read more about it [here](https://www.khanacademy.org/math/multivariable-calculus/applications-of-multivariable-derivatives/optimizing-multivariable-functions/a/what-is-gradient-descent). The very high level description is that gradient descent is an optimization algorithm used to minimize a function by iteratively moving in the direction of the steepest decrease, as defined by the negative of the gradient. We use it in machine learning to iteratively refine model parameters, aiming to find the set of parameters that minimizes the cost function, typically representing the discrepancy between predicted and observed data. \n", + "\n", + "So how does gradient descent apply here? We can reformulate our problem, that is, how do we take $\\vec{x}(t_0)$ (that's the estimation made by SGP4) and get it so that it is as close as we can to $\\vec{x}(t_{obs})$ (the actual position of an object at time $t_{obs}$). This is where gradient descent comes in. We create a cost function:\n", + "\n", + "\\begin{equation}\n", + "F(\\vec{y}) = \\left| \\text{SGP4}\\left(\\text{TLE}(\\vec{y}), t_{\\text{obs}}\\right) - \\vec{x}(t_{\\text{obs}}) \\right|\n", + "\\end{equation}\n", + "\n", + "1. **SGP4(TLE($ \\vec{y} $), $ t_{obs} $)**: The satellite state predicted at the observation time $ t_{obs} $ by the SGP4 model, using the TLE parameters adjusted by the vector $ \\vec{y} $.\n", + "2. **$ \\vec{x}(t_{obs}) $**: The actual observed satellite state at the observation time $ t_{obs} $.\n", "\n", - "In order to do this, we formulate the problem as looking for the minimum of a function of a free variables vector (i.e., $\\vec{y}$) $F(\\vec{y})$, where this function defines the difference between the given state propagated from TLE$_0$ at $t_{obs}$, and the state generated from the free variables that make a TLE which is then propagated at its current time: TLE$(\\vec{y})(t_{0}\\rightarrow t_{obs})$. So we can reformulate the problem as:\n", + "Our goal with this equation is to make $F(\\vec{y})$ as small as possible, that is, we want to minimize the function. Before I explain why, let me explain what $\\vec{y}$ is. $\\vec{y}$ is a series of parameters that are going to modify the output of our SGP4 algorithm. The way we're going to do this is through gradient descent. Imagine it like this: SGP4 includes all sorts of things eccentricity, drag coefficient, solar radiation pressure, lunar perturbations, and a whole bunch of other stuff but at the very end it is still wrong by a little bit. So you can think of what we're doing here is adding a parameter in front of all this stuff (it doesn't work exactly like this, but this is the intuition) and then we're going to use gradient descent to figure out which combination of parameters modifies SGP4 from guessing incorrectly, to guessing correctly across a broad dataset. Those are our $\\vec{y}$. Our minimization function's goal is to figure out which $\\vec{y}$ most reduces the difference between $\\text{SGP4}\\left(\\text{TLE}(\\vec{y}), t_{\\text{obs}}\\right)$ and $\\vec{x}(t_{\\text{obs}})$. Expressed mathematically:\n", "\n", "\\begin{align}\n", "\\textrm{given}: & \\ \\textrm{TLE}_0, \\vec{x}_0\\\\\n", @@ -26,31 +49,61 @@ "\\textrm{that minimize}: & F(\\vec{y})=|SGP4(\\textrm{TLE}(\\vec{y}),t_{obs})-\\vec{x}(t_{obs})| =|\\vec{\\tilde{x}}(t_{obs})-\\vec{x}(t_{obs})|\n", "\\end{align}\n", "\n", - "We can do this via Newton method, by updating an initial guess $y_{0}$ until convergence. Where the update is done as follows:\n", + "We can do this via [Newton's Method](https://math.libretexts.org/Bookshelves/Calculus/Calculus_(OpenStax)/04%3A_Applications_of_Derivatives/4.09%3A_Newtons_Method). We do not explain Newton's Method here, however the intuition and what is important here is that it provides a method to take some initial parameter estimation $y_{0}$ and minimize it such that the parameter $y$ modifies function $\\vec{x}$ to more closely approximate $\\vec{x}(t_{obs})$. When we run gradient descent, you can imagine this visually as us trying to find the bottom of a parabola as that is when some function is at its lowest value. Recall, that's what we're trying to do here, find the lowest possible value of a function, where the function is minimizing the difference between where we guess a satellite is and where it really is. Smaller values means our guess is more accurate. This is what it looks like:\n", + "\n", + "![gradient_descent.png](attachment:1f7a6234-cf18-4c3d-a0e4-cc66a5d5b53a.png)\n", + "\n", + "This picture also illustrates Newton's Method graphically. We start with some guess and we follow tangent lines until we find where $f(x)=0$. Here is what Newton's Method looks like expressed mathematically:\n", "\n", "\\begin{equation}\n", "y_{k+1}=y_{k}-DF^{-1}(y_k)F(y_k)\n", "\\end{equation}\n", "\n", - "with $DF$ the Jacobian of $F$ with respect to $y_k$. We can easily see that this Jacobian is made of the following elements:\n", + "with $DF$ the [Jacobian](https://www.khanacademy.org/math/multivariable-calculus/multivariable-derivatives/jacobian/v/the-jacobian-matrix) of $F$ with respect to $y_k$. Conceptually, this is a bit heady if you haven't seen it before. To wrap your head around this, you will need a good understanding of [partial derivatives](https://www.khanacademy.org/math/multivariable-calculus/multivariable-derivatives/partial-derivatives/v/partial-derivatives-introduction), but I will explain the intuition. A partial derivative gives us information about how the value of a function changes if we hold constant all variables except one. That is to say, if we want to go \"down\" it tells us which direction has a steeper downward gradient. This is relevant to gradient descent because we want to find the point at which the function is lowest and while this is simple when the graph is 2D, this becomes harder when the graph is 3D, and much hard when the graph is many dimensions AKA N-dimensional. The Jacobian packs all this partial information into a single matrix that allows us to plug in values for efficient computation. It looks like this:\n", + " \n", + "$$\n", + "DF = \n", + "\\begin{bmatrix}\n", + "\\frac{\\partial F_1}{\\partial y_1} & \\frac{\\partial F_1}{\\partial y_2} & \\cdots & \\frac{\\partial F_1}{\\partial y_m} \\\\\n", + "\\frac{\\partial F_2}{\\partial y_1} & \\frac{\\partial F_2}{\\partial y_2} & \\cdots & \\frac{\\partial F_2}{\\partial y_m} \\\\\n", + "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n", + "\\frac{\\partial F_n}{\\partial y_1} & \\frac{\\partial F_n}{\\partial y_2} & \\cdots & \\frac{\\partial F_n}{\\partial y_m}\n", + "\\end{bmatrix}\n", + "$$\n", + "\n", + "The takeaway is that it gives us the ability to efficiently calculate if we are heading toward a global minima or not. In our case, the Jacobian has the following elements:\n", "\n", "\\begin{equation}\n", "DF_{ij}=\\dfrac{\\partial \\tilde{x}_{i}}{\\partial y_{j}}|_{y_k}\n", "\\end{equation}\n", "\n", - "where $\\tilde{x}_{i} \\in [\\tilde{x}_1,\\tilde{x}_2,\\tilde{x}_3,\\tilde{x}_4,\\tilde{x}_5,\\tilde{x}_6]=[\\tilde{x},\\tilde{y},\\tilde{z},\\tilde{\\dot{x}},\\tilde{\\dot{y}},\\tilde{\\dot{z}}]$; and $y_i \\in [no_{kozai}, ecco, inclo, mo, argpo, nodeo, n_{dot},n_{ddot},B^*]$.\n", + "where $\\tilde{x}_{i} \\in [\\tilde{x}_1,\\tilde{x}_2,\\tilde{x}_3,\\tilde{x}_4,\\tilde{x}_5,\\tilde{x}_6]=[\\tilde{x},\\tilde{y},\\tilde{z},\\tilde{\\dot{x}},\\tilde{\\dot{y}},\\tilde{\\dot{z}}]$ (our state vector); and $y_i \\in [no_{kozai}, ecco, inclo, mo, argpo, nodeo, n_{dot},n_{ddot},B^*]$ (these are the values we are going to update with gradient descent).\n", "\n", - "Since we built a differentiable SGP4, we can compute the gradient of the state w.r.t. the TLE inputs quite easily. Furthermore, the initial guess ($y_{0}$) will be found by the simple inversion from Cartesian to Keplerian, which does not correctly invert from state to TLE, but it is good as initial approximation.\n", + "- $no_{kozai}$: The mean motion of the satellite, adjusted for the Kozai correction. This represents the number of orbits the satellite completes in a day, adjusted for long-term perturbations in the orbit.\n", + "- $ecco$: The eccentricity of the orbit. This value defines the shape of the satellite's orbit, ranging from 0 (a perfect circle) to values close to 1 (highly elliptical orbits).\n", + "- $inclo$: The inclination of the orbit, measured in degrees. It indicates the angle between the satellite's orbital plane and the equatorial plane of the Earth.\n", + "- $mo$: The mean anomaly at the epoch. This is an angular measurement that specifies the satellite's position along its orbit at the specific time defined by the epoch of the TLE set.\n", + "- $argpo$: The argument of perigee. This angle indicates the orientation of the elliptical orbit in relation to the Earth's surface, specifying the point where the satellite passes closest to the Earth.\n", + "- $nodeo$: The right ascension of the ascending node (RAAN). This is the angle from a fixed reference direction, typically the vernal equinox, to the location where the satellite crosses the equatorial plane going northward.\n", + "- $n_{dot}$: The first derivative of the mean motion. It indicates how the satellite's mean motion changes over time, which is primarily due to atmospheric drag and gravitational perturbations.\n", + "- $n_{ddot}$: The second derivative of the mean motion. This value provides a refinement on the rate of change in the satellite's mean motion, offering a more precise prediction of its long-term orbital behavior.\n", + "- $B^*$: The ballistic coefficient. It relates to how the satellite responds to atmospheric drag, with a higher value indicating a greater effect of drag on the satellite's orbit.\n", "\n", + "Since we built a differentiable SGP4, we can compute the gradient of the state w.r.t. the TLE inputs quite easily. The initial estimate, denoted as $(y_{0})$, is obtained through a straightforward conversion from Cartesian coordinates (the state vector previously mentioned) to Keplerian elements (the orbital parameters depicting the satellite's state). While this conversion doesn't perfectly invert the state to the original TLE format, it serves as a suitable starting point for the approximation. Said another way, the output of our gradient descent process is a bunch of parameters that we use to refine our state vector $\\vec{x}(t_0)$, but what we really want is a TLE (Keplerian elements). Specifically, we want to take our new state vector and use it to find a better version of TLE$_0$ such that it accurately guesses $t_{obs}$.\n", "\n", - "\n" + "## Code Demo\n" ] }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 1, "id": "95627789", - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-29T15:46:00.981855Z", + "start_time": "2024-04-29T15:46:00.979281Z" + } + }, "outputs": [], "source": [ "import dsgp4\n", @@ -59,29 +112,40 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 2, "id": "91d2d6c5", - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-29T15:46:00.984906Z", + "start_time": "2024-04-29T15:46:00.982562Z" + } + }, "outputs": [], "source": [ - "#we initialize the TLEs\n", - "dsgp4.initialize_tle(tles);\n", - "#we extract one:\n", + "# Initialize the TLEs\n", + "dsgp4.initialize_tle(tles)\n", + "\n", + "# Extract the first TLE to show gradient-based optimization usage:\n", "my_tle=tles[0]" ] }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 3, "id": "c376c651", - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-29T15:46:01.019529Z", + "start_time": "2024-04-29T15:46:00.985418Z" + } + }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "F(y): 1.0751821836135118e-13\n", - "Solution found, at iter: 6\n", + "F(y): 6.569747093430235e-14\n", + "Solution found, at iter: 5\n", "TLE(\n", "0 COSMOS 2251 DEB\n", "1 34454U 93036SX 22068.91971155 .00000319 00000-0 11812-3 0 9996\n", @@ -94,7 +158,12 @@ } ], "source": [ - "#when I propagate to zero, I expect the returned TLE to be identical to my_tle:\n", + "# This line applies Newton's method to refine the initial TLE (my_tle) by attempting to minimize the difference \n", + "# between the propagated state and the expected state at the same time (time_mjd equal to my_tle's date). The goal is to \n", + "# adjust the TLE so that, when propagated to the given time (which, in this case, is the TLE's own timestamp), the resulting \n", + "# TLE should match the original TLE (my_tle) as closely as possible. The expectation is that, since we are propagating to the \n", + "# same time as the TLE's timestamp, the optimized TLE (found_tle) should be identical or very close to the original TLE (my_tle).\n", + "# The results are then printed to compare the original and the found TLEs.\n", "found_tle, y=dsgp4.newton_method(tle_0=my_tle,time_mjd=my_tle.date_mjd,new_tol=1e-12,max_iter=100)\n", "\n", "print(my_tle,found_tle)" @@ -102,7 +171,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 4, "id": "e0f85268", "metadata": {}, "outputs": [ @@ -111,7 +180,7 @@ "output_type": "stream", "text": [ "F(y): 4.547475677268469e-13\n", - "Solution found, at iter: 7\n", + "Solution found, at iter: 6\n", "TLE(\n", "0 COSMOS 2251 DEB\n", "1 34454U 93036SX 22068.91971155 .00000319 00000-0 11812-3 0 9996\n", @@ -124,7 +193,13 @@ } ], "source": [ - "#I now propagate until 1000 minutes after\n", + "# This code snippet uses Newton's method to propagate the TLE (my_tle) 800 days beyond the TLE's original epoch\n", + "# (my_tle.date_mjd). The aim is to adjust the TLE parameters so that the satellite's predicted state at this future\n", + "# date closely matches the expected observational data. The optimization employs Newton's method, configured to run\n", + "# for a maximum of 20 iterations or until the change between iterations falls below the threshold of 1e-12. The\n", + "# optimized TLE (found_tle) is expected to be different from the initial TLE (my_tle) due to the significant time\n", + "# interval between the original and target dates. The final part of the code prints both the original and the\n", + "# optimized TLEs to compare them.\n", "found_tle, y=dsgp4.newton_method(tle_0=my_tle,\n", " time_mjd=my_tle.date_mjd+800.,\n", " new_tol=1e-12,\n", @@ -132,6 +207,19 @@ "#Newton still converges, and the TLE is of course now different:\n", "print(my_tle, found_tle)" ] + }, + { + "cell_type": "markdown", + "id": "94d5d1894caf5076", + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "source": [ + "Explanation by Grant Curell" + ] } ], "metadata": { diff --git a/doc/notebooks/tle_object.ipynb b/doc/notebooks/tle_object.ipynb index 2f0e53f..6c8c92f 100644 --- a/doc/notebooks/tle_object.ipynb +++ b/doc/notebooks/tle_object.ipynb @@ -20,12 +20,17 @@ }, { "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-29T14:44:40.355175Z", + "start_time": "2024-04-29T14:44:39.530116Z" + } + }, "source": [ "import dsgp4" - ] + ], + "outputs": [], + "execution_count": 1 }, { "cell_type": "markdown", @@ -38,8 +43,23 @@ }, { "cell_type": "code", - "execution_count": 93, - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-29T14:44:40.362812Z", + "start_time": "2024-04-29T14:44:40.356711Z" + } + }, + "source": [ + "#let us assume we have the following two lines for the TLE, plus the first line that indicates the satellite name:\n", + "tle_lines = []\n", + "tle_lines.append('0 TIMATION 1')\n", + "tle_lines.append('1 2847U 67053E 24063.46171465 .00000366 00000-0 27411-3 0 9994')\n", + "tle_lines.append('2 2847 69.9643 216.8651 0003597 77.7866 282.3646 14.02285835897007')\n", + "\n", + "#let us construct the TLE object\n", + "tle=dsgp4.tle.TLE(tle_lines)\n", + "print(tle)" + ], "outputs": [ { "name": "stdout", @@ -53,17 +73,7 @@ ] } ], - "source": [ - "#let us assume we have the following two lines for the TLE, plus the first line that indicates the satellite name:\n", - "tle_lines = []\n", - "tle_lines.append('0 TIMATION 1')\n", - "tle_lines.append('1 2847U 67053E 24063.46171465 .00000366 00000-0 27411-3 0 9994')\n", - "tle_lines.append('2 2847 69.9643 216.8651 0003597 77.7866 282.3646 14.02285835897007')\n", - "\n", - "#let us construct the TLE object\n", - "tle=dsgp4.tle.TLE(tle_lines)\n", - "print(tle)" - ] + "execution_count": 2 }, { "cell_type": "markdown", @@ -74,8 +84,32 @@ }, { "cell_type": "code", - "execution_count": 43, - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-29T14:44:40.368784Z", + "start_time": "2024-04-29T14:44:40.364047Z" + } + }, + "source": [ + "#let's print all TLE elements:\n", + "print(\"TLE elements:\")\n", + "print(f\"Satellite catalog number: {tle.satellite_catalog_number}\")\n", + "print(f\"Classification: {tle.classification}\")\n", + "print(f\"International designator: {tle.international_designator}\")\n", + "print(f\"Epoch year: {tle.epoch_year}\")\n", + "print(f\"Epoch day: {tle.epoch_days}\")\n", + "print(f\"First time derivative of the mean motion: {tle._ndot}\")\n", + "print(f\"Second time derivative of the mean motion: {tle._nddot}\")\n", + "print(f\"BSTAR drag term: {tle._bstar}\")\n", + "print(f\"Element set number: {tle.element_number}\")\n", + "print(f\"Inclination [rad]: {tle._inclo}\")\n", + "print(f\"Right ascension of the ascending node [rad]: {tle._nodeo}\")\n", + "print(f\"Eccentricity [-]: {tle._ecco}\")\n", + "print(f\"Argument of perigee [rad]: {tle._argpo}\")\n", + "print(f\"Right ascension of ascending node [rad]: {tle._nodeo}\")\n", + "print(f\"Mean anomaly [rad]: {tle._mo}\")\n", + "print(f\"Mean motion [rad/min]: {tle._no_kozai}\")" + ], "outputs": [ { "name": "stdout", @@ -101,26 +135,7 @@ ] } ], - "source": [ - "#let's print all TLE elements:\n", - "print(\"TLE elements:\")\n", - "print(f\"Satellite catalog number: {tle.satellite_catalog_number}\")\n", - "print(f\"Classification: {tle.classification}\")\n", - "print(f\"International designator: {tle.international_designator}\")\n", - "print(f\"Epoch year: {tle.epoch_year}\")\n", - "print(f\"Epoch day: {tle.epoch_days}\")\n", - "print(f\"First time derivative of the mean motion: {tle._ndot}\")\n", - "print(f\"Second time derivative of the mean motion: {tle._nddot}\")\n", - "print(f\"BSTAR drag term: {tle._bstar}\")\n", - "print(f\"Element set number: {tle.element_number}\")\n", - "print(f\"Inclination [rad]: {tle._inclo}\")\n", - "print(f\"Right ascension of the ascending node [rad]: {tle._nodeo}\")\n", - "print(f\"Eccentricity [-]: {tle._ecco}\")\n", - "print(f\"Argument of perigee [rad]: {tle._argpo}\")\n", - "print(f\"Right ascension of ascending node [rad]: {tle._nodeo}\")\n", - "print(f\"Mean anomaly [rad]: {tle._mo}\")\n", - "print(f\"Mean motion [rad/min]: {tle._no_kozai}\")" - ] + "execution_count": 3 }, { "cell_type": "markdown", @@ -131,19 +146,12 @@ }, { "cell_type": "code", - "execution_count": 59, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Semi-major axis [km]: 7264.027802311157\n", - "Apogee radius [km]: 888.5036731116484\n", - "Perigee radius [km]: 883.2779315106654\n" - ] + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-29T14:44:40.374357Z", + "start_time": "2024-04-29T14:44:40.370687Z" } - ], + }, "source": [ "#let's first define the Earth radius according to WSG-84:\n", "r_earth=dsgp4.util.get_gravity_constants('wgs-84')[2].numpy()*1e3\n", @@ -154,7 +162,19 @@ "#let's extract the TLE apogee & perigee altitudes:\n", "print(f\"Apogee radius [km]: {tle.apogee_alt(r_earth).numpy()*1e-3}\")\n", "print(f\"Perigee radius [km]: {tle.perigee_alt(r_earth).numpy()*1e-3}\")" - ] + ], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Semi-major axis [km]: 7264.027802311157\n", + "Apogee radius [km]: 888.5036731116484\n", + "Perigee radius [km]: 883.2779315106654\n" + ] + } + ], + "execution_count": 4 }, { "cell_type": "markdown", @@ -167,9 +187,12 @@ }, { "cell_type": "code", - "execution_count": 89, - "metadata": {}, - "outputs": [], + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-29T14:44:40.378261Z", + "start_time": "2024-04-29T14:44:40.375410Z" + } + }, "source": [ "tle_dictionary=dict(\n", " satellite_catalog_number=43437,\n", @@ -189,12 +212,23 @@ " raan=4.3618,\n", " mean_anomaly=4.5224,\n", " b_star=0.0001)" - ] + ], + "outputs": [], + "execution_count": 5 }, { "cell_type": "code", - "execution_count": 92, - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-29T14:44:40.381821Z", + "start_time": "2024-04-29T14:44:40.379062Z" + } + }, + "source": [ + "#now the same API as above applies:\n", + "tle = dsgp4.tle.TLE(tle_dictionary)\n", + "print(tle)" + ], "outputs": [ { "name": "stdout", @@ -207,11 +241,7 @@ ] } ], - "source": [ - "#now the same API as above applies:\n", - "tle = dsgp4.tle.TLE(tle_dictionary)\n", - "print(tle)" - ] + "execution_count": 6 }, { "cell_type": "markdown", @@ -230,8 +260,16 @@ }, { "cell_type": "code", - "execution_count": 95, - "metadata": {}, + "metadata": { + "ExecuteTime": { + "end_time": "2024-04-29T14:44:40.385322Z", + "start_time": "2024-04-29T14:44:40.382613Z" + } + }, + "source": [ + "tles=dsgp4.tle.load('example.tle')\n", + "print(tles)" + ], "outputs": [ { "name": "stdout", @@ -253,10 +291,7 @@ ] } ], - "source": [ - "tles=dsgp4.tle.load('example.tle')\n", - "print(tles)" - ] + "execution_count": 7 } ], "metadata": { diff --git a/dsgp4/initl.py b/dsgp4/initl.py index ec3d147..29b79d0 100644 --- a/dsgp4/initl.py +++ b/dsgp4/initl.py @@ -4,59 +4,104 @@ from . import util def initl( - xke, j2, - ecco, epoch, inclo, no, - method, - opsmode, - ): - - x2o3 = torch.tensor(2.0 / 3.0); - - eccsq = ecco * ecco; - omeosq = 1.0 - eccsq; - rteosq = omeosq.sqrt(); - cosio = inclo.cos(); - cosio2 = cosio * cosio; - - ak = torch.pow(xke / no, x2o3); - d1 = 0.75 * j2 * (3.0 * cosio2 - 1.0) / (rteosq * omeosq); - del_ = d1 / (ak * ak); - adel = ak * (1.0 - del_ * del_ - del_ * - (1.0 / 3.0 + 134.0 * del_ * del_ / 81.0)); - del_ = d1/(adel * adel); - no = no / (1.0 + del_); - - ao = torch.pow(xke / no, x2o3); - sinio = inclo.sin(); - po = ao * omeosq; - con42 = 1.0 - 5.0 * cosio2; - con41 = -con42-cosio2-cosio2; - ainv = 1.0 / ao; - posq = po * po; - rp = ao * (1.0 - ecco); - method = 'n'; - - if opsmode == 'a': - # gst time - ts70 = epoch - 7305.0; - ds70 = torch.floor_divide(ts70 + 1.0e-8,1); - tfrac = ts70 - ds70; - # find greenwich location at epoch - c1 = torch.tensor(1.72027916940703639e-2); - thgr70= torch.tensor(1.7321343856509374); - fk5r = torch.tensor(5.07551419432269442e-15); - c1p2p = c1 + (2*numpy.pi); - gsto = (thgr70 + c1*ds70 + c1p2p*tfrac + ts70*ts70*fk5r) % (2*numpy.pi) - if gsto < 0.0: - gsto = gsto + (2*numpy.pi); - - else: - gsto = util.gstime(epoch + 2433281.5); - - return ( - no, - method, - ainv, ao, con41, con42, cosio, - cosio2,eccsq, omeosq, posq, - rp, rteosq,sinio , gsto, - ) + xke, j2, + ecco, epoch, inclo, no, + method, + opsmode, +): + # Initialize a tensor for the fraction 2/3, used in Kepler's third law calculations. + # Kepler third law is: The ratio of the square of an object's orbital period with the cube of the semi-major axis + # of its orbit is the same for all objects orbiting the same primary. + x2o3 = torch.tensor(2.0 / 3.0) + + # Calculate the square of the eccentricity to evaluate orbit shape and perturbation effects + # ecco = eccentricity + eccsq = ecco * ecco + + # Compute one minus the square of eccentricity. This is used to describe how much an object's orbit differs from + # a perfect circle. If the eccentricity were zero, the value would be 1, IE: a perfect circle. + omeosq = 1.0 - eccsq + + # Calculate the square root of (1 - eccsq), used in later orbital calculations. Ex: calculating orbital radius, + # and other orbital geometry + rteosq = omeosq.sqrt() + + # Compute the cosine of the inclination, determines the object's orientation which we can use later to calculate + # how the object will move along the orbit + cosio = inclo.cos() + + # Square the cosine of the inclination, used in perturbation calculations. + cosio2 = cosio * cosio + + # Compute the semi-major axis (ak - half the diameter of an ellipse between its two furthest points) from Kepler's + # third law, essential for orbit scaling. + ak = torch.pow(xke / no, x2o3) + + # Calculate the first part of the drag term, important for orbital decay predictions. The equation effectively + # calculates the drag induced by the "oblateness" of the earth. IE: the fact that it is flatter at the poles than + # in the middle. + # j2: represents the flattening of the Earth at the poles and its effect on the gravitational field. + d1 = 0.75 * j2 * (3.0 * cosio2 - 1.0) / (rteosq * omeosq) + + # Estimate the perturbation factor to adjust the semi-major axis. + del_ = d1 / (ak * ak) + + # Adjust the semi-major axis with the perturbation factor. This improves the orbit accuracy. + adel = ak * (1.0 - del_ * del_ - del_ * + (1.0 / 3.0 + 134.0 * del_ * del_ / 81.0)) + + # Recalculate the perturbation factor with the adjusted semi-major axis. + del_ = d1 / (adel * adel) + + # Correct the mean motion with the recalculated perturbation + no = no / (1.0 + del_) + + # Recompute the semi-major axis + ao = torch.pow(xke / no, x2o3) + + # Calculate the sine of the inclination, used in orientation and perturbation computations. + sinio = inclo.sin() + + # Compute the orbital period + po = ao * omeosq + + # Define constants related to the inclination, used in gravitational perturbation corrections. + con42 = 1.0 - 5.0 * cosio2 + con41 = -con42 - cosio2 - cosio2 + + # Calculate the inverse of the semi-major axis, important for orbit adjustments. + ainv = 1.0 / ao + + # Compute the square of the orbital period, used in timing and synchronization. + posq = po * po + + # Determine the perigee radius, essential for closest approach calculations. + rp = ao * (1.0 - ecco) + + # Set the method to 'n' - standard processing mode. + method = 'n' + + if opsmode == 'a': + # Perform calculations for the Greenwich Sidereal Time in alternative mode. + ts70 = epoch - 7305.0 + ds70 = torch.floor_divide(ts70 + 1.0e-8, 1) + tfrac = ts70 - ds70 + c1 = torch.tensor(1.72027916940703639e-2) + thgr70 = torch.tensor(1.7321343856509374) + fk5r = torch.tensor(5.07551419432269442e-15) + c1p2p = c1 + (2 * numpy.pi) + gsto = (thgr70 + c1 * ds70 + c1p2p * tfrac + ts70 * ts70 * fk5r) % (2 * numpy.pi) + if gsto < 0.0: + gsto += 2 * numpy.pi + else: + # Compute the Greenwich Sidereal Time for standard processing. + gsto = util.gstime(epoch + 2433281.5) + + # Return a tuple of updated orbital elements and computed variables + return ( + no, + method, + ainv, ao, con41, con42, cosio, + cosio2, eccsq, omeosq, posq, + rp, rteosq, sinio, gsto, + )