{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Determining Structure from Microwave Spectroscopy\"\"\"\n",
    "\n",
    "__authors__ = \"D. Brandon Magers\"\n",
    "__credits__ = [\"W. D. Allen\"]\n",
    "__email__   = [\"bmagers@belhaven.edu\"]\n",
    "\n",
    "__copyright__ = \"(c) 2008-2020, The Psi4Education Developers\"\n",
    "__license__   = \"BSD-3-Clause\"\n",
    "__date__      = \"2020-07-13\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Determining Structure from Microwave Spectroscopy"
   ]
  },
  {
   "attachments": {
    "ethenone.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAALQAAAEZCAIAAAHGXah3AAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAAIdUAACHVAQSctJ0AAGR3SURBVHhe7d13kF3lmSZw/riSTkvtNB7P2jM1O1s1NVO7OzVhXTuzu7N27RpsTAaRk0UywUQTDBgw2SYnkYQQCJBEzlHknHPOGSSRc5AJvb8+z+2PS7daqGV1i6vtt26dOvfcc77wvM/7vO93blqs66ts5qRJzb3+bW6tPF9VH3Z2vtvZ+XRVXVZVzaNzsn5b+eMbb8zu7Oz67nc/7eyc1dl5e1Ud3Gg0X+tj/bbySFV93Nn5eWen4bw8evRNVTV+Plp5ecyYmWPGvNnZOXP0aC1eVFV7zkcrb3R23lpV91XVHVV1SVUdNWLE1f1D028rr3V2auiEqppUVfs3Gu91ds5PK693dprOyfWVj1XV+/Pdisfknlbmcyw3/O//feOPfvT7//yf7V/43//77T/+8aR/+qe81Nf6beX3ta2zzjr2d9555/333/9Xv/pVXuprX93K1ltvnVY222yz5mt97CtaOeCAA/bbb7/ddttNE5tvvnnztT7Wbyu/+93vDj300KOOOurggw/ed999t99++7XXXrv5Wh/rt5W9995bKxMnTiyt/OIXv2i+1sf6bWWXXXYxFzMyL/tbbbXVGmus0Xytj/XbChR++9vfQgS022yzzXrrrffzn/+8+Vof67cVTtl44415d5NNNjGXZZdd1nCar/WxflsxnX322WfppZdeaqmlfvKTnxx44IHz04qx7LXXXqussor9X//61wCaz1a4abXVVrO/3XbbzWcrZsTBq6++un1ImxES56W+1m8rf/jDH7QS726xxRbz2UrfsfB3Xuprc8OFj/bcc8899thjhx12mM9WjMWVhx12mK1BbbvttvPTiiFo4uijjz7kkEO08pvf/GbddddtvtbH5taKUWjCdvfddzeWKNYcrd9WdI4vQlFzO+20E31ZeeWVN9poo+bLX7Y5tPLOO+8cccQRaKZ/uNpuuumma665pjiA96677vrmm282T+2xL7Vy//33f/7553aE8i9/+Uudw4I4jR07VhP8TbreeOONGTNm3Hzzzbkk1t2KFz744INnn33W/ltvvXXbbbcZSP1qbyN9991338yZMw3nww8/NOp7773X8cXee++99O+1Sy+99LTTTjvyyCPrS/q1Cy644Mknn7Tz8ccfP/744xMmTFjs7rvv1r/hPPLII9OmTct5X2mPPfaY4ej41ltvxc/uGd1www1mcc455+y444456StNEzfeeONVV11lIJ420Z00adILL7yQ/Xm3zIvNwdOt9t7dd7/4VXXd3Jp4p66iXu3svLeqLlh++ebRPtZvEyqXrm99q+ub31RtPFVV581HqaBmUQ52/du/2VGIaeL3/RRic27i048+eqKq1CyvdHY+U1U3V9VJVbXrQJt4d8yYe6oKCjdW1WlVtUuj8duBNvFRXQvq/KiRI7dvNNTLWmm+/GXrtwl17f09EEL0hY6OATfBnQ/0NPFBZ+eAR/HJRx/d+eMfn/vDH+apyu2Kf/u3jf/6r/O0l825idmzZx900EEELU/JIu342c9+lqe9rN8mXEbHDj/8cFUBNaGsP/3pT5svf9nm1oRwto24U7YBN+H68ePHU1Y7RJQsD6wJidaVrqcLZqEJ01lmmWWaL3/Z5twE7TV+F0t6koOJSAsDg3PDDTeEpfnr3/VmoTBZfPHFmy9/2ebcxAYbbGAirmT6V0xIC8stt1zz5S/bnJtwGV8oqczf+M3IdAbWhIloQnLLUyPCi/kZRWki1FJw5mkvm3MTp59+uiZUdXnKL+qZFVdcMU972dyaULnn6fw0IUVqgl/zVBPK1vkchdIDNdSGPDLgJpSXqh8rEkOwkFDID6yJM844wyhUBsceeyx3oKmBDCxGhJnOjUI9ZzhWJKqfgTWh9CJTtoINFghiIYLjzZe/bL2beOKJJyiViYhuAQbRLbfcUsgsv/zyhqPdK664onlqj32piffff99WNaRbFddaa62lczGWtYx60PWvv/66suDcc8/NJWwxFZMiSuGhVrFvFCqW/sC/8sorXa+VV155RcX08ssvO3mxt99+22uqt4cffviWW27h/5zdnykm77zzzk8//dS+tq655prF1Dda0Z7KkPNy3lfarFmz1GnqtjPPPHMxLLr99ts1ltppHm369Ok33XTTeeedx19NOI8//vjszIfNmRfzYo984xsvVNWMMWNe6uh4sqpUEpdV1dtPPNF8eYA2P+NQfL3b2flZz2N2fafu5dGjH62qG6rq1P7rsrnYgMehdlHEfZKa8Hvf61p66c/rcczs7Hysqm6vqvOr6sgRI2Z/8EHzgnmzgY1DDXR37QsdK4BUZUqit+objspDFd41dYV3cKPRX5HYnw14HIrSO6vK1J+tKoWY7u2rzpSJl9ZOOazR2HMIxvFGZ6fK8oGODqy8sqququ9Vnl1Vp1TV+BEjdI+81/ZfKPZn84MHdxgNokyuqt80Gurk7RqNaSNGqHY9gHT1EIzDCDKO++p6vflCV5dgaR3HoPvF6injwIkpXx4H5g4RHopcSU8ClqikrtZ7FzvvvLN8nCJLYuuv5uzPFvA4mERqHP1V4P3Z/IzDAsbWOJRliiNblY1aa6uttip4DPo4IKGkA4naKmZYqfNUOArp3DMddL9kEEorSyFl2jHHHHPwwQdzh0pJvoeKQay77roqleY182YDG4fp8oLZG4pCzwjgYRDcxB3BY5NNNll77bUVTc1r5s0GNo71118/aymoAIApiZRN+LHjjjsahMoTGCuttFJ/xWN/NrBxWAoBA/JKRjw1JlWn7hWeRmCUa665pmXNXFZp/dn8jINT4AF8U1e6il4rTJVs7sRzkMAZ3HFkgWocwhUkuX8f8zQv0Q/j6G+x158NGA9xqzO0sERcddVVmy/0jANng4fFQfOFebP5HAe/6Lh1HIaVcaDtoOOhtA8/RIqO8/5BbPPNNy/U2WabbYRM84V5swGPo/hFpOSdmZjlZqHO0I0DHoIWHlOmTOEFkVzwyDj6W8j1ZwMbh8Vq8UvYwOyAR5hQFAOKX1ZYYYXmNfNmAx4HCTcUon7kkUfKL7JMyS+GAqTNNtuMogyuX84//3z5Ra/SyhFHHHHcccddddVVSTFlHCAZN27c4OqYuUolvBBKxikYY3CSTlLML3/5S+o+lzda52jzOg7dqDNgoDOTztu7zCrbvoM4YZRSDPISeMeN2MG5vH/cal8xjpkzZ+b2wocffgh/4Kv/hKhiR69cwCT6jTbaSKIBA3dIMUYGISeDatasWc8+++y999572mmnpc052hfjePTRR5955hnbG2+88Y033vjss89y/KOPPpoxY8ZDDz103XXXIeASSyzx8MMPO3733XcrMnRsywu6Z7lEMQYYYDz++OMvvfTSa6+99vrrr7/66qsvvPDC/ffff8sttzj49NNP33fffQiXS7rHoZs//vGPdnJLxwyc/d5777388svmMX36dHwsNxwHZJdffvnUqVNvuukm03v++ee1rP289PHHH5vtY489dsUVV0B6MXsG6ygAjOOdd9558cUXdW+kCCGJ5LI/0RSUJmNi0P3kk0/0pSNPH3jgAfMkAd143HXXXU899RQv2sL85ptvPvfcc2nDghpEMR7XPo8wjuagiy66aMKECTpq8kOv02o74YQTaJTcEU8tcOOL448/fvLkyeCBt45y/Cvi5Stt1qRJMyZNeqt+n/VPsfkcx6fvvPNsVc0YPXrGmDHPVdUjVXXT/N6Ris3POF7eeutX6yX1h52d79X3g17q6LDOvrmqps3vUOZnHLkF8lln5+ff+Man9Wgcgcrd9S2hIwe40o8NeBxP1fcXut91/c53ur71LaP5uIbkxY6OB+v7Y5Or6o89IjHvNrBxfEJdOjqM4/P6I4Fdf/d3Xf/wD/B4wziqKuM4ZeA3YdjAxgF8/b1T37fsHsr3vtf1F3+BInk7+d6quqIex+6DPQ5kfLqqzL7ctDSIbjA6OoTMbfVHFCdW1Y5DMI77q+qFjo7XOrs/T2gEAuel0aOfrMG4uqrOqD9wue9gj+Oxqnq0o+OBqsJWPiIh4IHEXVV1Xf0xgeNrp7Tev5tHG/A4yAYy6hgw99WxemuNhEGIlP0bjUtHjkSR5gXzbAMeB2bgxAU1Ja+q32O4sHYHWuzXaEwZORJaQzeOF+qOATCpqiZU1aGNxu8ajcc6OgiacUwd7HGUm7Xl5rXQ2L7R2LbRiMja8lc+DzsgG/A4Ch5YUvTqkw8/FDu5v2scQ+QXj143r40DEgTeOIaOH31vXhsHOSn3u4duHAS01ziMIHjIMoPOjwvqjz97XPk//scJ//RP5aM0EuyNP/rRXT/+se05P/xhPis9IBvYOMrN/NxEX3LJJXO81/3u1vvu82jzMw7W600F47B6s+K1gBu6cXyN8LDCDh4Wj1tttdV+++23EMYBjz322MNQ9qoNNo7/+te/zjis/Qd9HBbywQMPbHMLxH65yb/55psPHR68cPDBBx9++OGHHXYYDIwGDFDZaaedOCt3upsXzLMNeBziQt9GYClqfWwnYRJItt566w033HCgN8fYwMZhxroExvjx448++uiJEycaTes40PaXv/xl633VebQB8yPjOKI2q/WpU6cWv0BL7GQcg3t/3TiiV4aSjz/xkfDZe++9hQ+Gbrnllvyy6qqrlntD82gDG4ee4A8SfSOsrX0jK6Gb96NWXHHFwR0H2EFCPPTN8qYYJAxi++2332KLLYCx5pprLr300s5sXjNvNrBx7LDDDiat11133TU3LXfZZRfk3W677Xhk4403FrEB4+yzz25eM282sHGYPRoKzm233RYAzAjECPly3CBWXnnln/70p+uvv/7gjgM5kHSDDTbYZJNNSBZDCDCst956a6211tixY2U+gDlh0PEwDnq6xhprmD1K2hoBGJZbbrkllljCS9w06OMIHsIELThimdp+/vOfL7744nKvl4SPcQyRX0waQ7mjNThJSMaBxaJmSPFoHYcwzhCHdBzw6DUOWlLwGCKeUlJ4bLTRRr3GQeOHdBxficeg85SG6myOPFUg5qVANZevMs7RBjYO3WfSOutvHBIQng7FOOIXktofHoM+jvhljuPYc889Cx6D7pcyDjtzxMNLQ8GPXniQc7lXxj/xxBNb8Rgiv+jMjqHou3Ud5fiQjgMeQqZ1HWVAujcUZdiQ8lSdbOWSOjlDIfkyrdJw3LhxQ6Ef0U0jsHJhWdIZB6kFidrMIJQjgzuOrKN0nHXUsccem3UUWmQcXEPfVl999cEdh56MI4tb66jzzjsv38OIX7A167l80ad5zbzZgMeROsMiKusoY0IXzhI1Vjd4iqTWc4M7DusG8zb7BEsCJ4tKQ1TEb7bZZtb7+T5s85p5swGMA+wkS380Ayq6z1LKIMSqNYRXkVTZrGjlnQG9HT5P47jxxhvNmxdwou86ytOyjsqHHBXuCWOaO4/vi3/1OF566aW7775bfB5Tm7ig5ToGADe1rqMMwgrbEgZUxq1yNtA777xz8uTJzbb6t7mN4/7773/33XfzjeQzzzwTG8xSx73WUfazjiIb+SAIqPjLUOSdl19++bHHHrv++utt0+wcbc7j+PTTT1999dVPPvkkTz/++OOZM2fCn1yatKWKyKQQZR0lQKyjrCgNQrxgBhoR+9dee+2NN96YNWvWc889Z1blUyd9rTmO119//Zlnnnnqqaecfemll3744YeBgb3zzjtaAa8O+L4EwlJLLbVsbVb3COE4fjhu9HTMiO+9915tGsSbb76pfTPR/u233/7iiy8+//zzTzzxxE033fTee++lte5xePmDDz4waV4AYz52wxz09MEHH7zqqqtOOOEE8vDWW2/lpa80JD399NNvueWWRx55RK/QBYwthGbPnp0Pwzz55JPXXXfdSSed5PzFnGRQ+TQQbhv7008/baRm8Pjjj2sIM/AO6er2B2YyAHTRXJszZszQZjrKB29eeOGFu+66y3AdWQx0Rldf1eUk43BE92L14osvlkEQM6/On0EUluYN5lZ3wz6ReM455wCvexxeRkxn2HnllVfA4DLBhvZXX311LvsTTQAj6dtvv23OCMQ1HIQPt95662mnnYbUi3G/Q1BiBvHss88aIySOPPJI3mk2syBMl7DhGqPRI+7fd999V1xxxcSJE7s/H3TZZZc99NBDumdCHMmx8pRTTqHWzQYWkPH4PffcAwNc0ZH9a665Ztq0aWQGBN3xMn36dIQQUThxySWXTJkyRQ4THbl+ARrK8zhgyJpODUK6ThrqHgdaCB5+OvnkkxGT/gzGIGIKKCI7adIkfpcvSy78Qk8JC2K3fmVy8IyU9fqhhy/GMZT2x1mz7qs/2uTxSL19sP7Iwi3127EXdXR89NprzVOH1oYUjtcnTny6qmbV7/K+2dn5dr3zVmfna/UHe14cPfrZGpp7a1ymV9XpVXXzPP/SxwKxIYLj7TPPfLH+LM979Zvvn9QfPrP1+Lj+lNG7NS5AeaGqnqy/nnhn/UGG86vqhKq6a/z4ZkODbEMBx6N/9Vcvjx5ttk0gvvnN/HhJ97b+GNwf65fyuauZnZ048kTNkduq6vL6kzbHV9UBzh98G3Q4Hlt8cbwwTxT4NJ98Cxbf+U7Xn/9511/+5eff+lY+oAiOJkE6Op6q1eTW+hM3F9Qfnjii0fjNwD8KNlAbXDgeHj36oap6qf72sgknOrqx8Pj2t7v+/d+7Nt+8a8klwfFBZ+c7NTtmjB79XFU9Xn8kLB8GEy+n1N+n3aPR+N3AP1QyIBtEOD796CPCaVbPVxWfE05zphQQmV1vxYiHI+/3COrMMWOcXLTjxlpQBcuJ9Recf9to7D3IBBlcOKQJ6VMqNUmiwPlAoZriwgMdoOBggBBT+Zo3QgWLK3qoccyoUb+vP57W3y8LLSgbXDjM7fXOzjuq6uHa56LgpdGjyarJzxozZsaYMfZf6uhwHI/IJyAoqBiRUy6tqnPrrzpPGDVq/0Zj50bjpFGj5uMTvgOyQWcHOMTCdfWPgAocuDhIGjzMH16OOH5PzYib64+vyiaAOLX+7CTJ+EMtogCaPF8feB6QDToc+Zw9gXimo8M8ZYqb6ioLBcr2hloyyYQkcmb9EXyZ9aiRIw+sP0Q6ccQIUrIowMH5pCFw2GLK3R0das2z6t9f8Dinnj+xxAXZdGL9G7IHNRqUgkwcN3Kky+mOpAtHItL27CCTEkrKCnAoOoXGVVVlzofUj0MbDfv7NRr71L+MsUMtE7BLYZZLwIE+Jy0C7DCf5Fdz4+rAwdVznNsnH34YQrkkcChYCK1gAYdgGegvRQzUBp0dmVtWJZiC+RECobHbnODIJdjhksRXgsUli0KwmE+Yb2ueEm2Yb259Xd0KRy5xOXYsUsGSuWGHeYb5/aWJBItLUKPAYQmjHgHHHAm1YG3Q4YiUJrOAw2qVdpib+mqOcJT4yiUuL9oxx0sWrA0iHLNnz965/pTsAQcccFD9Iczf1z8K8pv6S+zrrLNO+ZB5sVzitANrc6H937V82bzvJQvWBh0O82mFI++lBo6+P1zjkp122imXBA5o5pLNN998jpcsWBtcOFrnNu/s+EPLBye+EsEFa0MKh3lmbv0xP3AkvgJHuWRRYEevYDE37CAEXL3eeuvtsssuv63NDuDy8eKtttrK/GHhksIOLwWORUc7ytxs4bLffvvZBqN8MGGPPfbwKnRwIW/KO1K0o7/4WrA26HAU5kMBBHH7wQcfXLaOMOcEFxDstttuQNl+++233nprYbXRRhsVOJZYYolm64NjQ8EOMzRVczb/+msYRxx55JHj6x9KPvzww/MZHq8Wpjh/1113de12220ndjbZZJP16x+oanvtoAioHiwOOeQQkwcEO6r+mTSIsMMOO8xLECkE2bP+LRcEAQeC4AWCwGLllVdu42C5//77TWnvvfdOjJiwmYcREDn++ONPP/30iy+++OSTT86HAFvZAQ7sSLyAI5/sWW211ZZZZhmvNjsYBBtEOMihlAEOkzTViIWZM4iwE088ccqUKfn0Xys1EEqwYBZ2bLHFFptuumnYETgG+umzAdkgwiHsEyz5eBREoqlFPpkjLFmmFxaFGhtvvHF+9ydfXWlXOEhgPqJlhuZptuZcLBAQWscxKEDk83c77rijC5FLWoGFCiXCsWz9W38D/c7bgGwQ4cBwvlVZmZ55mi1Thpn5XnvtZUsFPHWQUqQkk1AFiEu23HJLWJCMfBRPmMBCWgHHQL8aMSAbRDg22GADs+VeoHA18sNFFPB/MU8dVHR6NSWp6otYAEJyRYrVV199pZVWWnrppRdffHFYLLfccu0KB3bwvIigEfZN0lS53Zyhw0gDCwQgS2g4M1/QBUS+ZZaPaAo9+OJIe7MjcDBSKlny+YYbbsj5zMztO4057tW1116bZEIBIySRRIdG8AhqTl6k4CAZNMKsTHiVVVZZtTY7UgZbccUVzXappZYKCsBK0qEygYOOtDEcpl2CJXDIuPRS1jBVQJhzXxMUck3qFHC4HIiBA4PaWDvAwbEpSQsc2EEyRUp/XysAltMAkfLE5RJQgaPtg6XAYYfbsSNwhB3NU1vMq4GjXBU4aPCiAEdf7Zg7HK3siHZESsGhwUUHDnNLsEQ7KGh/cCBR2BE4WoNl0dGOAoe6EzvmDgcq5RJohh3ggG8bw2H9bjHWCod5OgIO6aO/b4pZwoqpclXYQXEUbPBVjwz0e24DskGHowSLiYHD8kQ9rvpUdPYHB04lWJjLwYFTqlicUpu0MRwmH9qzwAGgwDEXdrSGmJ0IMDiwo+3haJ1YCZa5s6M1WMKOBEt7s+OMM84w+cIOO/OiHRJqX3aAAzvko/bWjl7skFkcCTssW02SyRq2Vvp52wkL7PSFY1GQUpMvUlrYYeugeTIH7f/+97+HFFJY45CJ3/zmN9KtqElRGykFh0Tb3sHSCw77UkbuIWebHQeZV4Fi8q5CBxwBQW6jhx2LgnZwr0nyP0RMPu8ttNohhxyS1VphCkRcCILcPd50003xYs0112xvdqTu4Goz5HxYHH744UfVPw3Hjj766LztBJFD6zfiQhAnE4tEDYJQ0BBk3XXXVcguvfTSbQwHD5tb4sWETTuMOOaYY4499thp06adf/75p5xySn7wARyoRFmwAxyJl8BBesEhGVnCydDNDgbBBhEOPA8cHC4QSrDkbSegTJo0CWSgcdyrsCjUQCu5Rg6KfEixeasFHHNMzwvKBhEOqWHHHXcULOWNOP4PKBGLWEjhBLwoWMi+ocbmm2+uSFlvvfXEyMorr5w/NWp2MAg2iHBInOoIBIEI/ie/MpMv5ilGJMtKq1FQvJBoCxa5q7zqqqsuv/zy+cb/BwP8V8d5twUPh0maHi5ELCCijggoPG/OZo4FttEI5qUAgU0CxCVq0yhoeLHaaqutsMIKP/vZz2CBOM6B1DrrrHP33Xc3e11AtsDguO+++/gf86OXyRr59raJcXV528m0i5lb3nnyahghxGRW2imtRj7zdmRur+sC1okpqcdVzqdBV155ZXMcf5r9SXDMnj37qaeeeuONNz766KP8rgEIJNFTTz315JNPJhBhCufTws0228zQuZ068j8zeeapg3nnCQoYMW7cuAinzIoU0QsnaIol4gou+QGM++s/h7nwwgslrObg5svmB478X05+cSF/vVMMItJnEkTKbRMGROoo3kZ++9AxbVtmJeJIBAL/nUYmrE0wIu/XM3hRVlQSU3hBm2Ghlzdre+WVVwzp5ZdffvLJJx988MHbb799+vTpEyZMmI//XZkzHCb5em35cZQcfOyxx1588cV3331XN+UnMGKeOu7Ml156CV8uu+wyPk9GUCaonUxJQelM+6aK/+TASzEsyDtvy9X/GBSNiGkK48aOHQssQaRZASK+wP3CCy/MnDkTFvrNaF999VVH4MJbjzzyyF133XXddddxj4Pvvffe008/fccddyDRtddee8899+QvpXvZF3CY0uOPPz5jxoxAENQ1/cwzzzz88MO33norLvRCAWrvvPMOz4Sx+rj55psvvfRSwUI+8t6qEGie3cfMM3NuNUFUftull5EnyWWttdYSLNdff73ZOhMj4GLYhgERwzZ4U/CUb5577jmTcpC3DNWOl/KHS+LLpMTXJZdc0uygwBEJ+PDDD3neJNknn3xi/ppwHCk0arZgdjLJeOutt/QHKU4g7zfccMMFF1wABfUVMlP+BfVTKHM0UiVDHX/88Yraa665hs/FCH8gBVx4ES6GzUw+XM7/OZmUa/ObKUA0eIjwH8lPkuqGA8Cu7CUEacLFNAIKenrooYdcfMUVV4QvnGMo5513Xj7clWLBaJrXD4kZMHkiE2eeeeZVV10lEDCI54yWCznM/FsZbUZ8KUzghSMmDgVEc7k1lP3F0ImZs/O0HhRzMVyFHKSETNilP1haZZx77rnSm9xB4VgrjgvLaITok9eN8MYbb+RCAcLHrTNiIgC1scOsH330Uczi4LPOOktKkgEX42c087KLncowwhYQMCJOrcokvZ9zzjkTJ05Mkmv28DUz5L333nvFgsGbAi5kXnbMES9QOFgkzMlHlpGItpjouu222xCMHEIUwUATidaipl2GgXgBxZDCQoM66KDZ/9fPLr/88kSNSXGnSZmO6BBBXOs4puMFgSN5U6ZMUStlJdGtHRdffLEX7rzzTk3QJHxhDzzwgKcOShbCkmhNnToVo1RWsABqOv7a2tlnny2hqkEwxVzMyw6Dwk033WRGSEH7Md1igvChhquamYWrgSQQBB4WyPZ2LrroIkfOOOOMyZMnH3XUUYDImoI+5aqvudF+I4cLX5oLr9shE1A44YQTzMiSQl1nRhiUS76oO2IENqsvJihkeCwilq7Jb6+1o5lUbjKkXLZWkgTNqO/vT/WGY4jtvbvvnjVp0tO77PLSpEmvXXXVZwubdwsHjue33fbh+hcaHq9/sOGhL/9+ybUD/BPxBWhDDcdL22zzdFU939GR/2icOWbMy2PGvFj/5dpjNSh31L9ZcX5V3bTBBs1rhtCGFA5EKH99l592ebv+pqwHXF6o/40uNMkvIJ1UVZ/2/DTj0NjQwfFETQdAvNfzvesP668PexpQun/4p6NDBD1YVbfX3yM+t6qOHzWqtaYcbBsiOPBCdKADCPILSJ/Wv6SWXwwDCkRe7+yEV6Lm3vqXXi6rqtOq6uiRI5utDL4NBRxP/a//9dLo0fmtrNk9f2ba/ej5bbl83R5Y3f8c2dFBX4krEUGQ8+qfhTpgkH8fq9igw2EJSTvNEwVCis/ze3J/9mdd3/ymfUc+rl/Nb2W9PHq08+Wdu6rq+vr/G6fWv6H21vPPN1scTBt0OEgjapgqCoADHbp/TM7j29/uftQAgQNx3u2Jl2fq7HtPnXcvrn8K59hRo/Zx4eDboMOB+TRSHjFn0fE5CL773SYiCPKd7wBIBIHDOdhBYlrhIB/gOG7UqL0bjY97fh158Gxw4XjwP/5HuVOkSB/YQSa6f38RIoD4q7/q+j//p+vnP+/63vfAETVNsIDjkVpNb6rhOKv+oaj9Go2D/vZvm+0Omg0uHC/UQgAOgYAd4OgmSH5z8Pvf7xo3rmvLLbv+8i+LlM7q+f1FUnp3CxyTqGmjceiIEc12B80GEY77/uEfIgSzxozheXDkpwZFR5JL17/8S9d/+S+OoIZIiXA831OeRkovqnOtYNm/0dh3xIiTB/PTDGwQ4bj7Bz/ADmUV/r9RTxgLIqggAI0dD3HkJSegBtHFJsKRHyu9qq7Wp1TVMSNH7tNo7NFoHPDP/9xsfXBsEOG4/gc/eLn+QWPVt3hBENMu9aiHnWAhTEKN5+p13f11pCjD8mt74DhyxIg9Gw0EObB94biuhgPtOfzFjo5XepYqEEmdXoDwEgZB7ckeEVWkW7ZcUv/mnJXLoY3Gbo3GwSNGtDkco0dbj6mplN6cLxzkDqB4iA4PTx2MZCRM7qtVIyv9s+sabEL9R947NhpHtDUcN9RwPN7RwdvmSVP53xFlSB5KjG5SdHQACy+cI0zAJ6FQDSJ6ev27z4c3GiJlZ3DItW0NB2mUYqlAEJFBgUIgygMjAgQG4YV1ipOvrLEQJlYr1m+osVOjQVOxo42l9KYeOJ7q6Li1rjKt3EmDPEovPex4mt8rFSD0QozAQmEuTE6pf7j1kPo3OmnHBTVN2jhYZBbLU3opiVibmqoJmzamIAKJtQMjR+RUpLimTiWmre6STSaOGnVYo7FX/RuVWGNpC46D21pKCxwKLVMFyg09P99qhrZQcEQSIRbyiDnTi/ACI/atf8f28VGjnAym9pbSG+tgUW7AQloVNWYrEKQMEWH+HnYurafqJaSQRybXegGL39fZRNIRVmhFOxaRYAFHynDZRBRYoZIG87dVZUFBGY4RgEAKc6adClDyCUHZx5oYOyKl7c0O1Rd2pAwNHOamrJI+ZQ07ILCvsjimVgrrNAGyS51ZIehal2CH7NsNxyIjpWEH5kurXH1S/YPG5Zd9D6xDAwQYsV2jce6oUblEzZpgiXa0d7B0a0dHh3rc3FKYFziii1IGdbAFwfb1r8J7PNeDYAjlEoSSetoejm52jB5tSkU7rPSji4HjgH/8x+apPXb9979fLsmiRrBAMFIqWNq4DEuiJYdhPiFoZf4c4bju+99XtpdklEta4Wh77cjcPOASV5c0MWd29MSXS8KOkmjbHo5klghBNxz13aDm3EaM2L8fOAqhwIEdVjpBcHxbJ9qyZsncuuHo7KSLgeNwc5tTsEBQQgmhChzNS9qaHd0r2pr56o7C/C+EoD921FIaQmGWSwLHIpFZ+tQdVrGBA/P70w48cr6r7BR2tD0cuTkIhVKVqrhLoh0/cmRfOG5oYQdEuuOrJxl1E6qtE22CpUhp2FGC5ch+Ei3taIWjSGl3bm7ru2FTf/Sjq/7n/7z1xz++s37Yufbf//3if/3X0/7bfxv/D/+w89/+7W//7/9tntpjU370oyvrS+6oH7f1XHLqv/zL+P/6X7svWXzx5qmDY4MIx/bbb7/HHnvs1/Nt+v3r3/7Zddddf/3rX29U/xH4+uuv3zy1x7bbbrs9v/xjUvlezDb1r1W4ZINB/oTU4MLRd275EOMvf/nLsWPH9oUjCP6h5XvqQRAcEHRJ28NhSpkbXAoc/bGj1yUHHHBA2IFQ2NHecOywww5c3XduXwlHIVQQDByLGjvCfOwwN8EyF+1oJZSnRW5WWmmlRYcd5rb33nuDowjBHOHoq76LCDvAMXch6AtHLmmFw1OXBMH2ZkfSRCs7EizJmuAYN26cI55utdVWW265JWpsscUW/bGj7aW0PziiHTvWP4DsKct3rZ0faEqidYn4Kuxo77qjV7CUuf3mN7/Zfffd7TMAgaz7e/j1D3YABUcwxQkugUvYUUqV9mZHLyHIN+h+3/MPHDFP961/GxwikkgQQQeB49VCqEUhs7Syw06AoKnMkew47mAQQQpcCCKiZpNNNgk7EiyLTqLNnLOfr1LFDqoNKKEJ+iRqKIsAQZCNN944P2GQUmVRCBZAMHM280MOOeSwFuv1s1A4AsHd6v+bSMhsttlm8jEU2r7uSGaJz802WLT+usfhLX/OwgAHDtEBDgShuEiRH/5Zp7a2zyz5LRduN9vwAhb172R1/1CWfeZ4EAmP5hgvv/jFL1ZbbbX21g7upYIcHskIL4455hhAsICCI4Ug3RHVA0cRVHDk53/ys1BtDIesyQQL1ShwBIWActxxx9kpcEQ+WuHAjuQXcKy11lorrriieGm2Pjg2uOwwMTOMcJg2QwdAQOGEE0648MILzz77bE8LHLCjHa1wYEfUFDvA8dP2/W8nwW9iZlh0NNoRjkyZMuXqq6++9NJL8+89oUaBo2iHCrUECziWXHLJjwfzy4KDBcdrr73GvaSUdtBIEy6ISCjhyNSpUydNmuSgV0vd4RIhJtFSYpll8803j5QGjp/97Gfvvvtus49BsEGEg3ulTDM01cgHRMQFAwFL+BRqpOgokSJPE478ehg41lhjDZlliSWWaFc4zMfEJNoQxJzNnEUpYqlHGWoo0lODiRTUSKTkr2vWXXfd1VdfXWZZfPHF2xKOV1991Xw4mbfNEyLRVPOP2YeCg+FFFizBwlUpSRMphCM/OJhfmGpLOBjfihdCYJ5my/9mHi6wABG9IJ/OiYJSDbRKTiGiStISKfkZsbaUUobnZmV6fG6q5swQAQS2DBAhhVehFsmIglq2hRr5zcFVVllFpOT3F5utD44NIhxmhSBcbZJma86kweSpCVxs7RcsklkjGRQ0WGy44YbJKaiRXyFrYzhMD9tJAETsmzCamLlt2cGIBIgTEiOwUHe5EBbjxo1TjFKN8hOMbQwHUqCAuZkhpvC8I2bOAgGDQgEibMoPMZKMYLHqqqum+gLE8ssv38ZwmCeBEAtExDzRRL5gcFG/20IhQHgpQMAuMQILK3qrWFgkoRCO/n54fwHaoMOhuDBhkzRbyQJTTJ4RyxggKIVXBYgSQ1qlF/mnQFgstdRSgFCMerXt4ZBEwKEYRQqIMP43cwYC+4jDHA8QSJGfbM1P+eIFyVhuueWQKHAow5qtD44NLhwqi8Ch7sICs91ggw3EAmlgILAPBdkUI0RHSDF27FjaKZXghRRLdAsci0KwMHCAxhFZk/9jIGDKCig4jhGAWGGFFZAi2kl3tCABERpwOGcRCZbAocpwUOI0c4YItvwPBSYolllmmVJfuFDxCo7CDsuWNoZDWg0crMDhINVQVqGAyQsKW3pJLFN0MjGVq6xrwg5wkJhFLVjCDnCIESzwNPMvRk2UJFEcBg6AlmBpbzgQIXCUYFGkhx3goA7N81rssssuExpiJHC4qrBj0YEjEws7FGDgoKOio3lei/WCowTLoiClfdlRtCPB0jyvxQJHCZZWdkQ72rvu6A8O7JgLHGFHrtKCdLtISWmBw8QEi1J97uxYNLUDHHNMtIFjXqSUduxT/ydBgaONgwUR+pZhFvVzySzTp0/vy44ES9vDMce6o0jpXNiBFLkqdUcrHIuCdrACh4NzySzYscsuuzgzV6FJyjALYnAsIkV6CRZlmGChHf1lljkGS2uibWM4zIEQmlIvOLbq+Qv3a6+9tnlqi2FHKxwJFmXYZpttZuHbxm9ZB44yscBhSbLllluCY5lllukPjhIsRTuo8q9+9StwWOY1zxscG1w4UD3sMLfAIYLmAw7sCBxtnFlKsJiYrRAgpfPOjoQYOFyFHYLFmqW92QEOHs7EAgftCBz9aQe8emmHqwo72lg7ZMcES+bG54FDZll//fXnwo7AAcQES+DAjrXXXru9pRQcYUfgKIl2LsHS6/aPFiKlySyLgpSamG2BI8Gy9NJLfyUcrcHS9nCYQ6+6w8QGlFlETcksbQ8H7WjNLK3soB1fyQ4gggO/gJhgoR1tnGj7sqPAMW7cOJnlzDPPnDp16tZbb51PMGyzzTaXXHLJHKW0aEd7Z5Ze7DCxoh3Wcojgqa0I8tT5oCG0JVhclWAhQ4uClPZix5577mn+WACXvWtzBGVM2DotuLgKWZwGjl7a0d5lmDm0ZhZzgwIzQ9B4yux46rR8KgoogIBIPlHqhAJH269ZWoOFmblpmyFFyBHmVU8zbYhgigV+EMERi/q96k/RZQmHHW2sHZHSMm1wBIiDe77kVHYCCkScb/LUlJRs2/PhQfsJlkWEHSZsqsHC/qH1X3t3f9ep5dtOCZyk1YSMy6kMFDbaaCOghB2LgnaIArPt/vB1/U/nRxxxRD6WbqcVkRDEycSiEISCbFz/67tSpb3hCDvIAbebrWkHC3bUUUcdc8wxdlq/3lMIEjhkGXBEQdQpa6yxxuqrr97e7MB8iJin2SZGjqy/9hUDCoAcxw4WOCIfBBUcAAWHeNlggw0IxworrNDe7OBk3o5qYAc4MAIigGD5thNEwo7IbeTDheSjsGPDDTdUoY8dO7aN4TAZokgOCjsSLCHI8ccff84555x99tkTJkwowVKSSy92BI6VVlrpZz/72aD+heMgwqF2AIfpFXYEDhxBjcmTJ19ZG46EGuBADYVGtEOtUdhBRwPHkksu2a7fWDAfEzPDsAMFgkhAgcipp546ZcoUT0ONRArpTcEeOGQWWRYctEOw/PznPwdfs4NBsMGVUpyXOM2T8805oAQXJnyYg+gTaljCQDCR4nIruiRamUWWTbC0KxyChZPNMARJyAQR29QaYHIcXlGNXtSIcOTbThJtMku7wqGmlB14G0HMNhyJgYDZcZDBKwuWqIarUMPliRRZFhyrrrrqMsss85Of/KQt4Xj11VetOBIvqT7MOZNvNTC1YoEXwQI1tuz5IlxqMMKx1FJLgeOEE05o9jEINihwmL8YYdbppof/jkgZpg0XEMQChJcAUXiRMCmqUaghUixnwWHxIss0e1rQtoDhuOuuu8yK2+kCmQSH2oEu4ogJo0CQAk32QwqvuspprVighsXbL+pvO6HG0vXXOLzkHFuscUmz1wVnCwwO1ZGBpspIQpVKLUlIAE31Es+bNkMWlh3bkMKCzWkFiyxk11tvPfkVNZZbbjk5ZfHFFw98zneyYATWo48+2hzBgrAFA8fjjz+u0Ex9AQgoqLUcOeOMM4SDBGGeQQQFbGlEzNPcznACc2Z4keQqKPI7FcoN1IgAJb7gqEHsc/KJJ57YHMefbH8qHB9++OEjjzzy+uuvf/DBBwGCKcOnTp06bdo0OVXgcDVd5HaBw0zeTAJBGGFWXuVtuiCV4AUs1llnHUtYWEgoeIEmtfh+8V1TNAHonXfeeeutt55//vkK3OaY/gT7k+C44447Zs6c+c477/zxj3/87LPPZs+eDQvKjxTqTnnU0I0bHXjbVAOK8LEFQXbCCK9CzWmWJ0IAL4IFyRAmqEFxWnOTfS2Lymeeeebhhx++++67r7nmGg54//33m4ObL5tPOF577bUnnnjirbfe+vjjjwHRPNrVddFFF6GukAGECaC0WDBtCYL48Txum3nM/FlCw6vOQYpoJ71YccUV8SLZBHaaIhytCVsXs2bN4o8XX3zxqaeeevDBB2+//fZLLrnk5JNPbo5m4DZgOD799FP8fOWVV/jBfvNojxFUUmrQmEwXBAJGpHbgcFvOZ/mqE7PjKcmEQr7tFFIsv/zy+RYcLGy1A1aiq1mkgAi79tpr33zzTYgYzMsvv/zCCy+QsHvvvffGG2+0Vr755pubYxqIDQwOMqFXpBAXn9fWfKGryz6A3njjDWEcnUMKzgeBGearO3xuwuqIGAgKCoxqrrzyykiRb37RC1gwLWCQsCI3xCI0IahQoFl6ZPbhgiZixyBF8VVXXXXKKadQtOb45s3mDAeBBPyMGTMUl1lQOyJEHdRBX1IIGRgZ0LPPPivz0UuxwPmmqmRI8POYSoTbHTFzGOXbTgwKDiq0QIYUBQggEibQqErxiMSImtDkvvvue+mll0SKsAUHXOyEJs8//7zYuf/++2+55ZYLL7zw9NNP5yoTobgMdxDHhc2hf9l6wyECiYLmmJ2HHnooSIPmvffeA0SrUjBP4WU0POP8cNUg8H+VVVZJ7QSU5tldXUQXEGZo5lhgwrb2l1xySTJRgGCgb17T1WUAWCPcqIzoEzJAN0JUNTFuCEcMw7Q9BQqaiB2o3XDDDU8//TTsElD2TUq8y0RnnXVWs4Me+wIO7n3uuec0ZxxM03qCtEbvqa0XECwk0k36CHBYMHnyZHKYufXnB8BxdWZeDFNkpeYZfYwSO2H99dc3E37m/8cee8y0eQIvWAElNMnknWCQWMbKpFyogDba0047rfUnH5pwcL7r3377bTP8qDY7vdzemsPQxMn8oEuxqunrrrvu4osvPvXUUxUdYpv4kY/m2QvO+IyI5meUZFaphBuefPJJQWr8psqdxgyUxI55ZUY03taMMmyIGLbLL7/8cpmohH83HGDDJVs4eYEhgn0yoXuN6kyvQtHJ4lAHgck4EqIaRbyJEydKKyor2uGcuv1BsSuuuELVq8qQVgkBV5mb2DHJ6B06hCkC3FyKyQCc6lXzjY+F0tlnn61KSsvdcIROwCsgMfsQCcFQALvEmzMBxAnCKugihVQyZcoU46P5gABNs4lBNtXNcccdZzI4byQcI64jE/wPFHC0zogjPeWnEISPCSVfqpU4UjvOWcxsIYpF6NSaOAuWTgiW5AM/QZP9m2666bLLLqPbxx57rIpIdMCitfshMCNXmJ500kmmxDFiFosNj7f4jIPNonlqDYenDhaPhvKXXnqpFkzBOYs5JBycRPNzWbHAIQjBAXiF8NVXX00pucLOueeeq5WDDjoo60uh27xsyI1yH3PMMTQYMU1PNjFa48HuVjgYRKIA2BH5DxzYjWtCYTFRYMKmjR0ubiUIjPAKWNFLwQIFGuFigkwyrSyRgnw2L1h4ZuTWBMcff/x5550XiTVn4+81KS5P+AMLicQX6cEsM+JX+4tJ4JgjWEIQbGeu11AyC0KSDGF22223yXCcYOVuyaBGFB0YlJ6+Doa/ZmVudIGP6WgQiZkdryc5JNc6n5TSPussSaAbDgiJhVxMZkxPE3Y8TYxFR1NfUXIrdzoMC9rRHMXXzAxSyuNF+UFoMNMBBO8mTMzXCSKFgyUpfCfJqv5uONROaONiiAgNlss8dSUsxKErnYaEhEPO5wG8aHb+9bNwGUGMP5MynaReR7KuyT0BDhYpljbSotUQNVlMljJbSRtmWCQ/2ULH0yQRr7rS8lGVJdUfddRRqqyvMxzcSSCJHURMwVxwHBD0QjYxZ1jIjLCQGVMu7bfffpnRYnQYSWxFhFNdoBX42ScrjlNQcQGLkEqMiRSxl76/niZe0PmBBx7gTnFhOvxtxxwdNCNYTJ8+3YyUpKhhTdiEg/DSEhMGmFNxAa7IRlOQQnThhVTkSmFCNXarv6iXXr+2xv+SrsGbgomYFzPBzOj666/HC3WXMJEWqIYZ4ZQLu6tS1csFF1wgiaKJpZEsxQCksKE0XlJr0V5YSGauRLy606+1GTNdIAVyh7kwO+RPcjSjM888U9GkVpAizeiwww7LVc0lHFGQQa2LEAGuUMA3ORwppBLR5YK99tpLuaXiyCVffzNsk1IuwoUJAjMCkxlNmDDBjCKCrHlBgYO5QBmDP6BxgewNP2Ixfvx4ekF4VVyabp7dJibArSFszUXxCiBAoDlSKKYB4Ujz1Nq+gINJRWauXO3+uMEhhxxwwAH77LNPLmuFsL1MAfW7+mMmICAT9vOexhxn9CU4YioWiQqieIEv9psvtLlJujE1SPNQH5sDHIu2zZ4164P7739uiy3urP8x/O76cVe9vaP+l2iPG6rqmqq6oqoudHDttWdcdNFHr7zSvP7/J/v/ghwfPv74U0sv/fC3v/1oVT1ZVU9V1XMtj2eq6tn6oJceq6pH6scDVXVvD12ur4lyaf2//Rf98IezFhXx+EpblMnxwZ13PvH3f/9EzYCXOjpmjB79Smfnq52dr9Xb1+sdW/ses8aMmTlmjNNe7OgIV0KUh6rqvvo/Um+pquuq6vL6n0/PqKqTOzqe7vOuxCJmiyY5XtxqKzLAxy+PHv36mDFvdHa+VT/eqf/g9v36n25t32t5eNVp3Syp/0MYRV6oWRU5IST312nopqq6tqqm1xQ5s6pOHDXq0lVWafa6yNmiRo4Zu+0m6PmVDBAGLi//iv1xZ+fszs4/9jw89fio53+lnebh/DdrOaExWnixbkqD5Ofhmh+qk1trCZFoLqqqs6pqSlUdM3Lkhaus8slHHzUHsajYokOOT15/nf+6Mwi1qP8vnh7wOkJw/6f147Oeh/3POzs/6aFLWIJD+e9sFEnSoSJae76jQ13yeAs/UrReWVUX14XI1KqaVFUHNxpPL4jPYn19bBEhx6wDDyT+vCjioxayRrQBA7Ch65vf/NzjW9/CiS6Pb3/bEcc9DVdQBJM8XKiF8AM51CJSzPNV9XSPfihULW2UIPJLKUHox4SaH+cP8v9ODaUtCuR45bDDMOOFmhmSAsEgACQhtOhmwDe+0eWBEN/5Tte3vtX153/e9d3vIkf3wW98Az+Qw/l5FPHwIB7aVLjgR8TDeod4IAfxuLGqrqqqS6rqnKqaVlUnVNVRI0fu32hMXnbZz/t8CKYdre3J8d7NNz9YZ5MZY8aI9WiGHMHNyR3dzPje97r+7M+ajzADUbKtyeGBSWTGVak/lK4pPpSoWlZ5UA7FR1nlqjyQ45q6OD2/qk6vqpOq6thRow5tNPZsNK7YZ5/m+NrZ2pscH95/v3WmgCb7vWQj5OhOKMiRLdlAjr/4i25+/OVfdv2n/9T193/f9Td/0/Uf/sNn3/42cqT4iHIoWSxePDTbvb6t1z5WLoUclOPmHnLILMRDZjm+qsaPGPGHRuO3jcaL997bHGXbWhuT46k113yqo4OruI3zXqmXrJwq9CkHT4ccVKGbHLKJPIIWKPJXf9X1wx92rbVW12abdW2+edfii3f99V872VUe6JW0ojXK0X3zY/TorGnLmiVpxbJWWrm0Jodlyyk95JBZftdo2F61777NsbantSs5Xpky5emOjpfHjOEt1YCagP5zZ+5kZJFCDL4oO2xrCVGWdlPkH/+xa9llu37xi67VVuv613/97LvfLQWpFjCDCOVe2QyrlR7ZUNlIYffUNzwUpNKKBa2a47z6nkeTHI3Gfo3G7o3G4aNG7dVoPHLZZc0Rt6G1JTk+nT373r/5Gw57acyYh+qAphxCXKDzK+9KDSk7PHi9ECVaYqd75VIXHMoOR5ycO2NhhnYwQ4NKUbKhIwVHlrKyWN6FkVOsVijHRVV1blWdWlWTq+q4qjqsJseujcYfRozYt9EY/+MfNwfdhtae5Pjooxt/8ANx/OLo0ak5coNcfcCpgl5yKRTxSKXpYSdPc9wjalFSiYcWwozum2B1QgkzsFAKi2xgBtm4sqouq3OKalTBcaLVyogRhzQa+zQaOyNHTZGD/vmfm4NuQ2tXclz/gx9YVQprrhLQ4YfkMrN2rdAPRXg9t8kRxSP7eGDfDgJFKpyZtUl3kWHtqgK1fO15Qy6aIaFEM/Imi4QSZuQmWGTjiLrUkFO2azQOrMlxwDA5htiQ49of/MDCATke6ugQzTwnuXRXph0dVp6zPGpPq0L4PkTJI1TIQY+sR8obb91qUd/yCi2IEwomm9zVkk2ubqlDp9WL2DDj4EZj70Zjp5ocnu6GIsPkGGJDjhtq5eBIixShHH7wJY/y6wt1llFLevA394cBeWQ/Bz1UtbnNRXvUth4ayapVKrE2sXDNLXPLk7wxqwg9v4cZ3ZoxatSRPQmFWvym0Zg4ciRy7DJMjqG3kKO75ujo6M4gY8bcWr+rzpGFInFzcg0lsE116ZEdBMpLzknJSXuyJEG7vA2bCqPQoqSSc+pb5pihzsCMo0aODDNIxQ6NxnEjRypBLFtUHsNpZagNOW5qIYca4r1aPyIhKgN+lQv42DlcnpubtnkgTY5EJJQUTvNwSThBKjSSz3CEFvIIwbAwIRhocVq9cD2hJ5sc0GhYmIQZp4wcaWWLHEeMGEE5hgvSobYUpNwpdyCHGkKladHBowoCZQEhsY8oPO3B5RRFjrDFgDxy0CPnJHfQCZxwOU7kw4KWJNN7aKH2DC0iGFLJofXCda86g3hcO2qUS0hLIcdwWhlqi3J0r1Z6lCMrUmvU18aMIf5X1eFutcnNHgpJjkeXuD+q4GAeuRHugQ0uSb0ZTlxce/rcWi0sSQotjq5pEcFQZOzYaJwxYgS2aS3kOLlWFOQYTitDbb1qDuSIcnyU9eqYMY+MHn1hXTYiCjdzNpd7YMz19bbwwL4tMoUNHpYhrpUaUljk7XjrEUlkArVoNEILFcYeNS1OHjlS1SI9IQfO6YvGoBFy4M2wcgy1hRzdyjF6dG6ZF+UoNzAsU+/t6DizXlMIfc4W0LIDxtADD/u2eOC4R9jgZJfIHYpN0W8lMqmqjh058uiRIw+rV6pJIsoLS5KpI0e+Vr9nixwKF7mpKAcyHT68WlkolpqjNa3k5jdyoIh9R5BmpiVJR4e0ws28JZq5PI63PbXe5kEbco4tNkysP/mHEBxMJw6qb3fiRDLIdo3G5aNG5fZady81OZS9yhfkiHKEHMPKsRCstSDlHn7ySFqxjXI4PqsnprntllGj6L8SktcpwbGjRnlIE8coIOrj4+sHefDABonj9/UnM3avBYBObNtoTBg5ctbo0YWCCmHUpBzWw9bAepFWZKhScwyTYyEYclzXohyoEOVo1hw9bhPTL9RuI/jqULUFtxGJeF3RoJykB3b2rkngIV/8tpYHbLAuRQgPJzw2erReWimoF0dae+lFjijH8FJ2qK01raTmKOSIciCH4/kEV5SD4Mdt3aXiiBFIsO/f/d1d06Z5TF5ttTwO/Jd/yc6Jq6125f773zF58qXf+U56UcHoRRetvUhe3eTo7JS89NK6WtHLcFpZONar5kAFPuMwDztJK4lpgm8dcd+XqwHpAzn2/8d/bDbXj33y4Yc3fv/7WROFgmGGtIKLHjpCmvRSag76pJeiHMPkGGpDjuZqpSYHJ4UcSSsCmiOLciCHtELwv7gDMWKE3HHAPJDj2h5yIEGhoF7CjGbZ25K89GJJHHKgIHIM3+cYaisF6cs9ypG0kqWs/Vbl4LYI/helYu22eVGO62pyWDCHghGnQsG3WnqhHOU+h1UxfVKQ0qdh5Rhq61WQ9lKOxLTjvWoObmtVjgPngRzXf//7pRdtop0uCjnSS/eCue6l102w8XVaGVaOobaQ45GqmtFzEyzKEbdFOVqrgQh+yNFdKtbvesx7WkEOvYQcKFj0KWVvyFH0qZBjeCm7cAw5ro1y1OQgFXFbq+CnGkAO3u21jpj3tNJUjjqttBakdnSayqbUHFGOcoc0vQyTY6it1BwR/ChHIUeJaWnlxVrwH2pRjql1qTiPyqHmaE0rugg/Qg5H6FN3L6NHt9YcoWAK0uH7HENtyNG6WklMt9YcyOF4t9s6Oorgi2mCjxxiGjnmRTlKWunVC/55IEdRjiSvklaGlWOhWWoObssisxc5uM2R1By9yFGUg9vmRTmklW5y9FQ2vXrBQqTpJkdLL1GOcp9juCAdaps9e/b222+/00477bHHHvvWf895QM9f7DFP96v/erL8Xc0OO+zw61//+le/+tVGG220zjrr5K/01l9//WZz/Zhetqv/dOore9mt/p86vWyzzTa9/jhig3b+0n27koMndt55Z27jnl5us89tf/jDH7ht991357byxychx8orrzzv5Ci9lC7yV4a9yJG/V2ml4DA5Fo4Vcuy55559yRG3FXLsUv9LK7e1kuOnP/3pvJCDPulFO3rRZt9eKIpX6VPpBTk23HDDYXIsNOM2nuhPOVpjOuTYcccdQ44I/ryTAwWllf4oyFDQq4UcSStFOfQyTI6htqIchRxxWASfCx1xfO/6D0r7xvSA0kpqjrmQo69yDNccC9Mi+MVtRfD7VgPcVqqBklbitj+xIC0UpBytyWvzzTcfrjkWphXBn2NMF7fl16xbY3qg5EBB9Sz3I0evmqP00kpBaWWYHAvZ4rbWdUTc1ppWWgW/b80xR3L0+pOlKEd60Rpy1KxoWnrpS45ea6Jhcgy1hRytaSUOK+RwJOsIi8woR2upuNpqq22yySYOMrxh2aFGDCEwKX9+u8UWW7iQ70OOUDCGGY6EHOnFtdtuu23pJWXvMDmG2pCDL+d9tcL3qQbEN08zx8s25iVnSiJMy5gXovA3fmy55Za4omXtl8qGIY3kVcjRSsFh5Vg4VtKKamCOgs+Qo9QcnGeflhTzKmvdyR+AowvCOZ+hC6JgSUSlsMRTl+hCv63JKxQsySv3YYfJMdSGHLzFc71KxV5phb/jey50UKDbtlqOODkawJycC7VcWBItwQmMRBEZR7rZdNNNEcL5zizkwJ5CDgXpcFpZCFZqjigH19aS0bR4OpyI+x3Em16Wk1nNk26i5MJciyIREoYffN+aa5BAlsEDtYts1Z2Ket7BGVaOhWxRjqQVjizKYcupORLHH1xb2HDIIYccdthhhx56qCP5L6ocZ47k/LCkUCS5JhLSmmXCD/qhwiAh2DBu3Dh0QY7hmmMhW1EObuPFBH2SSGuKqYlxMDbgBDv88MPLzvjx47PPnJCTbQs/tBP90Kxewg/6EX5EJAo/6MeGG274i1/8AifWW2+93IcdVo6FY0U5+IzzKL8tlvQSjDj+iB5DCHbUUUcdeeSRdnIw57BckstDjogHU8wmuRR+KC+ykCn1x8Ybb7z++uvjx5prrmmpzIaVY+EYcnCPCOYnxWBSCV/GtXwsZfA3hWCocPTRR+OELUOLVnI4LfxIlkkLEQ+mzVKf6iiFp05TnLYmF+Ihm5ANmrHGGmusssoqyy677HBBuhAsysFDQpl4JL75klM5mJv5Gy34vvAjFDnmmGPCEluW404OP1wbfmiHRTyQg3LohXik7Ag5sJNySC4Wt6k8kIN4IAfxQI4VVljh5z//uRTTHHQbWluS4+OPP0YLThLKPMeFjDv5tZCDhRZ2QgLMOPbYYzFj4sSJJ5xwwpTa8nexXu2VViIbOBfl6F601JkFP7JsybI2yxbkyCIFOdZdd9211lpr9dVXX3HFFZdeeuklllhi5syZzXG3m7UlOS688EKxS+GRo1QbIQeLm8MPRj9CDiab4MfUqVPPP//8K2u7/PLLp02bhjEolctZkY2QIzfTWGqOkKPcPJVWsqZFjqQV5Fh11VUVpMsss4yyY+21126Ou92s/cjxwQcfEHPk6KUchR+YEf2QKVj4QUXshBzU4uyzz77kkksuu+yyc84558QTT3Q8V4UWsTpZNT8UEmaU1WwKUsNIQRpyJK1EORSkIceSSy4puRhzc/RtZe1Hjtdee028FuXYa6+9UpByZxEPVstH835GWMLwA0tSk8oyzI4TQotULdqJZoQZKUVbmaFrOaUwIwWH1UpWs3Si1BzSipr0Jz/5Sf4rv+2s/cjx0ksv8UqWsnymFGjlR8qFVpZwfFQhdGH2i+WcVlqEGbJVNKOUGpihR8zQdRKKUnTzzTfPUiUFB3Kk4LCOXW655ZZaaik1B3I89NBDzdG3lbUfOV599VVeoepFPLgQP7iTU7m2lSKxMKBY8+iXpSLbVBgRjLoGbd4778WMvgnFkrVVNuQUS1k5BTOGlWPoDNBcInBzk5Tn+C8UwQ/e7Y79HorYZseRUMGW5TjLyZGKMCO00CAjGGEGWpQVCmbkox6GIZswzBg3blwWsWRj7NixYUZyip3333+/Ofq2svYjB+MMes5JvMVtnBcJ4VQU4WCeZlzeywoVssWGcILw4ESrWhRaRDB0VCrQZJNWzcCMUodKKBaxWadENtDl008/bQ69rawtySG7843YpR98Zu3AhaEIpybosYSzeT0WxnQnjNocyQm22OB8hIhppNBCy5hRVq1FMFKBqjOSTbJ8jWYsv/zymGGFEmYwK6PmuNvN2pIc4pXPUgwKZZ4LRSg/j8a78XSIEvfXovCFOZJzWOrNWKRCa9SChRas3AlN+WkM6623HsFQZ6yxxhrRjKxQaMbiiy/epMYwOYbYOAYPOJjPeIvbQhGWWA9L4uxwhSVTFGslhPNZrlVyhhOSiPZTeIYWjGBYsvYSjFSglieYocLICoVZrWDJMDmG1PiGL5FDuWDLYZxHSApLCEninqe5nBXfRxiYp3k1Z9bU6r4wnJBBIhWSCIta6JelwohgrLrqqlIJwcjCtWgGiihOHRwmx1Cb2OVd5YKiMmsQHkURkd3KkrjZlsXxZVssrzqZuYolW7FIRTJIOFHUguNLHrEwQYsiGKoNTFIVURQvOTJMjiG1VnJYjjZvXxx4INd6iTsLUVicLTXE93Ycsc1OOSFlJislRdIHN5fFCE4wSaRVLVrzCIaRnJDDmU4YJsdQG7clrWRRGmbkdichoQTcyT1O4+NwhfFZq+FBjjPnMOczhFBpltzBUljIIKktsh5pTSLWrioYQ5KnQg6NSDrDyrEQjAtLzVGUI+QodzwllLiWj6UDxt+8jjTF/TH7XnUakzXCBnHPuyk26UTqTZxIyZm7W8yFxmBJbMmjmink0MswORaOldVKKzliuQeKHPvUH8JwGmerD7iqzglrcnxfy0tOY0Uh6AGRwAmEYHSChRNYstdee4WIreRQ2CIHTcI/DQ6TYyFYq3L0SiuFHHv3fKNJdaLwxCf1I2erFXid74t5mgKilQopJkriiGlBm+mOpSNjUP2EHKrdVuUYrjkWgoGey1Nz9JdWkIPPkIPPVCHKT1WFxEEYogFSQ+tOMkUpLVsNaTAgy6L0xezrRV8kqpBDWlHzUo6Qw4UuHybHkBrRjnJwTF9yFOVwAnI4k3JYkiCHEgE5cm87Tb3++uuP9NhNN93U3Hvkkeeff96r06dPdy3H66gvOdIRghZySCuFHNLKsHIsBEtaofCtaSXWlxxJKyEH5VBPtJJj7nbZZZfJFFYiyKHNVnIwHaXmQA7nlJpDWjHC4bSycEwBkbTCMf0ph4KxKEfSiqv6KsfcDTkQi+NRTUd9yaGjpBVMDTl0VNIKcgzfIR1qE5fIwR9973PwHy/2pxwWrvNBDo0kf2FDa0d6iXIYSWtaGS5IF6apOUKOOSpHyMGdreRoLUgHRA6q0yutpKPSlzHoqKSV1By5CTZMjoVgUY49699f6KscRe1LWolylLRijTrQtJKOeilHIUfraiX3OaIcw/c5FoIVcsxxKRtyJK3k7fiSVijHaqutNtCCNKsVHYUcsUIOY6BhrLUgLauV4ZpjqC1pRbzOJa2UmuM39fcWy1J29dVXH5ByZLWSjjSevph9XHEwytGXHMPKsXCsFKQc0xrNrJAjaYVyOFPdUMgx0KVsqTn6kiN9GYOOqAuJCjmG08rCtJAjyoEH8VYv5eiPHNLKUkstNY/kmD59ernPodleaaVVOVKQltvnKUiHybEQjO8liyiHfX7irb7k4DNpBT/4zCICOSxlCzmefvrpZnP921NPPbX11ltrZK+99tJmf+RQ/SStKEidTznQV82x/PLLqzlc22yu3awtyXH66acT8AQ0D4UcxRzhyKi9cwo5BPR6661XAvraa69tNjdXc6EWuL+XcrB0pCD1alQKZbepf20SOdZee+0VV1xxiSWWsNNsq92sXckhRqV5DOir9lGOkKOsMMu9KQG9zDLLzDs5tthiC17vpRytKqUjr1KOkINybFb/HMNaa6210korDZNjqK2VHMVncRjzlM9KWkkpwMdRDquVpZdeekDkKGlFs+kiFnJQjkKOpJUoB3IMK8dCsJBDyrBe5bOklRLNIUdRDj6zlC1pBTkGpBwuDDm0yWpWNE1HLGmlkCNpJcqBHD/96U+HyTGkVsiRRURfcmBMlCNpJXe1548cSSvc35pWYp7qHTl0hIg4ZFRWRtLKsHIsNOuVVkKOWCFHlAOBdqp/hafUHPNHjjmmlZAjNUdrQYocG220keJm7Nixw+QYaivkIA+9yME4zEEvhRzl9sPG9a9Bzgc5SEKUo1dHYWHf1UqUAzmGlWMh2BlnnFHSSi+f2WchR997U7lxOaCCdPP616tT3PRNKywFachhVFmtDC9lF5rNRTnsRzlKQZq0ktVKlGM+Viu0QbO9yMFCDixMQaq4KcoxXHMsHEMO/qYKAprPQo5ed0iRoxSkpeaYe0H6WW3NJz0Wcsyx5vCURTla73NYrag5hsmxcExaEaMpSLlnjsqRmqOVHGqOcoeU82gPX7Ide75UzTTrZCsOPpaJMIMM0J6Qo5dyeBpy9FrKpiBFjpVWWmnxxRcfJseQGuUIOeaeVnopx6abbuoqesORDHWyE3OmlziYTrgEITDG+dv2/M6TFpCg2U1tyKGjQg4XuiT3OUpaGb7PMdSWmoMvMaBEc3/3Ofibs4W+HOR8x7nTTjFHYnyc6oFpPFxxLWkJSzgeRVDNVemrV1qJcrSmlWHlGGor5ODRklYKOVghhy1Heuo0pIk5ofVpXuXjYmjE38lKoQghCUWScSxhMABdXOjyVnL0un0+TI6htta0wrucjRYx+/F3OJFXcwL2hEBlp1jOCWNcVShCS5J9pIxQROJIUSLRIIHywtZIck7SSmvNMZxWhtpynyPk4NHiezuhRTlYfpU2O4f1/HsGy/GY/ZokX3ArFNF+JCSJJuUICSn8ICFKmdxeIyQo20s5hlcrQ21JK72Ugy9LBnGkeB0PcCI/bJ0tK/+24Qi6dBOn5kr4wTSSBsOPlC+tVUhKkOiHRTKpsBRad9110SL3OXKHdDitDLUlrXAV5yFHUoBtq2DYlh87R4WY/fzNyvjaHOlmzZz4oamIB8OPFKqFH8iBneEH8cAG/ECIcePGocKqq65qZ3i1snCs9fY5n4nsUmFwLQdTi1Zm5Afw0eLo2uyzcCWnsSQal4ccLPlFyywlJyvFhwEgaMQjyeWX9V9qEI816n9qWqG24bQy1FZWK0xBgCLiu+QXPuby0CI8QAv7xx57LJYgR7bdf9VUHw9FWHce+nJmYUU5Ih59K1OL21J5bFD/BKWEsnrPn/EMp5WhNsrBMXI/necz/uPFMIN3kSOaweQOXrdFke5/0Kgl5Ljjjps0adLkyZNPOumkE0888fjjj8cVp7WmlchGLORoLUujHEkrZWWb926QAxvww1Jl2fpv3obJMaSWmkMEI4fVJnJwYatyhB/F6gKj+38hbfEDLc4666zLL788/9R05plnokuUw+Wt5NCynKKgiWzoLgsW5DCAohxZ0xblUG1QjrFjxy633HJLLrmk5PLxxx83h95W1pbkmDZtmsANOVQDdXg3b3MVcrDuVNHnn5qIBNlAiEsvvRQzLrjgglNOOQVj1BwhR51Mug0zGNnI3fFSjSJHubNOObaof7IyNQdypBRdbbXVQo58geq5555rDr2trC3JwR/cw0l8JqYFNy/iR5QjFn6kzCyrWYYluEJFqMWECRNKQnFJtActwrZoRlnHsiIbWcqSDSNhWa1YzWJG/ljDgkVaiXIgh76aQ28raz9yfPLJJ9K8wOUkodzKD36Ng8MSShDj+2xjORiridTMI4UWjGD0LTXIBlpENjDDMDAj1UbJKSqMslpRkKo5kEPLzdG3lbWlcpBxgVsyCxcKcb6MfoQiYUl8X0jAWvdZTrMNJ7SAZ9GMwoxohr7KIiUJxSKWZiShZB2LGSWnLL/88nKKpeywcgypiWOO4STe4rOiH60UYXF5N03qHaQJD5j9Yl4KJ2qx6P5PlnCCFcFIKkHHFKGYkTo0mrHhhhvmDgdmFNnITxxjBrvwwgubQ28raydy/PGPf7zjjju4RygjB1W3L5Q5jznInbntwSIAvN5q3ZSpqdCXEKktXG4bqQgtUmSgRSpQzNB1ySYpQjGjLFIwQ7WxTP1/PIvXv2GKInjj+MMPP9ycSZtYG5Djww8/POOMM3iI7zmSX5Mv0EIEc5idSEihSFhilRHH1wTo3sGGPLX1ajjhTJarWEkiaBFmlFVrUklZuLK8pdJaaqy44opZpKTaWHbZZcMqSuN8J5OZSy65pDm3r7d9fckxc+bMiRMn8o1Q5k65ACGUk1md1ncuxnM2/eA2DghFuJNrpQPWSpS43zbykKcxp5X0YRtCJIm00iJ3QjEjN0P5mJGELE/CjKIZYQYtKVTTYGEJbtEbl5x44onvvPNOc8JfP/vakeP111/nZh4S1pRfTUAnyu0KhhPZskmTJp188sngjoQky4Qi3TJSW7gSBsQcycHsO8f5nBdOsNYkErVAi9SeUQu0oBYRjCxcU2dghlIjzHB5LVjN2yRMd/rSi0F6VbPalJhwyKtvvfVWE4KvjX1dyPHqq6/ef//9Tz31FMF45plnik7kzkRuTtgixJFHHingTq/tpJNOchB7EtaFJSwOtuXsMMaOLQsVcpzlzLAhphGCVGqLsiQJLcbV/8oTWliY5O8grU1SgaozGB7U5U13ZROKEDksiW6FqUaid/166brrrpNrTj311ClTpjz++ONNUBa2LUxyfPbZZy+++OLdd9/9/PPPv/LKK0JHeTF79uxPP/30k08+sfxDhYhE/R7ZUThx5pln4sQJJ5xAS8oSFPQQ504u5Mvu21JbbOEpHzPRz+J1nqAu2S8WeXA+K7kjnMhiJFVn3jeRRwotSoWRVJJVq51QAS1YKMLy1HGvlloHS0zqiSeeEBj33Xff7bfffv3111988cVnnXXWhAkTVN9waIK1MGzhkAMcDzzwwLPPPvvaa6+1cuLz2pondXVdeeWVxx577OTJk3HitNNOs3/IIYdkIQrocCJJXRQSAB4V3NzJqSEK42yGLpQghgHMjoN5lTnTJbFwIjpRFiNZqbJCC2qRPBLBwAzHjceoUiexsKRJkB4LRdgtt9zycm3C4+mnn6YZDz30EJbceeedN9100+WXX64SN/1LL720icjQ2tCR46OPPsrkX3rpJYWFQuyDDz6wOsUJEtLKiVaTaOQRnKAT4YSwg74kjRMSBE4IfUHP02I9ss+dnGpf0IcrhTEskhASsGgDwwbmwlZC5L5WpIKNHTs2arH00kuXwjPmNKkq5Y6xFZYQiVChWzpaiPLkk0+KDfmUzZo1a8aMGVgiYJiXHn30USwhq7fddttVV111/vnny6HyzlCWJoNOjvfff59IkArzf/PNN8MJOlF/u6xfTjAvUZR33333ueeeC7itOpHsIPpFPAfjAXdyJP9ZMtjWQb4W7/JxqgRbqYHxvW04ZMeWOSHW68+aNKWqyB1P61KcaJWKYlhlJAgqQ+ErM0gsUVvUxWj3PbrCEsqBECIEILZvvPGGLQtLUES2pSVqLywRUZIOltx6663XXnutpDNt2rTjjz9eIm4iNWg2WOQwT5wQ9+HEe++9hxAyaHfa6J8QzDnYg0Ou0giZQQ4AAVdc4gT0JQWeEPE4wf1cmD9qFMrFZ1F+x+1wsHM4O8b3cb+tp+EBW3311Z3JXOJCJNNmcgdO5L95Uli0WvKLBYvWUI0CUSYsoWe9Vk9ojSWqKJNSd/OuCRIP5DBZ27DEQS9FS1gQgGR0F0vUIjfccIMCVmlisUZamtjNyWhzc2/g9tXkEN86QOF77rmHvxm5Y/fee+9dd92FzgZ69dVXY/Rjjz1mMk5DedMzYaEvm0gcrNncnAxddEFj3n77bQABDiLUVRrW5oMPPhh1lVDgjhNCPxUAL3IMTsRPPIc6zUZrI8JowbtOYxzJeD1GD7LjIB7EQgWm2bChl0gwB5HV7Jrd1EZydIFbWEKKsETmipYgiloYuUU8N8PH7GhDWIIH4EKL8CNyYos3XmLOAWzkBCZkOKUJqFUtV1xxxXnnnaew5ZFHHnkEYtEbp9lKTw8//LCTU8dcc801l1122Y033ghqnm0OvR/rlxwff/yxPljmYHCGyEwjA33hhRe85ATd85+BhiiOC32XI8RX6gQ5wQk6AQiNazZRok0sDCcU8OZ/wQUXCBSpWogL67wbXtxmaaOdZqP9GNQkJn5lvK5csC3mqeMa7EuFmFfVhubebK4fc4KFFWJhnqHKU+RN4lPzmoL1CBdyLS/yHE+bLwvCsGVFS7ozTU/SsS0skXeAjGGA0k6iFHSujRQxl9gGUo0nQzkZqk6++eab1TEXXXSRgpfXmkPvY73JgU2616sW9STyeI6BnqkAJAhbxzNcJxuoOCh0Nn+XN5ubk+kinEBecyicABaOI5nR44T8mg/jWKpY/XP/4YcfriZVe8ZbKCIImo3+CWZNFJs+fboGsw++5st/glleRagICf0/55xztE9lBS6UosQoy6AHQzhwPEgThKwwIxaXOx6ipDoBHSTjF2HZahzHUy5xfiiVwBPMhMQcg7CxKXWd3xx3j32JHJrAZb0ah/6cLSlwZGRANWCH/tt6KVyJ+uk4ioebqicd932TCSe05irDNclEQKG/q6RSugc7CJ577rlIrT7HiUMPPZSGS9i5SyF5a2HugvS1MhNHAtWoNdeECRN44uyzz5aFxS6hLSwBApbAkAEzEc8jLEoQLSksYdEGqPJISrpWc0TiQxquzLVYFcnRUYQZP0TChRdeeNppp02cOLHXmvkLciCvK/WHblrUnw7igzo5NM1UESUU4elwU69hcbIM/RD6zOCcbFsoHBohRMoruEAnnFBhhRPS59FHH33wwQcDNEWohQlJFyIZavuauVuZW9CaoGmefvrpHCNpYgl5J5mgAwvnQRJLTDksiZaEEyFHtuEEvzQ7qC1ucrzVTfghFWiHl7XMTQlIXZNMS+WpU6firuE1Wwk5dIC2YUYEw5FUDDmp1dIx02shpoGGlckvKSHVSliZ9GTrVS9F01LJ4kRoS28xN5zYf//96YQiHyfYCSec0Ox4kTPzxZLybgB5J/IKRrjxmUoiSQdoiIIliCW0eDdywtncPxc3eYmbuJK0FDe5kC9KGcBNHEHDcluWVI8fP17h72SNdJPDIFBJr2GGtjBujl0yx1lRrZIm9Jo0YSbp1STVDZTKTniaPIeqKTCpK7Yq4AGk0pQ4dqg/TSNxEBJdNLtcpI0LobTvvvuqqAj7tGnTyCcRxRKlCZbQkiQdGMIWwsKMpwTnXNwU03hxUwoAXkYOjbSSA9rcxB2nnHKKOkl8MpcvJqzVB/pL+akVzXEM0s2lY70aHCYZZfJFekUy+SL5DCfM0DwJg8iw6FBMkK9SYFrxC51STOy6666AaHbw/6WZfvdHFw88MKUJ3PiMuIJOjia3RBdLoJ28jxxRjv48xY9hRsn+Ja1oh5si4dxEOXSHHPxiALvttht1WSxST7KEvuvVmHplKFJnj26K9O3bsEomM9bUHKQvrSWZIb6JWS/RK2xASSKh78MOO0yg7Fz/4E7ek3RVs91hq03IiV1LM6BJOsIppQn28C7MuYnXs1wokZxts4m6EOZHNOKm5H0uVsHwumyFHDKXgu+6667L7XnVHkYec8wxuhaxYnsx2sJcwLvEQ8faIgn61jQS6N42O3hqQF4yuFZmpMtUo+GjpiMV5513nmWS0iG5I5xgJFR3zXkMWz/GwUJLYa4aI8YiPm7ibPhzfChSh/MXxlNZTjonSd8l0YzkfcwQwNwky5flYe7Kk/N99tlHzYeLi6U85FRxn3pHx5EQfaePmH2dITVOMP1lRR5myCZRqiQUfFRe4WM+dSECCIZKh1RQlObUh22eLRUbN/EuwCEvLHmBO3gqDorLio9wQtwybkULPkrSLzdMMUNCkbkEsAWB9SAKyil77rmn6O1WDpdZNVhScq2L5RelA/0JPYulMwNy3Ku5S2GUxpr7KuQnt+p0qcKyOsr9K0qV+1dJJRalumjOeNjm2ThFvPEoNwE8ZSkvcEeEId6JcRCLg5iwz72DRC9fq2AsC1KH0owwI1/7o1IpAbWzGBVyBn6IeA7mZhTBMi3q3piMoFg6izo5jUDldgpaGLfSJglMg9EMXbYWnrrEzeZ0h22ABtXcNBP3PMXZXM5Tgi1OseUvW75jySDxkfNdxcVokbtKolcdY4WSbELay40l5PCqHruXsoSIL11ASfiYsxmW6V6paAQ6sM8QQk9MMtKf09DQJS6ULEILgqG0sUZVSR111FG6lMOiGUxuqmc6bAM2YUyJrf5EIB/zdFgiTTBRzS98lB3u4yMOCicSuqrAvClhlZBUknvQoveggw7afffd0YKPPFXr6LF5h1SWKvd0uVlqsNyQIzKCGHnwNGzAX53hYNSp9IcWFia6tBjTh7pXl7l7IZP1XfUM20BN+PFU3qPBEp5CFMYp/BXLPgcJWj5ymtrCJfxLBVS4lgh8NGnSpBK95YYCOWn2VMgR05/LXCz6CQuiaJShW0iXbpxmjZrOnGl97JLQMP0deuihaIEN5VM5Lmz2MWx/solSqz9Zm6fhjy6Ck0f4K66xw3iQgwRtfER1ZBCh69rcJj/ggAOk+1133ZWPGHKoM5p91PYlcjDlbkpIeSFE0aiahSrEdIO5Dhqcc5wsfThfMWE9otZNf+X+N9JYWTVbH7YFZ7IGBysXxCRPMbLNL1xmJ4YNfOQEhJg4caKSMxmEj8qHpxhnWXU2222x3uQoJrnkawFGoFFcMw5mhzmCDePHj3cCXfpD/UOtOBFpIhVkw+ibbQ3bYJqIJdUWGnEW5Y5r7KOCg14KIThlt912K5yQSuZ+W6FfcsQ+/fRTyyTSlM9C5hOz9Wciu7+io8wsbGC//e1v5TmL3ubFwza09u6778o4MrsijxiUbajA+Ii0WNGobZvXzMW6uv4fDbAfmDDX0M4AAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "Scientists can learn so much more by analyzing the direct physical obervables obtained in a laboratory setting.  Microwave spectroscopy is an instrumental technique that can determine the rotational constants of small molecules.  From these experimental rotational constants, the Kraitchman equations can be used to determine the molecular structure to high precision of a reference molecule and isotope substituted isotopologs.\n",
    "\n",
    "$$z_i=\\sqrt{\\left ( \\frac{\\hbar}{4\\pi} \\right )\n",
    "\\left ( \\frac{1}{\\Delta m_i}+\\frac{1}{M} \\right )\n",
    "\\left ( \\frac{A^{\\circ}}{B^{\\circ}} \\right )\n",
    "\\left ( 1-\\frac{B^{\\circ}}{A_i} \\right )\n",
    "\\left ( 1-\\frac{B^{\\circ}}{B_i} \\right )\n",
    "\\left ( \\frac{1}{B^{\\circ}-A^{\\circ}} \\right ) }$$\n",
    "\n",
    "$$y_i=\\sqrt{\\left ( \\frac{\\hbar}{4\\pi} \\right )\n",
    "\\left ( \\frac{1}{\\Delta m_i}+\\frac{1}{M} \\right )\n",
    "\\left ( \\frac{B^{\\circ}}{A^{\\circ}} \\right )\n",
    "\\left ( 1-\\frac{A^{\\circ}}{A_i} \\right )\n",
    "\\left ( 1-\\frac{A^{\\circ}}{B_i} \\right )\n",
    "\\left ( \\frac{1}{A^{\\circ}-B^{\\circ}} \\right ) }$$\n",
    "\n",
    "Here $z_i$ and $y_i$ give the positive displacement of atom *i* from the center of mass in meters. The other\n",
    "variables are defined as follows:\n",
    "\n",
    "$\\Delta m_i \\Rightarrow $ difference of isotope mass from reference atom mass in *kg*\n",
    "\n",
    "$M \\Rightarrow $ total mass of reference molecule in *kg*\n",
    "\n",
    "$A^{\\circ} \\Rightarrow $ rotational constant A of reference molecule in $m^{-1}$\n",
    "\n",
    "$B^{\\circ} \\Rightarrow $ rotational constant B of reference molecule in $m^{-1}$\n",
    "\n",
    "$A_i \\Rightarrow $ rotational constant A of isotopolog in $m^{-1}$\n",
    "\n",
    "$B_i \\Rightarrow $ rotational constant B of isotopolog in $m^{-1}$\n",
    "\n",
    "We will use these equations to determine the precise molecular structure of ethenone with <i>C</i><sub>2v</sub> symmetry pictured below.\n",
    "\n",
    "<div><img src=\"attachment:ethenone.png\" width=\"100\"></div>\n",
    "\n",
    "Consider the molecule to be oriented in the *yz* plane with the *z* axis along the principle axis of rotation.  This automatically means the <i>x</i> coordinate for each atom is 0.  Assume the oxygen is displaced in the positive *z* direction.\n",
    "\n",
    "The following table contains experimental high-resolution rotational constants in MHz for isotopologs for ethenone.\n",
    "\n",
    "|  | A | B | C |\n",
    "| :-: | :-: | :-: | :-: |\n",
    "| H<sub>2</sub>CCO | 282101.185 | 10293.32117 | 9915.90548 |\n",
    "| H<sup>2</sup>HCCO | 194305 | 9647.0664 | 9174.6457 |\n",
    "| H<sub>2</sub><sup>13</sup>CCO | 282112 | 9960.9659 | 9607.1276 |\n",
    "| H<sub>2</sub>CC<sup>18</sup>O | 287350 | 9761.2368 | 9421.1236 |\n",
    "\n",
    "The mass of the atoms and their isotopes are needed.  The reference molecule (first row of the above table) contains the most abundant isotope of each atom type.  The isotope masses are listed here.  Source: [NIST](https://www.nist.gov/pml/atomic-weights-and-isotopic-compositions-relative-atomic-masses)\n",
    "\n",
    "|  | mass in *Daltons* |\n",
    "| :-: | :-: |\n",
    "| <sup>1</sup>H | 1.00782503223 |\n",
    "| <sup>2</sup>H | 2.01410177812 |\n",
    "| <sup>12</sup>C | 12. |\n",
    "| <sup>13</sup>C | 13.00335483507 |\n",
    "| <sup>16</sup>O | 15.99491461957 |\n",
    "| <sup>18</sup>O | 17.99915961286 |\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ==> Import Psi4 & NumPy <==\n",
    "import psi4\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1 - Using the Kraitchman equations\n",
    "Use the Kraitchman equations to find the displacement in the *z* and *y* directions of each atom except the central carbon (one connected to oxygen).  Modify the code below as needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define all constants and values needed for Kraitchman equations and convert to SI units\n",
    "h_bar = 1.054571817e-34 # Planck's constant in J*s\n",
    "RotConst_MHz = np.array([[282101.185, 10293.32117, 9915.90548],\n",
    "                         [194305,9647.0664, 9174.6457],\n",
    "                         [282112,9960.9659, 9607.1276],\n",
    "                         [287350, 9761.2368, 9421.1236]])\n",
    "\n",
    "# Convert rotational constants to SI units\n",
    "RotConst_Hz = \n",
    "\n",
    "# Define reference rotational constants for easier use\n",
    "RotConst_A_ref = RotConst_Hz[0][0]\n",
    "RotConst_B_ref = \n",
    "\n",
    "# Define mass constants\n",
    "amu_to_kg = \n",
    "M_ref = np.array([1.00782503223,12.,15.99491461957])\n",
    "M_iso = np.array([2.01410177812,13.00335483507,17.99915961286])\n",
    "\n",
    "# Determine total mass of reference molecule in kg\n",
    "M_total_ref =\n",
    "\n",
    "# Determine mass difference from isotopolog and reference molecule\n",
    "M_delta = (M_iso - M_ref) * amu_to_kg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define functions for Kraitchman equation\n",
    "# Each k is one of the 6 terms multiplied together within the square root in the equation provided above\n",
    "def kraitchman_z(M_delta, RotConst_A, RotConst_B):\n",
    "    k1 = \n",
    "    k2 = \n",
    "    k3 = \n",
    "    k4 = \n",
    "    k5 = \n",
    "    k6 = \n",
    "    z = np.sqrt(k1 * k2 * k3 * k4 * k5 * k6)\n",
    "    z *=  # convert from m to Angstroms\n",
    "    return z\n",
    "\n",
    "def kraitchman_y(M_delta, RotConst_A, RotConst_B):\n",
    "    k1 = \n",
    "    k2 = \n",
    "    k3 = \n",
    "    k4 = \n",
    "    k5 = \n",
    "    k6 = \n",
    "    y = np.sqrt(k1 * k2 * k3 * k4 * k5 * k6 + 0j).real # to avoid errors when computing value that should be exactly 0\n",
    "    y *=  # convert from m to Angstroms\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Use defined functions to solve Kraitchman equations and print results\n",
    "# Since the Kraitchman equations only return positive displacements we have to manually assign +/- signs\n",
    "zcoord_H = kraitchman_z(M_delta[0], RotConst_Hz[1][0], RotConst_Hz[1][1])\n",
    "ycoord_H = \n",
    "zcoord_C = \n",
    "ycoord_C = \n",
    "zcoord_O = \n",
    "ycoord_O = \n",
    "\n",
    "# Build xyz coordinates and replace y & z zeros as needed\n",
    "xyz_coords = np.array([[0.0, ycoord_H, 0.0],\n",
    "                      [0.0, 0.0, 0.0],\n",
    "                      [0.0, 0.0, 0.0],\n",
    "                      [0.0, 0.0, 0.0], # temporarily set as zero\n",
    "                      [0.0, 0.0, 0.0]])\n",
    "print(\"Incomplete xyz coordinates in Angstroms\")\n",
    "print(xyz_coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Part 2 - Use center of mass to find xyz of final carbon atom\n",
    "Now that the xyz coordinates have been found for four of the five atoms, that information can be used to find the coordinates of the final carbon atom.  The Kraitchman equations assume that the origin is located at the center of mass.  So the center of mass equation can be used to find the missing information.\n",
    "\n",
    "$$y=0=\\frac{\\sum_{i=1}^{N}m_i y_i}{total mass} $$\n",
    "\n",
    "$$z=0=\\frac{\\sum_{i=1}^{N}m_i z_i}{total mass} $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Rearrange and solve center of mass equation for coordinates of final central carbon atom\n",
    "# Remember to manually assign +/- signs\n",
    "\n",
    "zcoord_C2 = \n",
    "\n",
    "xyz_coords[3,2] = zcoord_C2\n",
    "print(\"Final xyz coordinates in Angstroms\\n\")\n",
    "print(xyz_coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Part 3 - Bond lengths and bond angles\n",
    "Compute all bond lengths and the H-C-H bond angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# bond lengths\n",
    "rCO = \n",
    "rCC = \n",
    "rCH = \n",
    "# H-C-H bond angle\n",
    "aHCH = \n",
    "\n",
    "# Print results\n",
    "print(f\"rCO = {rCO:.4f} Å\")\n",
    "print()\n",
    "print()\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Part 4 - Calculate Rotational Constants with Psi4\n",
    "Use Psi4 to determine the rotational constants from your determined xyz geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Define molecule geometry in angstroms\n",
    "# Example: \"H 0.0 1.23 4.56\" for each atom\n",
    "\n",
    "mol = psi4.geometry(\"\"\"\n",
    "H   \n",
    "H  \n",
    "C  \n",
    "C  \n",
    "O  \n",
    "\"\"\")\n",
    "\n",
    "# Set each atoms isotope mass\n",
    "mol.set_mass(0,1.00782503223)\n",
    "mol.set_mass()\n",
    "mol.set_mass()\n",
    "mol.set_mass()\n",
    "mol.set_mass()\n",
    "\n",
    "# Find rotational constants\n",
    "rot_conts_wn = psi4.core.Molecule.rotational_constants(mol).to_array() # in cm^-1\n",
    "print(rot_conts_wn,\" in cm^-1\")\n",
    "\n",
    "# convert to MHz\n",
    "c =  # speed of light in m/s\n",
    "rot_conts_MHz = \n",
    "print(rot_conts_MHz,\" in MHz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Part 5 - Analyze\n",
    "Answer the following.\n",
    "\n",
    "1. Compare the results from Psi4 to the provided experimental literature values.\n",
    "2. Why does the geometry/structure of a molecule and the rotational constants have a relationship?  Think about how a rotational constant is defined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your answers below\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}