{ "cells": [ { "cell_type": "code", "execution_count": 121, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "5.0 6.0\n", "-0.0036120437044643043 0.23377764703610254\n", "-0.0036120437044643043 0.1467266512161534\n", "-0.0036120437044643043 0.07449997714053234\n", "-0.0036120437044643043 0.03572736472611628\n", "-0.0036120437044643043 0.016086059610160422\n", "-0.0036120437044643043 0.006239606566377719\n", "-0.0036120437044643043 0.001313915095452941\n", "-0.0011490925370628443 0.001313915095452941\n", "-0.0011490925370628443 8.24117820105128e-5\n", "-0.0005333412219636843 8.24117820105128e-5\n", "-0.0002254648275609566 8.24117820105128e-5\n", "-7.15264792828104e-5 8.24117820105128e-5\n", "-7.15264792828104e-5 5.442670007391032e-6\n", "-3.30419178462913e-5 5.442670007391032e-6\n", "-1.3799629748509031e-5 5.442670007391032e-6\n", "-4.1784867445559425e-6 5.442670007391032e-6\n", "-4.1784867445559425e-6 6.320562179281153e-7\n", "-1.773238489790817e-6 6.320562179281153e-7\n", "-5.705696283779814e-7 6.320562179281153e-7\n", "-5.705696283779814e-7 3.070150662287306e-8\n", "5.005728244781494\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3dZ2AU1d4G8P+Z2dRNL5CQQgohQEIChNACCAIKCIgKIlwUrw3BhtgQL4qAIoIFr4pcQRH1vnppKqIIiAYBAyEk9JYQUgiEkpBedmfO+2E1YgrZJLM7W57fp80yzJwdwjz7P3PmHMY5JwAAAHslqN0AAAAANSEIAQDAriEIAQDAriEIAQDAriEIAQDAriEIAQDAriEIAQDAriEIAQDAriEIAQDAriEIAQDArllWEHLOn3nmGeO31+l0pmsM3Jgsy5Ikqd0K+4VffhXh5KtIkiTFZwZlFjXXqF6vd3FxMf6XrLy83M3NzaRNgqbodDpZlp2cnNRuiJ0qKytzd3dXuxV2CidfRdXV1RqNRqPRKLhPy6oIAQAAzAxBCAAAdg1BCAAAdg1BCAAAdk2B+401NTWpqakZGRm1tbWzZ89udJvq6ur333//+PHjMTExjz32mLOzc9uPCwAA0HYKVIS7d+9+9NFHv//++3nz5jW1zX333ffTTz/dcsstW7duvf/++9t+UAAAAEUo9vhEWlra4MGDKyoqGv5RVlZWTEzMhQsXvL29i4uLO3TocPz48fDw8IZb4vEJK4LHJ9SFEfwqwslXkbU+PpGSkhIXF+ft7U1E3t7e3bt3T0lJMcNxAQAAmqVkqDblwoULfn5+dT/6+/tfuHChqY1lWb7rrrvqfrzpppseeuihpjaurKwUBIz3MYdL1exsGRVU0ZVqVq5n1RLVSpKrwLWO3NeZAl0oVCuHaklkajfUblRVVYmiqHYr7BROvopaWhE6Ojo2u7E5gtDZ2fn63s6ampobDJZhjE2aNKnux+jo6BtsrNfrMe7GREpq6bdC/ttFnnqZHy4mDaPOnixYS/7O5O5ATg5MI7ByPS+qFg+VUH4Fz6nQFFbxrl4swZcGtGdDAijUDaloQjqdDr/8asHJV1eLgtCYYskcQRgUFJSXl1f3Y15eXnBwcFMbM8buvvtuI/csCAIqQmXlV/D12fzbHPngFd6vHRsUILzUk/X0ZX4N/tfrdLIscyenv36FKvV0tJjvv8x/Os9fSJV9ndi4jmxiuJDgh0RUHn75VWRFJz8tLW3JkiVqt0IZQ4YMmTlzpvAnBfdswiBMTk7u0KFDVFTU8OHD77///gMHDvTu3Ts1NbWwsHDYsGGmOy60Qq1MG7Ll1afkjKt8fJjwXJw4NJC5tPC3w1VDffxZH3/2eDfiJKZd4d+ckyf/IomM/tlZ+GdnwR/foQHM6/Dhw5cvX545c6baDWmrtLS0rVu3muiDKBCEhYWFAwYMqKmpqaqqioyMDA4OTk5OJqK5c+fecccdzz77rLu7+xtvvHHbbbclJSXt2bNnyZIlWq227ccFRVytoQ+Pyx8el7r7sEe7CmNDBScl7n0wot5+rLefuKg3/X6JrzopR6/Tje8oPNNdiPFGgQhgPuHh4RMnTlS7FW3l7Ox8/PhxE+1cgSD09fXdvn37X3v8s+v2q6++qnu2YcaMGSNHjjx58uTbb78dFhbW9oNC212toTcPSatPyXeECTtv03T1MlU+9W/H+rcTl9aIK0/Kw3/QD2gvzO8ldPdBHAKARVAgCDUaTURERMP3Q0JCrv8xPDy80WcHwfyq9PTOUfmdo9LEcCHjTk2w1hyZ5ONEL8YLs2KElSflW37UjwgSXksUQsxyaACAG7CO+72goI3n5G4b9OlXeco4zYdJonlSsI6LhmbFCmfudghzp54b9QvT5Wos7gsAqkIQ2pG8Cj52m/7lNPnTweK6YWKkh2rVmJsDLUgQD96hybjK4zbqf71gQatDA4C9QRDai49Pygmb9P3aCel3aIYEWkSHZKgb2zBcfLuveH+yNHOPVKFXu0EAYJcQhLbvUhWN3aZfdUpOHqN5qYfgYGH/5mNC2eG7NDqZem7S77+M0hAAzM3CLoqgtJ0FPOEbfU9ftnusCceFtpGHA308SFzWRxi/Xf/mYRlhCGCHUlJS3nvvvTlz5pw8edLMh0YQ2ixO9FqGPC1Z+nyIuCBBtLRCsKFxHYXU2zXf58q3b5Ou1ardGgAwr1mzZh04cGDlypVZWVlmPrTFXx2hVcp0dMd2aft5OXW8pdwRNEaQlu0crYn2pD7f6o8VozIEsEHr16/Pycmp+3Hjxo3Z2dlElJKSsnbtWl9fX/M3CUFog86W8f7f6Tu60Y5RmgAXtVvTQhqBlvYVF/QShv2g/z4XWQhga1JTU5cuXWp4ffHixWnTpqm+rKw5Jt0Gc/r9Ep+wQ3qll/BIFyv+lnNPpBDlye7YLmWVCk/FWvEHAbA0+y/zlw6Y7+ndezsJ90X97b/w9OnTe/fubZhrc9WqVePHj/f39zdbexqFILQpm87Jj+6RvhiiGRFkNd2hTUnwY3vGibf9JOWU82V9RcHqPxCARejmxV6IM99iil296r8TERHRt2/f//3vf9OmTVu9evXnn39utsY0BUFoO/5zUl6QLm8dqenpayOhEaJlu8Zo7tiun/qrtOYm0RGVIUCbuTnQcLW/KM+YMWPx4sUBAQEuLi5JSUnqNoZwj9BmLDkkLzkkJ98m2kwKGng50o8jNZV6Gr9dX4Un7gFswm233VZQUDBnzpzp06czpv4lC0FoC+alSWvPyL+NVXPWNNNxFmn9MNHPiY3+SV+uU7s1ANBmoig+8sgjZ86cuffee+venDhxYmRkZG5u7kMPPRQZGXnkyBGztQddo1bvhf3ST/n81zEaG172ViPQmpvER/dII7fqfxypcXdQu0EA0DY1NTV33323j49P3TuffPKJXv9Xt4+Hh4fZGoMgtG5zUqVt5/nPt2l8ndRuiokJjFYOFGfukUZt1W8dqXFDFgJYp/Pnz3/33XcffPDBrl27rn/f3d1drSaha9SKvZwmbc3jO0bbfgoaMKIPk8QYb3bbT/pK3C8EsE7FxcV5eXnr1q2LiYlRuy1/QBBaqyWH5A3ZfLvdpKABI1qRJIa7szu262uwkCGAFYqNjX399ddvvvlmtRvyFwShVVp5Uv7PSXn7aNGG7ws2RWC0erDo4cim/CJJmHkGwCZcvXr1nXfemThx4rhx4954442KigpzHh1BaH3WZcsL0+Vto8QOrjY4RtQYIqMvh4plOj5jD6IQwBb8/vvv6enpEydOfOSRR7777rupU6ea8+gYLGNlfr3AH98rbRulscknJYznKNDG4Zqbf9DPT5NeTTDfNBkA0EaffPLJ4MGDO3XqZPhxzZo1AwYMGDNmzJgxYwzvBAcHJyYm6nQ6BwczDYpDRWhNjhXze3bqv7pZE+9j1ylo4OZAW27V/DeLf3xSVrstAGCsrKysZcuWGV7n5+c/8cQT9eYazcrKCgwMNFsKEipCK3KhksZsk97pJw61nmWVTM3fmX4cKQ7+Xh/ixkYG47QANK/23IlrG1YQmemugnbgWG3fW65/Z/r06XFxcUuXLnV3d1+1atWECRO8vb3r/vTy5ctPP/304sWLzdM8AwShdajU07ht+oejhcmRKOL/ppMHWz9Mc8cO/Y5Rmu4olAGa4xAU6XX3E2Y7nManfb13QkNDBw8e/NVXXz3wwAOffvrp119/XfdHxcXFt95667333vuPf/zDbC0kBKFV4ET3JUuxPmxuD6RgIwa0Z+/1F8dtl1LGadpb2/qLAGbGHBwdQ6LUbcOMGTNefvnlgIAALy+vfv36Gd4sKSkZOXLkkCFDXnvtNTO3BxdWK/BKmnSpiq8ciCEhTZoUIUyLYnftwMOFAFbg1ltvvXr16gsvvDB9+nTDOxUVFePGjYuPj3/rrbfM3x4EoaVbly1/kck3DNdgEaIbe6WXGOjKZu5BEgJYOkEQpk+fnpubW9cFun79+l27dq1bt87X19fHx8fHx6egoMBs7UHXqEU7XPTHwxJ2+OB8SzGiNTeJSd/p3z8uP94N3xoALNq1a9cmT57s6elp+HHy5Mnjxo27foO6PzIDBKHlKqqhO7ZL7/UX8bCEkbQa2jRCHPCdPt6HDQrASQOwRNnZ2evWrVu5cuXvv/9e96ajo6Ojo6NaTcIXZwslc5r6q/6OMDYpAv9GLRDuzj4borlnp1RQiTlnACxRRUVFRUXF999/Hx0drXZb/oCK0EItSJeq9PRGIgbItNgtQeyxbsLdP0u/3KZxwLcIAAsTGxsbGxurdiv+BtcJS7TtPF91iv/fzRoN/n1a5cUegrcTzUnFwBkAaB4utBYnv4Lfn6z/71AxAI/EtRYjWnuTZuM5/k0OZl8DgGYgCC2LXqYpv0hPxIiDMdajbbyd6Oubxem7pXNluFkIADeCILQs8w9Krhp6IQ7/Lgro48/mxIv3/CLpUBYCQNMwWMaC7Czga87wg+M1AqpBhcyKFXYWyP86IC3pg2FHYI80Gs13333Xu3dvtRvSViUlJd27dzfRzhGEluJKNU1Llj4dLLbDrUHlMKJPB2t6btKPCOLDg/D9AuzOpEmTunXrpnYrlBEeHm6iPSMILcVDv0mTI9kIXKyV5udMa24SpyVLGXdo/DBBD9gZR0fHhIQEtVth6XAvyiKsPCnnVfBFvdF9ZxLDOrApkezh3/A0BQA0AkGovtMlfN4B6cshIqbVNp1FvcWccr7qFIbNAEB9uPSqTC/Tvb9K83uJXbzQKWpCjgJ9MVScmyplleJpCgD4GwShyl7LkH2caAZWSzC9bl7spR7itGRJQhQCwHVw/VVT2hW+4oS0erCIYtA8nogRnERadhgdpADwFwShamokmpYsvdNP7OCKHDQTgdEng8W3jkhHi1EVAsAfEISqeeWg1MWLTY7EP4FZdXRjixPF+5MlPcpCACAiBKFa9l/mn52WPxyA5yVU8EC04O9MbxxCEgIAEYJQFTUS/TNZWt4fk8iogxH9Z5D43jF0kAIAEYJQFQvTpS5e7G4sPa+eEC17PVF8YBdGkAIAgtDsMq7yj0/JHyShU1RlD0YLHg70zlF0kALYOwShWelleug36c0+WHRXfYzo40HikkNSJh6xB7BvCEKzeueo7ONE90XhtFuEcHc2t4c4fTf6RwHsGq7I5nO2jL95WFo5EI/PW5AnY4QyHa05jQ5SAPuFIDSfR3dLL8SL4e7IQQsiMvp4kPhiqnSpSu2mAIBKEIRm8kWmfKWaZsXghFuceB82LUqYvQ+LNAHYKVyXzaGohp7fL60cKGpwvi3SK73EvYV8+3ncKwSwR7gwm8OcVGlCuJDoj05RC+WqofcHiDP3SNUoCwHsD4LQ5PYW8h/ysPq8pRsdwnr4ssUZSEIAu4MgNC29TDP3SG/1FTwc1G4KNOfdfsKHJ+QzJeggBbAvCELT+vdxuZ0LTcJsatYgSMvm9hAf34uiEMC+4AJtQgWVfHGG9D6WmLAeT3QTLlbRumw8VghgRxCEJvTcPvmRLkJnT4yRsRoagd4fID67Ty7Xqd0UADAXBKGpJF/gewr53B4oB63MoAB2UwBbhFEzAHYDQWgSepke3yu93U9w1ajdFGi5N/uKn5yST2HUDIB9QBCaxIcn5EBXujMMp9cqBbjQiz3EWb+jKASwC7hSK+9SFS1Kl5b3R6eoFXu8m5BbTt/lYNQMgO1DECrvpQPSvVFCVy+MkbFiDgIt7y/O3ifXoCwEsHUIQoUdvMK35Mkv90Q5aPWGB7E4H/Y2lrAHsHUIQiVxoqdSpIUJoqej2k0BJSzrK7x9RCqoxKgZAFuGIFTS11lylZ7+2Rln1UZEuLOHo4W5qSgKAWwZLtmKqdLTnFT5nX6igJuDNmRuD3FHAU+9jKIQwGYhCBWz7Ijctx0bFIAYtCluDrQoQZiVIiEJAWwVglAZBZX8vWPSkkScTxt0X5RQI9H/zqKDFMA24cKtjJcOyA9HC2HuKAdtkMDo7X7inFQZy/YC2CQEoQLSr/Kf8uUXMa2o7RocwBL82Lt4lALAFiEIFfBMijS/l+iOpXdt2pJE4a0j0qUqtdsBAEpTZk7oS5cuLVq06OTJk4mJiS+++KKbm1u9Dd5+++1Tp04ZXrdv337BggWKHNcSbM6VL1fTg9H4SmHjIj3YvZ2E+QelD5NQ+gPYFGUu32PHjq2oqHjppZeOHDny0EMPNdxgy5YtjLGEhISEhITY2FhFDmoJ9DI9v19e2lcUcXPQDszrKW44J5+4hgGkADZFgYpw7969mZmZe/bs0Wg0MTExQUFB+fn5wcHB9TYbPnz4hAkT2n44i/LxKTlESyODEYN2wduJ5sSLL+yXv7sFRSGA7VCgIkxLS+vTp49GoyEiPz+/yMjI9PT0hpt9+umn991339KlS0tLS9t+UEtQpqOF6dLSvrgm2pHHugnHivmvF1AUAtgOBSrCwsJCb2/vuh99fX0vXrxYb5uRI0d6e3s7ODh8+eWXq1evPnjwoKura6N7kySpZ8+edT+OHz9+9uzZTR26oqKCc9UuSYuOaoa2YxEO1WVlajVBTTqdTpbl2tpatRtibvNihWd+1/wyolbdfoDy8nJVj2/XcPJVVF1drdFoDKWXMZydnR0cmhnKqEAQurm5VVdX1/1YWVnp7u5eb5vnnnvO8GLKlCldu3b99ttvJ0+e3OjeBEFYtWpV3Y+BgYEN91aHMdZwYI55FFTy1Zn69Ds17lpnVRqgOkMQOjk5qd0Qc5sWQysy9T9c0t4TqfIIqRv81wBTw8lXi4ODQ4uC0BgK7CskJOTrr782vJZlOScnJzQ0tKmNHRwcwsLCLl++3NQGhjE1bW+Vqc0/KD/URQjR4u6g3WFEb/YRH9gl3RkuOGKwMID1U+D/8ZgxYzIzM1NTU4lo48aNbm5uffv2JaKUlJQtW7YQUVVV1fnz5w0bp6am/v777/369Wv7cVV04hr/NkeeE4+7g3bqpkDW1Ys+OoHn6wFsgQIVoaen5/Lly2+99dbo6OisrKy1a9eKokhEmzdvPnz48G233VZaWtqlS5fg4GCNRpOTk7Nw4cI+ffq0/bgqmpsqvxAnemHRQTv2Rh9x+A/6+zsLHphIAcDKMaUGm5SUlGRnZ0dFRWm1WsM7tbW1siw7OzsbXmdlZcmyHBkZaXinUXq93sXFRafTGXnQ8vJy898j3FvIp/winZyocbbvgtBu7xHW+ecuKURLCxLU+T0oKyvDbSq14OSrqKWDZYyhWBAqwiqCcPD3+gejhWlR9n53CEGYW857bdIfneAQ4KLC0XEtVhFOvopMEYT2fjVvqe9zeXENTe2E8wYU6sbu7ywsTMeaFADWDRf0FpA5zT0gvZ4oYEI1MHixh/i/s3JWqQV1qwBASyEIW+C/WbKHA40NxUmDP/g60VMx4stpGD4KYMVwTTdWrUyvpMmLE+17hAw0MCtW+OWCfKgIRSGAtUIQGuvjk3IXLxoUgF5R+Bs3B3oxXnwpFXcKAawVgtAoFXp6PUNe1BvlIDRielfh2DXaU4iiEMAqIQiN8u9j8qAA1tMX5SA0wlGg+b2EuSgKAawTgrB512rp7SPSggScK2jS1E7C5Wr6KR9FIYD1wcW9eW8dkcaGCp09UQ5Ck0RGCxKEfx2QkIQAVgdB2IzL1bTiuPxyL5woaMZd4YJM9M05PEoBYGVwfW/GG4ekKZ2Ejm4oB6EZjGhhgvhymiyjKgSwKgjCGymo5J+dll/EcktgnNEhzMOR/i8LRSGANUEQ3shrGfID0UKgq9rtAOuxqLf4arqsRxQCWA8EYZNyyvnXWfLzcSgHoQWGBrIQLa3NRBICWA0EYZMWpsszugl+TS6eCNC4hQniwnS5FlEIYCUQhI3LLOXf5sjPdEc5CC02oD3r4kmfnEISAlgHBGHjFhyUn4wRvRzVbgdYpwUJ4usZcjWmmgGwBgjCRpwq4T+dl5+KxcmBVkr0Zz182ccnURQCWAFc6xvx6kF5Vqzo4aB2O8CavZogvHFIrtKr3Q4AaA6CsL7j1/jOAvmJbjgz0CY9fVnfduwjFIUAFg+X+/pePSg/0110QzkIbTa/l7D0sFSJohDAsiEI/+ZoMd91QZ6JchCUEOfDktoLK06gKASwaLji/82rB+Vn4kStRu12gK14pZew7LBUgaIQwIIhCP9ypIjvKZRndsU5AcXEerNBAcKHx1EUAlguXPT/siBdfra76IpyEBT1Si/h7SMoCgEsF4LwD4Zy8FGUg6C0GG82OBB3CgEsF677f1iQLj+DchBMY15P4S3cKQSwVAhCIqKjxXz3RXkGykEwDcOdwo/MUhSeK+Nb8znWBgYwHi79REQL0+Vn4lAOggnN6yksM/0zhTnlvP93+uf2SS/sxzynAMZCENKxYp58AeUgmFZ3H5YUIJh69tH5B+VHu4q/jtF8elo+V4ayEMAouPrTogx5dnc8OwgmNztW+I8pg7Cklr45Jz8RI/g60X1RwurTGJ4DYBR7D8KT1/gvBXh2EMxhQHtWI9OxYlMVaptz5WFBgo8TEdF9UcK6s6gIAYxi7wGwKEN+KhYzi4KZjApmP+SZKp9+yOOjgpnhdZwPK9NRVimyEKB5dh2Ep0v4tnz5McwsCuYyPIj9XGCSHktOtLNAHh70RxAyomEd2M8FCEKA5tl1BryeIT8Rg3UHwXwGBggpl7hkgng6U8KdRdbRjdW9c1Mg230RQQjQPPsNwrNlfEue/GSM/Z4BMD9fJwp0Zaa4Tfj7Jd6/Pbv+nf7t2e+XEIQAzbPfGFicIc/sKng6qt0OsDN9/FnqZeXzKfUy7+v/tyDs4skuVfGiGsUPBWBr7DQIc8v5pnPyU7Gi2g0Bu5Pgx9KuKB+EB6/wBL+/BaHAqKcfO2iCYwHYGDsNwiWH5Ye7/DHQHMCcevqyjKsKh5PM6Ugxj/dl9d7v6csyihCEAM2wx8fICyr5V1nyyYkYJAMqiPNhR4u5zEmoH1utl1nK/Z1Zw2FfcT7sFwwcBWiOPVaEyw7L06IEf2e12wF2ydORvJ3YuXIl8+loMY/zaSRXu/uwoyZ7fh/AZthdEF6qos/OyM/F4e4gqCbGW+H5ZY4VUzevRt7v6sVOlZjkaQ0AW2J3QfjOUWlypBDoqnY7wI5182Inrim5w+PXeDfvRipCrYbau7BszL4NcEP2FYRFNfTxSfn5OPv61GBpunixE9eUDKeT13hXr8ZvOXbxImWPBWB77CsS/n1MHh8mhLopN0oBoOW6eLJTyoWTzOlMCY/2bCIIPdmpEqUOBWCb7GjUaJmOPjgu7R1nRx8ZLFO0FztdolgQ5ldwLyfW1MTxnT3ZQaWf1gCwMXZUEX54XB4RJHTyQDkIKvN3JomTUnO+nCmlzh5N/mlnTyVDF8Am2UsQVunp3aPS3B728nnBwkUpl09nSninJvpFiSjKk86gaxTghuwlGFadkge0F2IaG1kHYH6dPFimQosFZpbyG/RzBLmy4lpeoVfkUAC2yS6CsFampYdllINgOSI9KKtUmV1lllKnprtGBUbhbuwsVugFaJpdZMPaM3KMN9WbkhhARZHuTKnl47NKeeQN73xHerAsPEoI0DTbD0KJ0xuH5Jd6YCoZsCARHuysEuHEibLLeIT7jYIwwoPOKlR9Atgk2w/Cr8/KQa40MADlIFiQCHdSJAgvVZGrhtxvOIF8hLsyoQtgq2w8CDnR4gz5pZ4oB8GydHBl12qpqs1jWLLLePgNy0EiinDHLGsAN2LjQfhtjuws0i1BKAfBsgiMQrUKrEFhTBCGu1N2WRuPA2DLbDwIX8/AYFGwUGHudK7N+ZRdRuHuzR6I5ZRzlIQATbHlkNh+nlfq6faOtvwZwXqFuSlQEZ4r52HNVYRaDWk1dKmqjYcCsFm2HBKvZ0hz4gUF1wEHUFCYOzvX5lt358p4mBGTyCtyLABbZbNBuLeQ55bTPRE2+wHB2oW5UU55W3eSU04d3Yw5Fstpc/UJYKtsNicWH5Kejxc0Nvv5wOp1dG9r1ygnyqvgxiwr1lGJ0AWwVbYZFIeKePpVuj/KNj8d2IaObpTbtiAsrCJ3B3I1YmGxju6oCAGaZJtRsThDfjpWcMLTg2DBAlxYcQ3VSK3fQ04Z72jcKtOhWkIQAjTFBoPwTAn/5YI8vYsNfjSwJQKjIC3Lq2h9PuUa1y9KRKFuLBddowBNsMG0WHJYntlVbGrBbgDLEaqltuRTbjmFao3asqMba2M3LIANM+L2glXJr+DfnJNP340YBCsQ+kc+tfIRn9zy5qeVMfB2IolTqY488D8DoAFbqwjfOiL/s7Pg46R2OwCMEOpGuRWt/+t5FRRqxLMTfx4LRSFA42wqCK9U09oz8uzuGCQD1iFE26Zwyi039h4hEYW0rRsWwIbZVBAuPyZNDBcCXdVuB4BxQt1YXtuCMERrbBCGurVpYA6ADbOde4SlOvrohLzvdtv5RGDzQtwor7Vdo1V6KtdTOxejj6VtU+gC2DDbqQg/OiHfGizceKluAIsS0obHJwqqWJArM/7XvS2hC2DbbCQIqyV696g0J95GPg7YCXcH0jAqqmnN382vZCHGPTthgMEyAE2xkeT49LSc6C/EeqMcBCvT6lt356tYiNEjZYgoRIuKEKBxthCEepmWHpZfRDkIVihES3mtGsyZX9GyijBYy85XYHlegEbYQnh8dVYOc6N+7VAOgvUJbu1twpZWhM4iuTvQZSzPC9CAMkG4ffv2YcOG9ezZc/78+ZLUyCzC+fn5kydPjo+PnzJlSkFBgSIHNZA5vXFIfrEHnh0EqxTixvJbFYT5lWT8sxN1x8ITFAANKRCEOTk5d9111yOPPPL5559v3rx52bJlDbeZMGGCv7//unXrfH19J06c2PaD1tmcK7uINCII5SBYpVZ3jZ6vFFrUNUptG6QKYMMUCMLVq1ePHDly0qRJsbGxixYt+vDDD+ttkJaWdvz48aVLl3bu3HnZsmVHjhxJT09v+3ENFh+SX+xhCx28YJ9aXaWdr6TgllaErQ1dANumQIQcOXIkMTHR8DoxMTE3N7ekpKTeBnFxcU5OTkTk5OQUF5u2aG4AACAASURBVBd35MiRth+XiJILhdJaGt8RQQjWKtiV8ls+mLNMR3pO3i2cU7fV9yMBbJsC87BcunTJy8vL8Nrb25uICgsLPT09G93AsE1hYWFTe5MkKTw8vO7Hu+66a/78+U1tXLT5s2+upl48hv/b6uCcsxY80g3ketdjmrBu17/jxSm/wqmsvLxF5/FUKQt0FsvLW1bf+WuEtEtieXn9ATPVP35We3x/i3Zl5zjnZfjNV0lxLV0Y8nC//j2M3N7Z2VmjaSbpFAhCT0/Puv+QhhfXx55hg4qKv770lpWV1dvgeqIo/vzzz3U/+vn5ubk1OcH+sMl3t9Pc2eqWQ1vo9XpZlh0dHdVuiNUo2bTSsabC9e+/z25E7g66KtHN+MnSiKiohIdoa27wX6NRUb7802zJza1+IVldXOh56xSnqLgW7c2eVVRUaLUtvEMLSrhaTZN/1P3Yp33DX+O2UCAIw8PDz5w5Y3h9+vRpd3d3Pz+/hhsYqgfOeWZm5vU1X0MRERFGHtpD66Jp4bUAlMJ1OibLGicseWUs5tR41hl6LNu5tKDCyCvnwS2fXD5Y21Q3LBc9fDS+gS3eo70SHMs07u5qt8IeLU+VkqKl9lqFb4cpsLupU6euX7/+/PnzRPTee+9NmTJFEAQiWrly5e7du4no5ptv5pxv2LCBiNavX88YGzJkSNuPC2AbWvEERV4FD9a2+I5AkJZdqOQy7iSAdSqppY9PybO6yorvWYGKsH///jNmzIiJiXF1dY2IiNi0aZPh/Q0bNtTU1AwcOFCj0axdu/bee+999tlndTrdF1980WyPLYD9aMXkZ3nl1MuzxYHmKJCPE12s4h1ccX8LrM+HJ+TRIUJHt0YeVW8jZQLp1VdfnTNnTkVFxfWdotu2bat7PWzYsNzc3KKiIl9fX1HEw+8Afwlu+QJJ+RV8bKs6MkO0LK+cOmDNTrA2VXp676j0820mKaIU62l1cXGpd2uwHo1G065dO6QgQD2hLV8gKa+Cglxb00GEyWXASq0+LfdrJ3TzMklnBrooAVQWomV5FS1LtbyK1gyWIaJgrEEBVkgn07LD8v+GmaqOwqPoACpr6YQvxTWkYeSmaU1hh3XqwRr9N0vu5EF9/E11bxtBCKCyDlp2sYpLRsdTXgVv0boT18OqhGB1ZE5LTLyyAoIQQGWOAvk6sYuVxiZhbjm1dLrtOrhHCFbnmxzZw5GGdTDhUGcEIYD6QloyXiavgoe2tiIMxbzbYG3eOCTPMfG66whCAPWFaFmu0bfu8sp5S1cirBPgyq7W8Frln0gGMInt53mlnm438coKCEIA9YW6Ua7RFWFuBYW0dmJBkVGgKzuP3lGwEoszpDnxgqkngEAQAqgvtCWDOfPKeWhrK0IiCtVSLnpHwRqkXOLnyumeCJPnFIIQQH0hbi0Ip5xyCm3DVPMhbiwXFSFYg8WH5OfiBI3pYwpBCKC+UKPDSeJ0sYq3dG36vx0LFSFYgyNF/MBl/s/O5ggpBCGA+jq6GTtYpqCS+zkzhzb8x+3o3oKBOQBqeeOQPCtWcDbLpJwIQgD1+TlTpZ4q9M1vmVNGHdu2BGdHN5aDIATLllXKt5+XH+1qpoRCEAKojxGFGpdPuRW8Y2sfIjQIdaMcdI2CZXvzsDyjq+DuYKbDIQgBLEJHN8opa36zNo6UoT+7YVESgsU6X8HXZ8tPxppvqSIEIYBFCDOuIjxXxsPaVhFqNeSqoctVbdkHgAm9dUS+v7Pg62S+IyIIASxCR3d2zsggdG/r48VhbkYdC8D8rlTT2jPy7FizZhOCEMAihLnROSO6Rs+Vt3WwDBGFubNzZQhCsETvHZMmhAtBbXhAqBWwMC+ARQhzZ9llzcwBKnPKq2hr1ygRhbvTOYyXActTqqMVJ+R9t5s7mFARApgNI2qyDgs3omv0QiX3dCCXNl8lwtxYNipCsDwrjssjg4WINnf+txSCEMAitHehCh2V6260TXYZRXgocI0Id0cQgsWp0tO7RyVTr7jUKAQhgEVgf/SO3iifzpbxcCW+LEd40Fkj7kcCmNOqU3L/9kKMt7nLQUIQAliOCHd2trkgjHBX4EBhbiyvnEt1h+JEKlx8AP5SK9OyI/LcHupEEoIQwFJEeFBW6Y02yCqlSCW6Rp1E8nduwcJPAKb2+Rm5qxf19lPnGxmCEMBSdPJgmaU3CqesUh6p0DiCTh6UecPQBTAbidOSw/JLPcw3lUw9CEIAS9FsEGaW8k6eygRhpAfLwngZsAxfn5UDXGhQgGod9AhCAEtx4yqtpJaqJWrvosyxojzZmRIEIaiPEy3OULMcJAQhgOUIc2MXKnmN1Pifni7hUZ5Mqe/MnT3pDLpGwQJ8c0520dCtwWqO10IQAlgKjUBh7k32WJ4u4Z0V6hclos6e7NQ1VISgvtcyVBssWgdBCGBBoj3ZySby6VQJj/ZU7ECR7iynnOuamdMNwLS25vNamW7viCAEgD918aIT1xr/oxPXqItyFaGTSCFuLOuGY3MATG1RuvRSD0H1p1gRhAAWpKsXO9FERXi8mHdTdNKNrl7sOHpHQT2/XOBXqmlCuPoxpH4LAKBOjDc7XtxIOOlkyi7n0cpVhEQU40XHixXcH0DLLEqXXuwhiKrXgwhCAIvSzYudKuH6BrfuTpXwjm7MSdER5rE+7EhjoQtgBnsLeXYZ/SPSIjLIIhoBAAauGgrWstMNbt0dLuJxPgp/c+7uzQ4XIQhBHYsypDnxgsYyIsgyWgEAf4r3YRlX6+dTxlXew1fhIOzixXLLeaVe2b0CNO/AFX60iO7vbCkBZCntAACDXn7s4JX6QZh2hfdSOggdBIrxbiR0AUxtUbr8fLzgaDH5YzENAQAiIurtx1Iv/y2cZE5pV3iCCSbm7+3H9l9GEIJZHSriqZf5Q9EWlD4W1BQAIKI+7Vj61b896n7iGm/nwvyclT9Wv3Ys5RKCEMxqUbr8bJzgrObcovUhCAEsi4cDRXqwtOt6R3+7yJPam2SM+cAAtrsQQQjmc6yY774oT+9iWdFjWa0BACIaGsh2FvyVTzsL+NBAkwRhhDvTMKrQmWLfAI14LUN+urvoqlG7HX+HIASwOKNDhB/y/ugbrZVpZ4F8S7Cp/quOCmGXqlEUgjmcvMZ/LpBndrW43LG4BgHATYHsdAnPr+BEtP08j/VhAQotQ9jQnWHCxSpT7Rzgeq9lyE/FiG4OarejAQQhgMVxEGhShPDhcZmIPjgu3Rdlwv+nN3dgNRJlYvZtMLHTJfynfPnxGEsMHUtsE4BtYkRGx81zccKqU/KsFCmrlKZ2MuH/U5FRuDv9NwtBCKb1Wob8ZIzoYXnlICEIASxTqBv7ZoSGc9o6UjT1c8chWnamhKfjyXowmcxS/mOe/IRFloOEIASwWAPas+X9xXB3k0/OLzCaEC4sTMcqvWAqr2XIj8eIno5qt6MJCEIAoFHBbN8ljunWwBSySvn3ufKTlloOEoIQAIjIUaTn4oQFKArBBF7LkB/vJnpZajlICEIAMJjeRUBRCIrLKuWbc+WnYi06ayy6cQBgNi4aFIWgvEUWXw4SghAAiDgRoz+LQgwfBaVklvItufIsyy4HCUEIAHVcNPRCvPDqQRSFoIyF6RY9WLQOghAA/vJIF+HAFZ7WYGVggJY6XcJ/zJOfsuDBonWsoIkAYDbOIs2JF+YflNRuCFi9henyU7FWUA4SghAA6nk4Wjh0lbByPbTFyWt823mLfnbwetbRSgAwGyeR5vYQXklDUQit92q6PDtWdLfImUUbQhACQH0PRAsnS2gvFq+HVjlazH8tsNCFJhplNQ0FALNxFOhfPYSXURRCq7ySJj8XJ2otbBn6G0AQAkAjpkUJOeX06wUUhdAy6Vf5vst8huUtQ38D1tRWADAbjUCv9BLmoSiEFno5TZoTL7hYTzlICEIAaMrkSOFqNf2Uj6IQjJVyiR8pooejrSxZrKy5AGA2IqNXE4R5aRKSEIw0L036V0/BSVS7HS2EIASAJk0IF/QyfZuDSdegeb9e4OfKaFqU9cWK9bUYAMyGES3sLb6cJsuoCqE5/zogze8lOFhhqlhhkwHAjG4LYe4O9NVZFIVwI1vyeEktTY60ykyxykYDgDm91lt8JU3WIQqhCZxo3gFpYW9BYGo3pVUQhADQjCGBLNydPj2NJITGrTsrOwh0e0drDRRrbTeA9WGMuLXeanutt7gwXa7Sq90OsDx6mealya/1Fq2zGiRCEAKAMRL9WR9/9sEJFIVQ35ozcoiWhgdZbw4iCAHAOIt6C8sOSyW1arcDLEm1RAsOyq8nWtuTg3+HIAQAo3T1YqNDhGVHMOka/OWD43Jvf9bH34rLQUIQAoDx5vcSVhyXC6vUbgdYhpJaevOw9Fpvq88Rq/8AAGA2oW5sWmdhYTqKQiAiWnpYGhMidPWy7nKQEIQA0CIvxotfn5WzSq11+Cso5UIlfXRCnp9gCyGizGdIT09PSkry9fUdOXLkuXPnGm4wY8aM3n+66667FDkoAJifnzM9HSv+Kw3DR+3dqwelB6KFEK3Vl4OkSBDq9frx48fffffdOTk5cXFxU6dObbjN6dOn77777pUrV65cuXLBggVtPygAqGVWrPDbRZ52BUWh/TpVwjeek1+Mt+7BonUUWDxx27ZtRPTkk08yxl555RV/f//jx49369at3mYREREJCQltPxwAqMtVQy/3FF7YL+0YbVWrr4Jy5qbKz8WJ3k5qt0MhClSEJ06ciI+PZ4wRkVarjYqKOnnyZMPN5s2bFxMTc8899xw7dqztBwUAFT3QWThfSVuxZq9d+v0SP3CFPxFjC3cHDYz9QvfFF180fLNPnz6dO3cuKipyd3eve9PT0/PKlSv1tnziiSc6duzo4ODw2WefDR48+NixYwEBAY0eSK/XGzLV4NFHH33zzTebalVFRQW32jmrrJ1Op5NlubYWz1cbS6fTUXW1VFamyN7Ky8sV2Q8R6XW6qupqXQsb9kqs8FyKQ/9baqx0nuW2UPDkW6NnfnecGyPpKqt0ahy9urpao9FoNMaGl7Ozs4ODw423MXZfO3bsaPhmhw4dOnfu7OPjc/To0bo3S0pKfH196205fvx4w4ulS5f+/PPPP/zwwwMPPNB4gzQanc7Y08sYc3NzM3JjUJYhCJ2cbKVzxPR0Dg7Ozs6u131rbCN3hXZVo9G4uLg4t3Bv93ShDzP1my5q7+9sO5WB8ZQ6+VZn0zm5UpIfjnVS6wuQg4NDi4LQGMbua82aNU39UadOnVasWGF4XV1dffbs2aioqBvsSqvVoowAsAFL+4gTf5bujhBcca/QPuhkmpMq/3uAaGPdAAp8lRs1alR5efnatWtlWX7rrbc6d+4cFxdHRBs3bly+fDkRlZaW/vDDD5WVlTU1NR9//HFaWtrw4cPbflwAUFffdmxAe/bOUTxKYS/+c1IOd6dbrHl+7UYpEISOjo7r169fsmSJm5vb5s2bv/zyS8P7ubm5hlEzer3+1Vdf9ff39/LyWrVq1aZNmzp16tT24wKA6l5PFN49KmHSNXtQUkuL0qU3+9jIIxPXU6ZHY+DAgQ3Hgs6aNcvwwsfHZ9++fYocCAAsSoQ7mxYlvJImfTTQBq+PcL3Fh6TbQoU4H1srBwlTrAFAG73UQ/wmRz5WjPHbtuxcGV99Sl6YYJtfdxCEANAm3k40t4f47D7MxG3L5qTKT8aIga5qt8M0EIQA0FYzugrZZXi+3mbtLeR7C/kz3W02L2z2gwGA2TgItKyv+EyKpMcAUpvDiWbvk15PtOWHZBCEAKCAMaEsSEsfnUQS2povM2Ui+kcnWw4LW/5sAGBOb/cTF6VLRTVqtwOUU6GnuanyO/1EGxwqeh0EIQAoI9abTQgX5h/EqBnbseSQNDiQ9W9n2zmIIAQwH0Zk48NJXk0Qvz6LRylsxLkyvuKE/Eai7ceE7X9CADAbXyea11OclYKi0BY8t19+KkYMtok16G8MQQgASnq0i3Cxkjadw6gZ67azgB+8wp+Ns4uMsIsPCQBmoxHovQHiM/vkKr3aTYHW0sv01O/SW/0EZ9ucSaY+BCEAKGxoIEv0Z28eRlForT44LndwpfEd7SUg7OVzAoA5LesrvH9cyi7DqBnrU1hFr2VIy/vbRzFIRAhCADCFEC2b3V18OgVFofV5fr/0z85CFy/bHyNTB0EIACbxTHfhxDW+JQ9FoTXZfZH/UsDn9bSjcpAQhABgIo4C/XuA+NTvEkbNWAu9TI/tld7qJ7g5qN0U80IQAoCp3BLEevmyNw7hsULr8N4xOcCFJobbXS7Y3QcGAHN6p5/w4Qn5dAk6SC1dfgV/45D0/gD76hQ1QBACgAkFadlLPcTH9qIotHSzUuSZ3YQoTzsaI1MHQQgApvV4N+Fq9R+r+YBl2pLHjxTxOfH2WA4SghAATE0j0MqB4nP7sUKTharQ0+N7pQ+TRDuZR6YhBCEAmFyiP5sQLjy/Hx2kluiVNGlQezasgz12ihogCAHAHF7rLW7L58kXMGrGshy8wr/MlN/uZ6/FIBEhCAHAPNwd6P0BwiO7pWqUhRZDL9PDu6UlfUQ/Z7WboioEIQCYybiOQrwPW5iOJLQUbx+V/Zzovih7DwJ7//wAYE7vDRBXnZIzrqKDVH1nSvjSw9LKgXbdKWqAIAQA8wlwoTf7iA/skvR4mEJVMqeHfpP+1UMMc7ffMTJ1EIQAYFbTooT2LrQEqxWqasUJWeL0RAwigAhBCADm959B4vKj0rFidJCqI7uMv3pQWj1YFFANEhGCEADML0TLXust3p+MDlIVcKIHd0nPx4vRdjmbWqMQhACggoe6CL7O6CBVwQfH5RqZno7Fxf8vOBcA5sIYcXQG/oERrRokvndMwghSczpTwhcclD4dLIqoBq+DIAQAdQRr2dI+4n3JUg0eLDQLidO0ZGleT7EzOkX/DkEIAKq5L0qI8mDz0pCE5rDkkKx1oMcxUrQBnBEAUNPKgeKXmZiD1OTSrvDlx6RPB6NPtBEIQgBQk58zrRokTkuWrtWq3RTbVamnqb9Ky/uJwVrkYCMQhACgslEhbGwom7EHHaSm8sw+qbcfuycSF/zG4bwAgPre7CMeLeKfncHTFMr7Jkfels8/SMKcok1CEAKA+lw09N+h4nP7pNMluFmopPwK/uhu6cuhooeD2k2xYAhCALAI3X3Yqwni5F/wNIViJE7/+EV6Klbs1w63Bm8EQQgAlmJGVyHCnT27D0mojPkHJSeRXojDdb4ZOEEAYEE+HiT+kMc3ZONmYVttP88/Pc0/H6LBzNrNQhACgAXxcqT/DRNn7pUyS3GzsPXyK/i0ZP2XQ8T2Lmo3xRogCAHAsiT4sfm9xAk7pCq92k2xTjqZJu2UnogRbwpEMWgUBCEAWJwZXYXuPuxRPFnYKs/uk/yc2Zx4XN6NhTMFAJZo5UDx0FX+wXHcLGyZLzPlH/P5ZzdhKrUWQBACgCVy1dDGEeLCdGn3RdwsNFb6Vf50irRxuOjlqHZTrAqCEAAsVIQ7WztEM2mnlFeBLGzepSq6c4e0IkmM9UY12DIIQgCwXLcEsWe6C7dvkyoxcOaGamWa8LP+3k7srnBc1VsMpwwALNrs7kIPX3ZfsiSjLGzao7slf2f2agImFG0NBCEAWLqPBoqXq/hLBzCItHFLDsmHivjaIRgg00oIQgC7xzkxi76EOgq0cYRmwzm++hQGkda3Llv+8IS8+RZRq1G7KVYLZw4ArICvE225RRz8vT5Yy24NtujYNqfdF/nje6VtozQdXHFOWg8VIQBYhyhPtnGE5r5kfdoV3C0kIjp+jU/4Wf/FEE28D1KwTRCEAGA1+rdj/xkojt2mP2P3yxbmVfBRW6VlfcURQUjBtkLXKABYk9s7Cleq6dat0m9jxCCtnWbA5Wq69UdpVqwwtROKGQXgJAKYCSOy9ypGIQ9GCzO7CSN+lC5Xq90UNZTU0sit+jvD2NOxuIArA+cRAKzPs92FCeFsxA/6ohq1m2Je5Toa/ZN+YHu2qDceGVQMghAArNKCBPGWYHbLj/piu8nCCj2N/knf3Zu92x8pqCQEIQBYqzf7iEMC2Ygf7aIuLNPRqK36aE/2YRIenFcYghAArNiyvuLNHdjNW/SXqtRuiildq6Vbf9R382L/GSQKiEGlIQgBwLq92UccH8Zu2qK31UUqLlXRzVv0fduxFQNRC5oEghAArN78XuIjXYRBm6UT12wtC7PL+KDv9eM6snf6IQVNBc8RAoAteDpW8HOmm7fo1w/XJLW3kchIv8rHbpPmxgszu6FoMSGcXACwEfd2EtYO0dy5Q//1WVuYm3tLHh+5Vf/+AKSgyaEiBADbMSKI/TxaM26bdPIavdzLioeVvH1EfueovOVWTW8/6/0QVgNfNADApsR6s323a3YWyBN2SGU6tVvTclV6uu9X6auz8r7bRaSgeSAIAcDW+DvTjtGaDq7U51v9sWJrGj6TWcoHbNYLjHaNwcpK5oMgBAAb5CDQvweI/+oh3PyD3lqW8/0qSx64Wf9oV2HNTaIzpo4xI9wjBACb9Y9OQqI/m/yL9GM+/yhJ9HNWu0FNKKmlJ3+XDlzm20Zp4rC4oNmhIgQAW9bZk+0dq4l0p/iN+k3nLLE03JrP4zbqPRzowHikoDpQEQKAjXMSaUkfcXwYf3CX9Hkmf6+/EGwZCxkWVtEz+6TfC/mng8WbO1hEk+wTKkIAsAv927H0OzRxPtRzk37xIblaUrMxOpnePSp336ALcqUjd2mQgupCEAKAvXASaX4vcd/tmgOXeZd1+s/OyJLZh5TKnL4+K8ds0G87L+8ao1nSR3RFx5zazPovoNPpHBwczHlEAIB6ItzZhuHinkI+N1VanCG/2EOYEik4mL4okDitz5Zfy5BdNbQiSRyGKtBiKPCPn5WVNW3atPj4+NjY2Ka2SUlJ6dy5s6+vb3R09P79+9t+UACAtkhqz5LHaFYMFL/IlCO+1i8+JJtuIaeiGnr7iBz1P/0Hx+UliWLKOA1S0KIoEIR6vb5Xr14zZ87Mzs5udANJkiZPnvz888+XlpY+++yzkydPlmVLHLsFAPZmaCDbPkrz/a1iVinvsl535w5p4zm5Sq/Mzmsk2pwrT/5Fivxal3GVf32zuGuMZlQIItDiKBCE0dHRTz31VO/evZvaIDk5ubq6+oEHHiCiBx98sKKi4rfffmv7cQEAFBHvw1YNEs/d4zA2lK04IQf+V3fnDumjE/LJVi3qlFnKV52S7/5ZCvyv7q0j8k0B7Owkh7VDxER/RKCFMsc9wszMzC5dugiCQESCIHTp0uXMmTM33XRTU9sXFxfXvXZ1dXVycjJDIwFMjjFeWy1XliuyM15VISu0Ph2XFKqArJyHA/2zs/DPzkJRDW3Nl7fl8yWH5bJa3suPdfdhUR4sxI11cCVPR9LUMn0NEVFxLS/TUUEF5VXwMyX8aDE/eJU7CmxIILstlH2Q5OBvqY/ww/WMCsK8vLyvvvqq4fsPPPCAr69vs3/92rVrWq227kd3d/fro64eSZIiIiLqfpw6derixYub2riioqLZo4OJ6HQ6WZZ1Oiuc1VglkrN75ZY1JVvWKLI3znkZU6rCYDXMQV+uTELbAEeice1pXHsiosJqdqiInSgV9l+kTVXsYhUrraXiWgfGap1FchHJy5HaOfNgV95RSw9EyMt7yR1c/9yRnnBSFVddXa3RaDQaY6s4Z2fnZjc2al86na6oqKjh+5Jk1JM4fn5+paWldT8WFxe3a9euqY1FUbxBTDbk5uZm/MagIEMQol43ntsdD9MdDyu1t7KyMnd3d6X2Bk1xc6NIv/pv4uSrSPMnJfdpzEYRERE3KMua1bVr16NHjxqendDpdMeOHevatWur9wYAAKAgBQbL1NTU7NixY//+/ZIk7dixY8+ePYb3Z8+evWHDBiLq379/SEjIwoULr1y5snDhwrCwsD59+rT9uAAAAG2nQHVZVla2ZMkSIho0aNCSJUsCAgKSkpKISKfT6fV/3ITfsGHD448/3qtXr9jY2PXr17f9oAAAAIpgnFvQqpV6vd7FxcX48RcrVqy4/fbbO3ToYNJWQaP27NlTUlIyevRotRtij0pLS1euXPncc8+p3RA7tXjx4lmzZrm4uKjdEHu0adOmkJCQGzyw1wrWPdfo119/ffLkSbVbYad27969Y8cOtVthp/Lz8z/55BO1W2G/Pvroo8LCQrVbYad++umnlJQUZfdp3UEIAADQRghCAACwawhCAACwa5Y1WEaSJCcnp44dOxq5/YULF7y9vZ2dMYuRCkpKSiRJ8vHxUbsh9kin0xUWFgYHB6vdEDuVl5fXoUMHURTVbog9unr1qqOjo/ETGkyZMmXhwoU33sayVoQURfHs2bN1D100q6amBjObqEWSJM65svM7gPHwy68inHwV6fV6xpjx30ICAwOb3cayKkIAAAAzwz1CAACwawhCAACwawhCAACwawhCAACwa1Y55K+8vDw9Pf3UqVPR0dGDBg1qdJv8/PxPP/20tLT0rrvu6tevn5lbaNuuXLmyevXqwsLCUaNGjRgxouEGn332WU1NjeF1ZGTksGHDzNtAW/Ptt98mJycHBQU98sgjjY4aP3bs2JdffikIwtSpU7t06WL+Ftqwc+fOrVmzprKyctKkSQkJCfX+tKqq6vPPP6/7MSEhoeE20Do1NTWHDx8+evSor6/vuHHjGt2mpKTk448/LigoGDp06NixY1t9LKusCJ988smZM2cuWbLkiy++aHSDK1euJCYmXr58OTAwcNSoUZgSU0E1NTVJSUnHjh0LDw+fNm3a9VeBOk8//fTBgwfPnj179uzZS5cumb+RtuSdd955+umnIyMj9+7dO2zYMFmW621w7NixAQMGGJbh7tev36lTp1Rpp00qKChITEwsLy/39/cfNmxY3RpzdcrKymbMKSZ/dAAABg5JREFUmHH2T40uYA6ts3z58ilTpixfvvzNN99sdANJkoYOHZqSkhIZGfnUU0/9+9//bv3BuBUyPME2e/bsRx55pNENlixZMmrUKMPr5cuXDxs2zHyNs3VffPFFfHy8LMuc8/Xr13ft2tXw+nre3t6ZmZlqtM7W1NbWBgYG7ty5k3Ou0+nCwsK2bt1ab5sHH3xw1qxZhtePPfbYjBkzzN1K2/Xyyy9PnDjR8Pr1118fN25cvQ0KCws1Go3Z22UXDNf5VatWJSUlNbrB999/HxERodfrOefbt28PDg7W6XStO5ZVVoSC0Eyzd+3adcsttxhejxgx4rfffuN4XFIhu3btGjFiBGOMiEaMGHHixIlGa77PP//83Xff3bt3r9kbaFNOnz5dVFQ0ePBgItJoNEOHDk1OTq63TXJycl0H9YgRIxpuAK1mzLnlnL/33nvvv//+kSNHzNs6G9fsdT45Ofnmm282PFk/ZMiQS5cuZWVltfJYrftrFu7ChQv+/v6G1+3atautrb1y5Yq6TbIZ159bDw8PZ2fnCxcu1NsmKSmpqqoqOzv79ttvx5p5bXHx4kVfX9+6STTat29fUFBQb5t6v+0N/zmg1eqd25KSksrKyus3EARhxIgRhYWFhw4dSkpK+uijj9Ropp26/l9Ho9H4+Pi0+pffQgfLTJ8+ffXq1fXe7NmzZ2pqqjF/XaPR1M3TZnjh6OiobAttWHp6emJiYsP3d+/e3a9fv+vPraH7ouG53bx5s+HFo48+Ghsb+8QTT4SGhpq0zbbq+rNNRDqdruHMXvV+2/GrrqB655YxVm9OQT8/vx9//NHweuzYsf/4xz8efvhhzEFqHg4ODpIk1f2o0+la/ctvoRXhypUr9Q0YmYJEFBQUVPfF+fz581qt1tPT02SNtTU9e/ZsePL1er1h8O3157awsFCn091gKr+uXbt6e3tnZ2ebqek2p0OHDkVFRVVVVYYfz58/3/Bs1/tt79Chg1mbaNPqnVt/f/8bXGqTkpLKy8sxOsxsgoKCzp8/b3hdUVFx7dq1Vv/yW2gQtkJ1dfXOnTtra2uJaOzYsRs3bjR8WVi3bl1bhtVCPWPHjv3+++8NHUTr168fNGiQt7c3ER06dMgQeNXV1XV3ZPfu3VtaWhodHa1ig61ap06doqKiNm3aRETXrl3bvn27YRx5UVHRrl27DNuMGzdu3bp1htf4bVfW2LFjN2zYYBipe/25TUlJMaxQX/cdhYg2b97s5+cXEBCgSlPtR3JycnFxMRGNHTt227ZtJSUlRLRp06Zu3bqFhYW1cqetG2OjrjVr1iQkJLRv397Pzy8hIWHFihWcc8NV+MKFC5zzysrKPn36DBw48J577mnXrt3Ro0fVbrLtkCRp9OjRPXv2vPfee319fZOTkw3vjxw5ct68eZzz7777LioqatKkSePGjXNzc3v33XdVba/V++677/z8/O6///5u3brde++9hje3bt3q7u5ueF1QUBAWFjZmzJjRo0d36tSpsLBQvcbamtLS0vj4+CFDhkycODEgIODMmTOG9yMiIj7//HPO+bJly+Lj46dOnTp8+HBPT89Nmzap2l6bsmvXroSEhI4dO7q5uSUkJDz77LOG911cXHbs2GF4fc8993Tr1m3atGl+fn5btmxp9bGscvWJixcv1lXERBQYGNihQ4fa2tqMjIyePXs6ODgQUW1t7c6dO0tLS4cNG+br66teY22QLMu//vrr5cuXBw0aVNcXcfr0aa1WGxQUJEnSoUOHTp8+7erqmpCQEBQUpG5rbUBOTs7evXuDg4MHDhxoGK9bWlp65syZume3y8vLd+zYwRgbPny4VqtVtbG2pqam5ueff66srBw+fLiXl5fhzcOHDwcFBfn6+tbU1Bw4cCA3N9fb2zsxMRGXGgWVlJRkZmbW/ejl5RUZGUlEaWlpUVFRHh4eRMQ5/+23386fP5+UlNSWgQhWGYQAAABKsZ17hAAAAK2AIAQAALuGIAQAALuGIAQAALuGIAQAALuGIAQAALuGIAQAALuGIAQAALuGIAQAALuGIAQAALuGIAQAALv2/ySPAw5T6fV1AAAAAElFTkSuQmCC", "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/html": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 121, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# - Solves the Schrodinger equation in one dimension, for a potential given \n", "# by the function 'potential(x,xmax)', for x in the range (-xmax,xmax). \n", "# - The integration of the wave functipn Psi(x) starts from boundary conditions \n", "# applied to Psi(-xmax) and Psi(-xmax+h), where h is the integration step. \n", "# - The boundary condition at x=xmax is taken to be Psi(xmax)=0. \n", "# - Bisection is used to find the energy (preceded by a course search \n", "# for two energies bracketing a solution). \n", " \n", "function potential(x::Float64,xmax::Float64)::Float64\n", " v::Float64=0\n", " if abs(x) <= xmax/10\n", " v=-1000.\n", " end\n", " return v\n", " end \n", "\n", " function normalize(n::Int64,h::Float64,psi)\n", " norm::Float64=psi[1]^2+psi[n]^2\n", " for i=2:2:n-3\n", " norm=norm+4*psi[i]^2+2*psi[i+1]^2\n", " end \n", " norm=norm+4*psi[n-1]^2\n", " norm=1/(norm*h/3)^0.5\n", " psi.=psi.*norm\n", " return nothing\n", " end\n", "\n", "# This function solves the SE for given energy (ee)\n", "# The wave function is stored in the vector psi (function argument)\n", "\n", " function numerov(nx::Int64,xmax::Float64,ee::Float64,psi)\n", "\n", " h::Float64=xmax/nx\n", " h2::Float64=h^2\n", " h12::Float64=h2/12\n", "\n", " psi[1]=0\n", " psi[2]=1.\n", "\n", " fn::Float64 = 2*(potential(-xmax,xmax)-ee) \n", " q0::Float64 = psi[1]*(1-h12*fn)\n", " fn = 2*(potential(-xmax+h,xmax)-ee) \n", " q1::Float64 = psi[2]*(1-h12*fn)\n", "\n", " for n=3:2*nx+1\n", " q2::Float64 = h2*fn*psi[n-1]+2*q1-q0\n", " fn = 2*(potential((n-1)*h-xmax,xmax)-ee) \n", " psi[n] = q2/(1-h12*fn)\n", " q0=q1\n", " q1=q2\n", " end \n", " normalize(2*nx+1,h,psi)\n", "\n", " end \n", "\n", "# This function starts with a course monotonic search for an energy\n", "# bracket, starting at e1 and using incremental steps de. In a subsequent\n", "# bisection procedure, the energy window is halved until the difference is\n", "# less than eps.\n", "\n", "function matchboundary(nx,xmax,e1,de,eps,psi)\n", "\n", " numerov(nx,xmax,e1,psi) \n", " b1=psi[2*nx+1] \n", " e2=e1\n", " b2=b1\n", " for i=1:200\n", " e2=e2+de\n", " numerov(nx,xmax,e2,psi)\n", " b2=psi[2*nx+1]\n", " if b1*b2 < 0 \n", " break\n", " end\n", " e1=e2\n", " b1=b2\n", " end\n", " println(e1,\" \",e2)\n", "\n", " while abs(e1-e2) > eps\n", " e3=(e1+e2)/2\n", " numerov(nx,xmax,e3,psi)\n", " b3=psi[2*nx+1]\n", " if b3*b1 <= 0 \n", " e2=e3 \n", " b2=b3 \n", " else\n", " e1=e3\n", " b1=b3\n", " end\n", " println(b1,\" \",b2)\n", " end\n", "\n", " return (e1+e2)/2\n", "\n", " end\n", "\n", "# Main program\n", "\n", " nx=10000 # number of x intervals on each sode of x=0\n", " xmax=1.0 # hard walls at x=-xmax and x=+xmax\n", " e0=-20. # initial energy in search for bracket\n", " de=1. # step size in search for bracket\n", " eps=10^(-6) # energy tolerance \n", "\n", " dx=xmax/nx \n", "\n", " psi = zeros(Float64,2*nx+1)\n", " x = zeros(Float64,2*nx+1)\n", " v = zeros(Float64,2*nx+1)\n", "\n", " ee=matchboundary(nx,xmax,e0,de,eps,psi)\n", " println(ee)\n", "\n", " for i=1:2*nx+1\n", " x[i]=(i-1)*dx-xmax\n", " v[i]=potential(x[i],xmax)\n", " end \n", "\n", "# The potential and the wave function are plotted in the same graph\n", "# - they are normalized to occupy the same scale for better visibility\n", "\n", " v .= v./maximum(abs.(v))\n", " psi .= psi./maximum(abs.(psi))\n", "\n", " using Plots\n", "\n", " plot(x,psi)\n", " plot!(x,v)\n", "\n" ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia 1.9.3", "language": "julia", "name": "julia-1.9" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.9.3" } }, "nbformat": 4, "nbformat_minor": 4 }