{
 "cells": [
  {
   "attachments": {
    "3D_Spherical.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAALcAAACpCAYAAABpsqx0AAAABmJLR0QA/wD/AP+gvaeTAAAACXBIWXMAAA3XAAAN1wFCKJt4AAAAB3RJTUUH4goPFSAdaeQxPgAAIABJREFUeNrsvXl4HMW1PvxWdfdsWm3J+77I+8KWELLfS27W+0u4BDsQLpCAQxJyjbGNrRn1zPT0TPeM5AVMDAkXArkEAsSGJGSB5EtIIAECIQl43+VdXmRbsixptu6q74/pIUKxZmqMAWN0nkfP4xmfqqmuevvUqXNOnSPh/UujAXwPQBrAvwN4CsAeAFvQT/30HqdJAGY6//4NgP8PgLt/WvrpfKIlAA4CGNw/Ff10PlEdgG4Ac/unop/OJ1IAvAzgAeezBKC8f1r66Xyg5QD2AvA6n6929PB+Ok+IvE+fexyAPwLIAtgNoAJAJYAPATjVD4t+cL/Xn9vdaw4YcmbBfuqn84KoKGNTU5Pn3nvvFeLVdZ2Ew2FPKQOJx+PCZshQKEQMwxDuf/ny5VixYoVLlP+WW26hoVBIeDyGYVBN0+QSxk/D4bAiyh8Oh6VS+i95cc9Tull0Djjn00+ePCm003V3d1NJkupKGQjnfHyJh+EposzJZBJdXV3DRPk9Ho+iKMoE4W2QELfL5Rokyl9ZWel2uVzC45FleQCAgfnPmqb1qyV9kaZp0HUdLpfrIb/fv8LtdrMiTWxK6QzG2FYAXOAnJAATAGwVfHlsQkgd57xZFNyEkPGc862ia0gpHc4YOyjCTwhxARgpOh5CiJsQUsUYOyLYvxdANef8kGD/1QBkSumxkydPbk8kElmhl+L9CG6XyyXff//9dbfffnu5bdsdjr5dEHyMsU4AHSLg5pxLnPNOSulJB+gi4O7knHeIjJ8x5pJluZMx1kEIEQK3bdvVhJAOEfBxzj0A8s8LQX7Z4ScC47dlWZacuYegsHDZtn0qnU5b/dr0acg0TUSjUa9pmjOcr74vqpY0NjZevGzZMiEg1dfXS5qmzSplbLFYbGoJOqtiGMZsUX5d16Hr+hhR/oULF7p1XZ9Wwrx6TdMUVjOWL1/ujUajwuOJRqM1mqYNLlUtOa917hdffPFNny3LqqSUTuju7t76vjINEVIyf6lt3knSdV1MVz+fF/UjH/kIdF0f0d7e3lpWVlZJKR0WDAY39O9h7w86760lkiRVDBgwYIKiKOWSJG04deoU5syZ07/y/eA+Dx6Q0qmMsWOapu1Jp9OoqKjA2rVr33cLLcty/iBc/VbO4gA+jlz8e9lZHuKgfnAXP/gBAK655hpqmuZU27b32bbdWYqudr7RwoUL8eKLLw4BsArAjDPspgzAzwH8G4DPIBf/Lp3FYX4FgHlWX+jzbSH9fj+amppq0un0MNu2tymKUpdOv7+96hUVFRdIkmQCmA9g25lodwB+CmA7gLyE2AXgUgAvnaVh3u38LQewFGL+hPNXcp/OJBQOhwdlMpmqrq6uzeFwOIt+GgHgiaFDh953hsAGgLEAPgXgoR7ftSEXaFYyXXLJJWzAgAF2HwD/OoCa973k1nUdixYtIsOGDXMvWbIkpWnaBJ/Px+rr65v7Mf0GRQFsv/766/+eSCTy6sWnkHPSPAvgiwCmArgDuSjJ09EVDpjX9/iuvET9fQByXtu/f/rTn77O4/GMAbASQHsPnmbkHGqXAvj1+17nrq2trchms6NisdhYAF319fW783q3IIluf3zp0qVCvG1tbYxzXuq2Ksxv2zYX5B8K4KsA7ho6dGheN0sD+DKA7wL4IYDPOjxvSNJUKsV7jf9zjtSvBjDKAelQAC3O/zMU9vJ+FcCvAFgA5kuSxFwu1wgAz+Gf8fRwxpAF0DvOhjPGOABwzs//2JL58+dj9erVMAyjmlL67wB+19DQ8C+x2IlEYlpXV9duwzCSeVUmf7B0u93/d/vtt2s+n88uhicAFwDYCHH3+yRK6SZBAcI455MJITtEnp0x5uKcT4rFYoNuv/32P/t8PjuRSHxwwYIFfy0vL+cAIEkSdF2/LpVKha666qoPXnTRRRWMsX3l5eV8yZIl91mWdcEHPvCBG77zne9sjMfjn/n617/+B8t6w7PtcUC848iRI67vfve7L5aXl3e5XK42AMSyLG9HR8e0iRMnfvW22257/uTJkz7G2ABKaUtvTP3yl78c+fLLLz9TV1d368KFC3//P//zP2tvuOGGBdu3b7/0xRdf/O4HP/jBT/zXf/3XXsYYNm7cWPHYY4+9MnDgwMa77rrrwX379uU1i4GUUpkx1ipJ0gG/38/Pa8m9evVq3HrrrW7O+fRsNruuoaHhVDweF1JlQqGQ+/7775+STqeT27ZtOwTgcKE/VVUPSZJ0vBhf/i+bzR62LOvYoEGDDonwu1yuQ5zzY4QQIX7btg8RQlqz2SxdtmyZX9f1G5PJ5IR4PH6AEHIYwOEtW7YcSaVSH+WcN8uyvD2bzR6RJOnwggULWrLZ7ATO+R0rVqx47uDBg8fnzp37qG3bb/Tf2dl5yLbtY7ZtH37kkUdSACqSyWTg+PHjH5w3b96HkslkDIDd1tb2+/Ly8iMAjlBK/2XsnPPDr7766uUA0rZt/zwQCLgZY25Zlps7Ojq6AWDjxo0namtrDwE4/NRTT8kAyjs6Ov7x3//930fy/TDGjqTT6cOEkBa3281xvlIsFssfJmtM05xlmmZlPB7vMzwzkUhMCwaDXtM08+1rDcPIx33ci/d+bMkaZ0dx95KaiqNK/NR5qfOxHC4AxwF8rK/+Fy1a5I5Go/nYkpEAunocHokzb1vyv1cgtoQAeBTAP5x5vhzA2gcffNAzZsyYsKP398w6cJWzS05806GhV2zJealzP/300wiFQojFYhM8Hk+1qqrrCSFCb7KqqjBNczKAymAwuL5UPfdcI845AfA/Dkj/AqCp1/PIyF2f23PixImeTcc4L8LG05j7TkenkMsQYPV4aT4L4MH877W2tpJMJkP6OEekeoxrNoCXv/71r2fb29v/wzEj9hzchxxpve9szNE5D+68pJ49ezY+//nPIx6PX5rNZjsCgcAu0T4sy2KGYcxkjHWFQqHzwpLCGCMApjugvhfAL3utJ3E+tz/99NNvssQ5ADrZ47t5yCUmOp2X8CSAPyHnmQSA25C7d7oqj6G77rrrXtM0/X0M9QkAk5FLozERwO8lSfrSyZMnpwD4To+XBs6L+mwBq835Be5QKATTNCd86UtfchmGMYkx1hyNRltvvvlmUQknl5WVTZIk6VAoFDpwvsSVODvWUwCedKTfs70sFlln21fKy9+UsWIcgN/14h3omAdn9vFz33Gk7v8iZze/qgcACed8mGVZ38DpU2M8DeBGAEEAnwbwTUrp56dOnfpfAHoeoCc5VphlZ2tHPafB3cMaNUiW5emc8wOc81YAuO+++/ps94lPfAIAEI/HFQB1brf7oGVZxwCcN3EllFLuSNv9jiTtTVnk3OWzvvSlL/X8Pg7g27147wDwZ+RyJZ6ODgO4AcA3ASxwdPaelqSbALxMCOnso/0aAF8DsAnALffee+//XHPNNft77TKNAO46jbp0/oFb0zQQQhCNRgdTSofLsrwpFAp1h0Khgu3mzJmD559/HpFIZLBt29MJIbtOnjyZLNbuvUTt7e1EURQR6RYBUEMprauoqOjL3FnrqDZ3IedEKZU+IMuyPmnSpPoipv3avCrU1tbWWz9vdF6gxrMqAM5VYOu6Dk3TKgkhozjn2z0eT0ZEP1+7di1uv/12n6IoI0Kh0OsAMudy4P2Z0Pe///3rduzYUdYDNAP6YO0EcCPnvOmPf/zj0D54kgDucdSbM6FDl1566ZKbb7652CFwEoCdvb6TkfOg7gFwy9nStc9ZcCcSCei6jmg0Oq6srMwXCoX+zhhTRBx+oVAI0Wh0QllZWZmqqq/h/KRxlNLq8ePHd/ewfhS6qb7d7XZ/ta6urq2P/+86Q4mdpwNf+MIXutLpdDEJ8rKzQ6CXShNH7rrf2VfdzrWVCwQCME1zOmPs1NKlSw8L6uYcAKLR6DRCSIeu662iLtr3IH24urr6CcZYzzUsCKxUKpWqqqp62y7WCkYaZPGvSY/ypsK3hc6JwKmFCxe6BwwYkE2lUrLP55uVyWQ26rou+tCcEOI2TfOSVCq1JRaLdfV0sZ/NNRQNnX2b7yD+eNGiRUo2m80HLbUUE1JLlixBRUVFKYfVksFdyvPatl3y5PTsX3R9zwllNBaL1SaTScvj8QyTZXl7Q0OD3fMBDMO40Ov1vrZ48eLTqTFDs9nsJLfb/VJ9fb01Z86cN1lEeseWOP1BVVV3dXX1IwsXLlwlcDizCSGzOeebIBZbQgFMIIRsE9wdGSFkAud8dyHrCCEEtm0TzrksSdIEzvlWzjkhhBSTnoQQMphzfhhieUVkzvkY5GK2RcgtSZLXtu12QX4XIaSWc94iyF8lSRJnjLUlk8mNsVis65zfX/M2Z8Mwpuu6PtLv95O806YnGYZx4cqVK9/4PH/+/NwRu7FxqGma06LR6Li+FjfvfgeApqamfLvBDz/88Fjk4oeFqLGxcXYpzxYOh6eV+IIXyzB7CXLBTLjllluoYRgl3aiJRCIjRXm/+tWvyrquC2e8NU1TMQxD+JrYqlWrXLqujxDl13W9OhgM1vSU3Oe8zj1ixAii6/poxli12+0+4na7uYjJbtasWUTTtBFOopl9kiRRkW1x7969ME1zRCaTGXjdddftQc7lLDoHcimxJYSQUlW+Yle2/h9yNm3U1NRIPfhHOYfKPikajYIQInwlbOjQoaWOXy6FP5vNltQ/IUSSJEnqAfZzV+c2DAOUUmrb9jRZlg/btr2NMVZ00L/61a/wwgsvkEOHDk3zeDzHAoHA1ng8XlEM2PmsTEOHDp1h23ZrOBx+r+Ut+XcAr/Txf0McVWMv+undl9zJZNJj2/Y4Qsj+QCBwzAkCKkhjx47FCy+84PZ6vXWMsb2BQOBICW++Eo1GpyuKsjscDh95D1pSXkDOjX1On53el+BubGzE8uXL8/qc5PV6p4ZCoV3BYPCkIDj5qlWraEVFxcWapm2PRCKdxdrkAWzbNnG73ZeGw+FNtm13lbKtnUNUyIG1E/3lBd89cGcyGZdt29WxWKx28uTJM9xu97p8fLUIQAEMuvzyy6cB+KthGKKHEEQikWGyLE9kjP3FMAwEg8H32voMBvCBIjzt6K8G8e6Am3MOxhgsy5rKGKvSNG1dZ2cnU1W1YLsf//jHeRf8OEKIJxKJbMxkMpYoQHVdH0sIqeScb0+n0/Z7ENhALpdHph+m5yC477rrLhBCoCjKdNu2j2qatkvEAK9pGq699lqYpjklm82e4JzvX7lypVA7B9jTbNvu1jRtm6OPvhcvJQxwxr2uCF8dSkhE3w/us0T79u2TTNO8SJKkAy6X66SovpvNZmXTNGdLknQ0Ho+fFDlwAoDX61Wi0egMt9t9PBqNHn2Pr00WwI8E+KqQu3HTT73obTMFxuPxMs75eErpdsuyspIkVYq0MwzDyxibkEqlmsvKyoR0Sc45b2pqcmez2cmyLO/0+/3d58HadAry8X4YvwOSO3/YMwxjGGNsSjqd3hwIBIQWiVKKxsbGgQBmt7W1bYnFYqf8fr/Q7xJCfAAuOHDgwKaGhoZukcNqD2AwQV5bNG9JU1OTzTm3S5y+nvwFnTLt7e12j3FvAlAwLXM4HEYp4zl+/HhJ42eM2b3GX5AURSmpf875m/jf0bwljlt7SH19/ZFIJDKIUjogHA5v7wF6tyRJlYFAoPV07WOxWK1lWR5JkjzhcHjnaV6af4ktyevuhmFM4JyPDoVCfzxd371jS5544glcddVVeOaZZ2quuuqq1UuXLjVFauIwxmY7eUiEYksYY3WSJG2GYNkQ27YnybK887HHHrv04MGD4xYtWvSTvpxTnHPZtu06WZY35x1UxSLzbNseLsvyQREVz7ZtF6V0NCFkpyAEPIyxckppqwimstmsR5blgYSQFsGXp5oxRimlx06cONG8YsWKd64mTldXF5FleYhpmm5CiKuhoWF7Ke0lSRoly3IyEAgU9RzmQe1YUsZzzsslSdpdrN2aNWswd+5cXHXVVdA0bUI6ne7q6upq0zRti4j0bmxsLKeUbhGR3rfffjv1+XxSNBoVnodYLCbV19dvR+46lj8QCPTJGwwGFY/HQ/1+v7CnVdf1blVVhbyYt912m3vAgAEZZ26KkmmaHgDVfr9fKER5xYoV3q6urkGapgndco9GozW2bUu6rh/lnGPFihXvjFry2GOPIRqNckLIxGw2297Q0LBTJDnOz372s/ykz2aMnQSwy8llV5Aikcgb0jydTh+hlO4SkY5z587Ng3R2JpM5csUVVxwu8flJCS8rOcOQV4rcTfZialjPj1OQy/VXCBylDYJSlPgApJT5cUJkz2iCSml2xuDO6z379u2jpml+lFJ6IBKJdGiahoaGhqLtHn30UZpIJC72+Xw7OOcdnHMUklZ5mj9/vhyLxabIsry7sbGxS2SSOOfcNE23aZpjU6nUvkQi0YlzjJzHsND3Jd2+qBxnPxH8+/tA6agF5alUasJDDz30om3bKRFTn9Ouetq0aeOz2ew/lixZ0i361sdisQHDhg2bRAjZ4ff7RWOH4fF4Bti2PVmW5QORSKStVEn2TlB7e/tbSeTeH1tyNsCdl7yxWGy4oihDFUXZdcMNN5TSfoSiKLVer3d3OBwWNWOxVCo13Lbt2mQyuS0YDBY9aS9ZsiR/GBlICKnNZDIb6+vrrbz14Bwj99133/2fZ9h2M85iOoT3Jbhvuukm+Xvf+x7RdR2maU7IZDJKMBjc6XK5mOj1K13Xx7hcLiUYDO5MJpOlmMpG2LbtjkQiO7q6uoq2i8ViWLZsGUKh0DBK6YCTJ09uO3r0KDuH1+GGqVOnHjzDtt3I3WDvpzMF96hRo0YcPny4wjTNiQBOxWKxvY5loCjQOOc8m82OcrlcKVVV94ioL/kc25qmDaeUZiZPnrz78ccfx5133lm0nXMLfqLH48kC2OV2u6VCSXzeZVIAdF533XWv9MPxXQK3LMseRVHGptPpw6qqFnVt59OdeTweL4CZhJBWkRjsvNrT0dHhTiQSl3o8nlbGWOuBAwdw9dVXFx2nz+eTE4nEJzjnh1VVPfYeWIMsgEc7OzvP1Cw7HcCsfiifIbij0SghhEzlnO+PRqNFLQ2BQAD33XcfFi9e7E6lUjPLyso2ybKcEgG2rusIh8OeysrKGYFA4JVAIJAVHaemad62trZZgUDgeU3T3i2LiPTggw+WHzp0yON8LkPfHkfaw1pypm50D3JpifvpTMBtWRYsy1pnWZaQ3ppIJBAMBodUVlaO55y/2t3dLbRwjiVlsNfrHXv8+PF/iD4E55w1NTUNppSOcblcrztgebfmtLK1tTXxgx/8YDVyCSOfA/CTPub6pnfb6nU+k9BWGI1Gua7rwpIlkUiMyGazcjgc3qJpmlAejP3799uNjY1jMpkMAoHA1hLq2tiEkLHJZHJ3R0fHtoqKiqLjzO8QhBDptttuq66uri54SHViJ8pTqVRVsXlQVTXp9XpPZLPZkdXV1dMvvPDCq15//fVr/v73v7sffvhhFwBUVFTwpqamsclksvyee+4p37VrFwghFbquVwmZjhhzAaiIRCKVzc3N+8rLy9mwYcOqeiTq6U2EUlqp63olBMyGtm17OOeljMfr9N8l0n8mk/EBEO6fUlpFCJGi0Wh648aNHWvWrOFnDdylkGmaHyWEbAmHw8d7xH8Umxw+cuTIfwewLRwO71+8eDFEgqa+8Y1vwLbtiymlJBKJ7Jk/f75QSC3nHA8++GDd7bffPqS8vLysWGwJIcS2bdvjcrnKUSS25IEHHvBRSj/udrv/+Ktf/eruF154wfOf//mfP1qzZs1At9sNADh16hTt7u6+7uKLL37k4MGDVV6v17Zt2yvLcrngTuWyLMurKEr5tGnTRLYoks1myz0eT7lIbEk2m/VIkuShlIqOx2PbdoXH4zklGLtSJsuyT7R/xlgZY0xRFCU5ZMiQjrMquQV+HJRSVyQSmQpgs9/vP5GP5ShmSUmlUi5FUeoIIRv9fv9hkcsM0WgUnHOXoijjOOcH8vrq6tWrC7aLx+PIZrOyJEkTP/WpT+07ceLEgVgsdgACwVCNjY3D/X7/AUHrx4RJkyY9/tGPfvR4Ab7gc889l37uuefyc1Hd0NAgZA5saGiQfT5fpaqq+QJLRccfiURofX29UKDSrbfe6qqpqXFpmnZQUKB5CCFd9fX1orElnq6urpRo/9FoNGVZFg0Gg609d953RFeTZdkHYJTL5WpWVfUEgKLAdtp5XS7XZMZYC+f8eDETYd4C4/P5XLIsT+3q6tpHCGmFYJlmxpgiy/L0dDq9Z/To0d0QyLPXU/oJ5i0ZTQgZMGbMmGIL19s5UFLsSo+PswFcVEwYlBLLIcsyeZtjS0rt/03jF73gfcbgzgNN1/VqRVHGMcb2NDQ0pIupIAsXLgQALFq0qDyVSk1saWnZxDkXulxw3333ob6+3pdOp6d1dnZuME0zKTKpTiSZjxAyS1XVdbFY7G1Lvggn4+qwYcP6yqr6ceSKKJ0160z/gfIsglvTNNx3332IRCLjZFkeo6rqRjjB6sUu4t55553QdX1kbW3tNF3XN/zv//6vkAVm8eLFME1zSGVl5TRVVV9PJBJC7VhOWRucSqXGMcZe6+rqwttcOuRiAB2VlZV9vUCXAzh4Fn/PgviFi35wF6ORI0eSSCQyRlEURVXVdeFwWHiLicViEwCQhoaGv4qExvZoNzmZTBJVVf9Wwk0bEELGA+DBYHATAFZWVvZ2lw55wOv1fnjQoEGnc4n/PwDP4OxeDVuHXCm8fnor4M6rPUeOHLmAMXZCVdXtokDLZrNWNBqd1dnZeVTTtP2xWKxgaGwvHXM8Y6zFMIzDfr8fxdJC9AjummRZ1knOeaumaXiHSoccXb58+R6n2lhv+gOAV/thd46Bm3POJUlyG4Yx07bt7dFo9JSmaUWB5ph+3JTSmZIkNTc1NZ1auHChENAMw/DGYrHJnPMjmqad0jQNIvZvRVF80Wh0RiaTOexyuU7Ytv2OZpnat29fXztZF0q4ayhIbuevn0oFdz7tLyHEpyjK7C9+8YsbPv3pTxdNS5a/02eaZoUsyzNkWd6uqmpnXu/u0z7m6OyNjY0VsixPSKfTOxhjXSJ5SwAgEAiUEULGpNPp7bFYrEM0LcTbTNVw0g+/DTQD/bElZwbu+vp6NDU1lRFCqtPp9N9mzZqFyy67TEiFSSQSYznnVZzz10RuO8+ZMweGYcAwjNGMsZqTJ09uNgxD+LAUiUSGe73eUZdeeukW0zTPlUxNBLliRt63sf/+9A5nKrnr6+u7OOcHWQH/bm9JH4vFBtu2LQWDwQOCp3nyqU99imiaVmNZlruhoWGPqEUEANF1fTCltCwcDm996aWXSpkDUWDwJUuWCPFalsX5P6+jD0KuCOr2s7VovW66dzt/fZKT2kH4BbBtm5fC7zjRSuFHqS9kz/GIpnYo6qGsr68vST93JP3kdDrdFgqFdpWQLpgdO3ZsZllZ2eGlS5fuKMac91I5O8IMWZZfVFV1R34x+6KFCxe+oRZJkiQvWLBgUEVFRcGXiDFmc86rIpHIIE3Tii4KpVSSZbk6Go0O+Otf/1rm8Xienj59em2h/mVZrtY0rVZwrlyc8wG6rtdkMplWWZY5ciX7Cp1DasLhsNCdU0mSPISQAaLjyWazHkJItaZpWQj4HZLJpA/AQE3TRJMnDXS73bKDpdZ3Jfn8LbfcIo0cOXJaNpvdGw6HO2KxGEKhEAQys0q2bV+QTqe3qap6ytktCjbQdR2qqkqc89mEkB2qqraapln0gHvnnXdizpw50mc/+9kJ3/rWtyRFUWyPx1Nsh7AAcI/HY4tIHMYYtyyLl5eX25/85CfbOOewLIsUkGS2bdvM6/VagrEZFiHEppRaHo+Hcs6L5i3JZrO2z+ezBWNLsgCY87xCwp5zbrlcLqH+GWMWY4x5vV5bZINwhAtXFCWbTqepoCbw1sGtaRoYY+Cce4cNGza5oaHhjexHhSwieckbCoW8lNJpbrf7T6FQyOq1W/wLJRIJBAIB+P1+b1lZ2TRCyBbGWA2AgsDOA1/XdZ8kSVNuvPHGf9x0002ZpqamNojlLTlFKW0TyVtSX18vVVZWti9evJgil2K4KMVisQ5VVYV4Q6GQ4na7O5z85nldvs+rZk4sTqXopepFixa5q6urOwKBQJsIv2maXgAQ7X/58uWpZDKp+P1+of6j0Sh18pa0l4LNt+y21XUdlNJKWZaHh0KhzfmDoUi7cDg8yO12j2aMbchkMkJWjUAggMbGxhE+n28MY+wfnHNbJExBVVUYhjFaUZTRnPPXz0TvE6Xrr7+eP/744/8O4D/egXPTNORu4/TT2QL3Qw89lJek02RZrrBt+42ybiIewEQiMV1RFG8wGNxWyiEqHo9PTqfT3nA4vDUYDAqB07ZtFovFRkuSVN7Q0LA1FAq9re7qyy67zLNr166LAfz8HbLG9NPZBPf27dtJNBqdmclkTjoWESEKhULUNM3JmUymKxQK7RNsxru7u2k0Gh2XzWbTmqbtLOGgSsrKysZalkUCgcDmEtpR5BLelPdhpeiT0um0b+bMmT/HWa5l3gd14V8jDPvpTHXu5cuXy6lUajohZEcwGOy2bZvIcvGuvva1r0lut3sG53y3pmnCQeeccxnADMbYgUgkcuLnP/85rrjiioJtHn74YVx33XXgnNfZtv1XXddbSqgsTADc6VggKgD8tqKi4s+dnWLXMhcvXtzmdrt3vvLKO3KhfSv67dxvXXJzznljY2NVKpWaqChKczAYFDLlEEK4ZVnVdXV1kwA0B4NBYWBHIhGfLMtjbdtujkQiJwAUBHZeMm/YsKFc1/WLOOctmUzmRF7PF6Q7katscC2A6wAsfvbZZycJpjBWqqqqeB+xJW8H9QP7rYKbc46ysrJaADNDodDWZDJZNDF8PuovnU7ffEkFAAAgAElEQVR7CSEzVFXdcurUqaLt8qbDefPmyYqifNCyrG1VVVVCFcycgyqtrq6+QNO0f1BKO0sMjJ8C4NsAVjmfTwKwN23a9CHB+fxmfX39OxmCWobc7Z9+KhXc+csFlNIaxtggv9//gogUzEf9GYYxRFGU0ZTSF4B/ejALHDQRDAYRj8eHTZw4caSqqs9RSmURfdcB9giXyzXc5XK9eIZz8p8A2pArgffG5tPR0VH2+uuvF2t7K4A/NzY2vpOXB6YAuKAfymcA7jvvvBOLFi1yE0LcXV1dm0Q6bW5uRigUQiwWm8g5dxNCtogKz0AggHg8PoEx5vL7/XsED3Lc2SnGybLsDgaDB9Lp9Jlu1/+GXF3HpHMmkQEMcLvdmQsuKIihcuQi/ta1tbW9kxaMfmvJW5Hcd9xxR5ox1tKz/nYh1WD8+PEwTfPCrq6uVsciIpzB1DCMC23bbg0Gg3tFY8U55zSRSFzU1dXVHg6HmwEIheL2QTMATEauau9Lzt/AwYMHHy3ygnYCuPtdWMM2R3U6e0q8gMfz3aJSqj8XNXEUy83XmyRJkuLx+EWZTGZzY2NjV89SH4Xm89SpU0osFpsIoDkUCnWIxop3d3e7XS7XcEmS1pum2T5nzpyidnbTNNHQ0IDy8vKqJUuWzHC73ZwQgo0bN3offfTR0SNHjvR/5StfeRYA+clPfvKBffv2fWDmzJmpK6+8cubpFr2ystLq6OiQe75snPO6eDxuE0JEVBQbwKREIiFkvbJtWwYwKR6P25IkgVIKy7JmFmrDGBuWSCQqIJa3RJEkaUwikRAVSm7GWFkikaiBWN4StyzLtYlEolKw/2pZlmEYRm1ZWVkzgNRZAXcplpR4PF7W3d1dk0wmX6usrLQAYOXKlX3GluSBKElSuc/nq6qsrNx+8uRJu5hOP3/+fKxevRrxeLwMwEjbtg/U19e3+/3+gsDuYQr0PPjgg2O6urpORCKR9T1YhgDI7tu377fLly/PK9hfA7Bh2rRpv1u8eHFfhZWuQC+HTTgczkaj0c2i8xeNRq2edYQK0U033UTGjRuXbWhoEE5dHIlETkQiEaG7m1dffbU0derUTk3ThMYTi8UUSml1XzWPTiMwlY6OjsGiqR10Xa+2LEuKxWLHAWDFihVFE7CeFXBrmgafz8cZY0MIIcmqqqpdIqGheaAFg8EhlmWNJYT8ZcGCBUUDdfKxJaFQaBDnfBjnvFmSpNpiOn0+tiQSidTatj38pptuWo9/3hzPWzcs5Bwi3T0sEVcAuEtRFLJs2TJyGnNgXe95rK+vlwghJc0tIURYdRs+fLjcQ9Ub5EjLowVenLyvQLj/UsZPKZVLwZJlWSX1TwiRKKVvzI8IsEsyBfZ1+LvttttIJpMZa9s29fv9O0WB3drailgsNs7tdnsIIeslSWIi7QKBAEzTHOvz+SpUVV1PKeXFtkJCCFFVFfF4fJyiKJXhcHh9H6zHAewGMML5/EVHZfhhd3d3X4sxB8Cz76IaOgbA2P7j41kGd3l5OVm1atUE5ErZHRJtd+zYMTJkyJBxhJBkKBTaK3LgpJRCluV8JbLOQCDQLPp7mUwGmqZNSaVSJ1VVbS4SgmsAqAegA/imI7kLbbcrnUPdu0n9jpyzBe6tW7cikUhQxtjknTt37mlubt5fiklqyJAhF6TT6b3BYFAo/RYhhDv29g9ls9k90Wj0WCmXICoqKj5MKW2ORqMnNE0rlltlLXJeyR8A+DyKp014t+M6DgLY1w/lswTue++912vb9kXbtm3b/sgjj1ijRo0SATYzTdMXjUYvTqVSm2bNmlVUDckXZuKclycSidkA1ucD6EVc6ZFIpNIx672mKEqmj3bsNFKwFcB+FL6+9VH0CKp6F+kQgCP9UH7r4OZNTU3DBg0aNC4YDP5t6tSpTNQ5wxgbbVnWuHA4/PcLLrggU6xKwqpVqxAOh6Hr+hBJkgZ3d3evk2W5W7RYk2mag2VZHgdgVyqV6u4jR4rXUTsWAmgAsNjRoccVOSB5kPNkdvZD6DwAN+fcdrlc4zs7OyW3270lkUgIJ8cBMDabzVrhcHiTpmlFk2QuXLgQt912G3RdH+PxeIht27sHDhxYtE5lvoJZLBYba1mWLxgMrkPhWzZe5NzsPwFwP4BHAax3pP1CAAsADDxNuy84vOcCjcTblzbi/Ad3JBIhsixPAdAWi8UOHDlyhIsURH3hhRdIPB6fRgjJeDyeQ4QQIXWirKyMRKPRSR6Ph+TT4gpcvIfX6yWRSGQ8IcSladoegUc7gVypuxZHFTkEYBuA3wBYDuD/c0D+HwBQU1OTb/cqzp3yeEOdv37qRUK2RrfbzW3bPpjJZNoBYPny5UV15XQ6Lfl8vqmc88Oc81Oc86LeKNM0oSgKSafT0wEcrq+vL6lgk6IoMymlB1RVPWEYRtGknAK0BYAG4LMrVqy4zrKs1xyd/Fw7wPVbS85UcgcCAdi23SX8xshypcvlugjALtGKYk5pEZ9lWZdyzneGQqFSgK2YpnlZd3f3DlVVT5imKQpsjuKXgxmAp0eNGrWxvb39VuQuLxSkpqYmWyQJUW+1T5T32LFjNv6Zlq0ZOdt8n+TkLRHuv7W1taTxc84tlJAmTpKkUvu3ndItAADRkjKleNGETo5+v9/HGBunaVpJCR87Ozt9Xq93ZENDw8si/I2NjfD7/bBtu4wQMjYYDP4m/3+FzgJ5DydjzF1RUTFgyZIll7hcLi4AvExzc/PnL7vssj0f//jH9xd7VznnE+LxuAIBGz7n3JYkaWI8Hi8XnGeZEDLBNE0XpRROOoWxxdSXeDw+GGL5zGVCyGhnPCLgc1NKy+Px+DGR/lOplEuW5Zp4PD5ECHiEVLndbsTj8ROU0i1+vz99tsFddIFM0xxq2/Yg5yAnXt6B0sFut1tRVfV10XZ+vx+NjY0jGGMj4GRzKhY0NWfOHAQCASQSiRE/+clP5K6urvYnn3zytS9+8YusiLrEr7jiikkdHR1rf/Ob3/z84osvJkXmQlIUJaPr+jpVVYsutmma3DCMTGdn5xaBwEtIkuRyuVzpYDD4OuccO3bsID/84Q/7LKxFKSWyLI/WNG2PyHiy2aynvLy889SpU5tdruJVAL1er5cxVq2qakux/hljqKqq8mUymUGpVGqviLXN7XbXMsakZDJ5pKmpSVgFe8vgdnRbbhjGeM55Wzgc3hCNRvNmvKISNBqNTrYsS/J4PBtF2nHOQQjJt2sHsJlSOgQofut+7dq1iMViE1OpVOc111xzEAA2bNjAN2zYUHTCLMuyPR7Po8lkEqZpFuSvr6/niqLwPHAFhQNvbGwU4nVu73NHqo1xVKf9Bc5Ab6RHExnPokWLGAC+fPlyofE4fQr3v3z5cgaAx+Nxof6j0SjjnJNSgF2SKfA0Cwggl5XVMIzLCCFtwWCwpb6+vmA6s7x50JGgF3DOj0qS1GLbNhGxYV922WUwTbOOUnoyGAweoZQKqUvhcJgsW7ZsFCGkW9f1w6U+78c+9rE9jY2N3Tj3aBBy0Yz9dLYkd1NTE2699Vb3oEGDRnV2dq73+Xzj898XIqeQqpdzPty27W2apiVjsdigYr/HGINhGD7G2GjO+T7Ry8kA4HK5vISQiel0ujkUCnWJpF3rfWDlnCOdPiczKJRUbKnfWiJA8Xi8avDgwXWpVGp/dXV1t+gEO4U+J/h8voPBYDApqtJnMpkKxth4t9u9WwTY+UsSnHMvgAnpdHprKBTqKnbgPJ3QBvBhSZLOVXPbLrz5vmc/vRVwR6PREZzzoalUapNpmunubjEhahjGYFmWR+i6vtHv96eKJcjs4eGsATAqHA5vrK+vT4vU0lm5ciUaGhqqnerC6w3DyJZyRSl/VgJwNYA/2bZ9rkrHE85fP70VcHPOWSQSGSFJkktV1W2GYRSVZrFYDJRSZlnWYEVRylRV3ZL/v0K26MbGxvxt9oGUUo+qqptvvPFGAChYSyf/wkQikRGKolQA2JK/+3kGpUMGA/ge+p0k5z24qcfjGS3LclpV1d3FmPNSsry8nDLGphBCUvX19UXb3XXXXQCAQ4cO0cbGxgu8Xm+WMbZ/5cqVePDBB4sOsrm5mRiG8SGPx5PRdT1vQThTcO4FsOkcX8OJyKV36KdSwa1pGlasWAFZlmcyxo4Fg8FjAvm2oes6Lr/8cnR3d19s2/ZeSZKEskwtWLAATz/9NBk6dOiFR48e3RoIBE5B8Pb8/fffL40dO/aCvXv3vu73+1vfJ2tYjXMj9Pa9B25d13Hq1CnJsqxN2Ww2WUydyFMkEqn85Cc/OZMQsk6SJOGKvaZpVvztb3+brijKxjvuuEOoHeechUKhipaWliler3fr/fff/1YrBJcjd8XsvDcMnM8kZArUdd3Wdd0StYgsW7ZsYDabHaSq6oZYLAZCiEgpOd7U1FSbTqdrKioqNrW3twvXqqGU1paVlbkJIVs7OzuLxizcddddWLBgAQCQuXPnKjNmzHijTSgUsiil8yorK38TCoXk7u5ubN++3QYgpdNpORqNijh8JEKIpOs6Ebn4m06nbQByNBoVTe2gAJD9fr/0s5/9bIfH47Gvv/56uaurz/AfQgiRo9GoBDH3uAJAikQicl9ez17CRSaEyJFIRBLxO+THL/q8ABRCiBSPx+W2tjarWOBeSeAuhRobGy9Mp9NHy8vLtztgKVo2hDHGZFm+xLKsE11dXdvD4bDw7XkAkxljg9va2v62YsWKou3i8TgWLFiARx99dNT8+fOHT5w4cSQh5A33+5VXXjlYUZRBfr8/nUwmRxFCMGnSJJsxNgjAKKe4UeGVUBSJEDLYsqzRVAAdHo/HppQOYYydEgEfpdQFYEhZWdmo6667jgBAd3c3CrmyKaUjGGOWiJT3eDweSZKGEEJEL2O4KaUDKaVUJE9LJpPxSZJULVJAzBn7AEqpy7ZtNyFkt+gZ6qyA2/HsSvF4fGwmkzkVDocPipjdYrEYJEmilNJRjLGTqqpuL6VdLBYbyTm3CCEbRIDd1NSElpYW6lRY6D5+/PjeeDy+qxfbbgCv9I5Xb2xsrGloaBC6lHzllVeSGTNmeKLR6F7ROYzFYr5QKLRHhLehoUHy+XxuUX5HTbQikch+Ed5vf/vbytChQ7mmabtF+E3TdDPG2sLhsNB1tzvvvNPt5C0RGk80Gj1pWRaNRqPHegm2t19X83q9lBAyTlGU4+FweGcJZjeSyWRmWJbVDccRUahdHvihUAi2bc/IZrMdhJBdEMxyWl9fj6FDh87IZrNtc+bMOeq83D3noJDUpMuWLRNSy+rq6qhgpikI/vabSJIk2oN/BoDZRcCBUsbjrGdJlrRSsGRZVqn9k578oibdMwZ3HmiRSMTHOZ8G4MCSJUtOFYu1nT9/PgBg0aJFLkrpTI/Hs51zLuSE0HUdsVhMTiQSFwDYEY1G20VBUV9fL8disYvdbvd2XddPl1tPAvDl9+C5yYP+FMZnD9yBQAC6rsMwjGpZli9QVXUD5zwF5EJRC9Hq1avh9/sra2trP6yq6nq/358SAejixYsRCoXKOecf+cMf/rA+FAoJu+4ppRVVVVUXEkLWP/HEE6k1a9acjm8q3psxGvb7HMMDkTOH9qQBADxnpHM7l4OHMcYGh8Phl8LhsDAoNE0badu2t6Gh4bkSfpIlEomxqVSKhcPh50taedsezjnPDh069NW8h/M0FYYJgBsBhN+Di7vhfQ7wLyOXGOkLAP7sfPcAgOMlSe78aTwWi00lhCAcDq/L5xYR0LMswzDqKKWWYRg7SgleopSOymQyGV3X94kcOPM8sVhsFCGESpJ0YNOmoo7GRrw3UzVYeH+HB9yP3KXuj/T4bj2ALaWAm3POqWEY023bPmoYxiFN0yCYR4QSQqbLstymadphTdMgkns7kUhITlrjlKZpwgWbXC6XZBjG5Gw2a0mSdMCyLLJy5co+n8v5O/oeXVwJ/U6cvyMXvUkcNeUKAD8WVksYY4rb7f4A5/y5SCSSWbhwYVGgaZoGj8fjtm37Qs75fr/ff0wEoK+88gpeeOEFVyqVmsQ538U5L1+5ciWK5frmnONb3/qWYtv2xHQ6vd80zc5EIjGgUBMAlzhbe+o9urDTHYC/9j4G918BfMaxfi0CcAeAI0JvfFNTk8fj8QyzLGtdMBjMAGJJ6SVJqs1kMiMsy9ooy3JXMTNOvvLw7373uwFdXV1jPB7P1nA4nBSpJ+6YvKpGjx49OZPJbDNNU0TFqAZwDYDMe3hhXSihcsV5Sq8hFxb9IeSCyH4kZC1pampCfX19Kp1O77Ntu6h0W7ZsGQDAMIwKSulwTdOaZVlOi1givvzlL2Pp0qU+zvmYcDi8Y/HixZbo0+m6Xk4pHaOq6kav1ytaTWw8gPsAuJHLxd3Xn/vQoUOeSCTiK8JXBqDshz/8oe/QoUMe5Mx0ZQJ/7gMHDngBCPW/du3aspaWFq/Tf8+/0/JHIpGyXbt2eZGLTy/a/29/+1vfgQMHRMde9utf/9q3efNm0f59v/jFL3w7d+4Uft7XX3/dt2nTpjy/5zTruNPZeR8EsCz/ZVG1JH9XUrTc3dKlSxGLxcbZtm1rmraecy6kX1NK7R07dtQNHDgw7ff7X29ubsb48eNFdgfLtu06WZYzY8aMWR+NRkVv2nwJuQq/fuTya7v6OJiRmpqa6jVr1vi6u7unCBzeeGdnp++pp56qclSGovM2ZMiQwU888YQNwapkhw4dqn7iiSfcyNV9l50xTepLoxw8ePDwJ598sgtAsTKJBAA/fvz4kCeeeCIrOB579+7dI1999dXjAjsgAcB27949+vXXX28TOMATAPaLL744vK2t7YSzXocB/LoX3ylH797oqChi4C6VTNOcIknSSb/ff8iyLBBCisaWUEqJbdszOedH/H7/4cbGRiFgh8Nhks1mpxFCjqmqut80TZEDrgTgvx1AfAzATQA6ADx1OubNmzf7XnrppWHz5s3bJfL8zzzzTOXBgweHzJs3b4cI/+9+97vKvXv3Dpw3b94eEf6nnnpqQGtr60DR8axdu7a2ra2t7OabbxYKBfjhD384GID761//upBr/OGHHx5BKcW1114rVALkoYceGkYI4ddff73QJe0HHnhgtMvlylx33XWF+H0O8BNvwtVbBXM++SSllESj0Rkvv/zyDr/ffwgACpXMzkvzSy+9lACYQQjZGg6HDwOFHUGxWAwAcO2111KXy3UxIWQvIaQNELobSQAEkLuE8KMeNtEXTsccDAaH/vSnP508b968XYVMkPfff39+bGNfe+21wfPmzduRSCREBMGkv/zlL6PnzZu3p1D/+bmKxWLjNm/e3BvY9HS7g6qqiMViF23ZsqXi5ptv3itSeMswjLoDBw6UFQJ2fiyf+cxnEI/HL9u/f3/62muvPVjsYrjTdvrBgwflQsDuGfFnGMalO3fubO8D2O4ez30bcolJ3+TpfsuSe/ny5ViyZInH5XLVZTKZn/3yl7+0RSqKOWU8Km3bHgFgB4qnNQOAfH3LgZzzIbZtb1YUhQKoKtSmrq6OOW+37ujYvaXq8TfQT8gwzvmhxsbG8dlslgSDwdd+9KMf4frrr++z/2984xuIRqOXAtjZ0NBw3DCMPjPSRqNRjBo1iuzfv3+CbdvJcDi8vVBew0996lP5ufpwZ2fnpng8fjKfbcuhWQ7A/5E/I2UyGYVzPpZS2qKq6uE5c+agL1PoQw89hBtuuAGmaV7GGFsfCoW6YrEYQqFQn+u2bNkyT3d39yjG2LZgMHhC07Q31NfTPa/L5XJbljXr5MmTG5ctW5aMx+N9XhVcsmQJ7r77brmtre0CAJsSiUTnwoULT2fA+BpyaS0qHPXm8d4MZwzuZ599Fpdffjni8fhIy7J8nPOdhmHYzlYoIiWGM8aqjxw5smXIkCG1IpmWnMkaTCkdvmvXrnUPPvggj8fjBXP3+Xw+dvvtt89A7qJvAkUu03LOvzZ58uTUwYMHH1m9enUrgILANgxDymaz4zjn+8Ph8HFH4vfJ73a76b59+2a1tbVtXbVqVaoY/2WXXSZ98pOfnGpZ1uZ4PH7yNDvbmyR3JpNRJEmatG/fvh333ntvpth6rFu3juq6Pi2TyfxN1/VsXoAUOLj7urq6Ru3atWvHI488wopZwI4fPy7V1tZOCoVCb6TXK3QH1jTNyvb29pru7u7Xy8rKrAKWuZ0APoCcV/Kh0zGcEbhXrlyJyy+/HNFotI4xZluWtUOW5Zmi7TVNG5/NZpXy8vLNd9999xuqRiFavHgxDMMYTgipdLvdr4vcp3Qk2SUej+djsVjsB5IkjXcsJIVo5y9+8Yvxzz777Je/+93vvtLZ2SkVMcON9ng8BzjnqXg8PrzYudmyrAmyLO8ePHjwjGK3+CmlPtu2hxFC9lJKJ/bmlySJ//rXv56USqXwla98hSeTSTeldBwhZNfo0aNnFemfM8aqKaU1yCXTnF2MnxBSTggZYdv2zmnTpl1UjJ9SWssYqwTQHI/HLym2VpzzakJIDYC/JBKJYpayZ1Gk0NYZgbuyspKYplkHoF1V1aOlxJZEIpFJkiR1hEKhUrI+kUQiUWfbdkpV1a0rVqwQbfe1lpaWL0iSdE8oFHqVc041TaOW9eZ5O3HiBBszZsxAxti4L37xi0994xvfsIYNGzZn6dKl5S+99NKLF154IQ8EAlL+3kE2m+UDBw4cZFnWIM75z7u6urKSJPU5By0tLXZdXd0wSukI27af7O7uLliR4siRI2zUqFHDnWy5r7S1tdl98XPON4wYMcJOJpMjCSEDOjs71zDGUOiOxGOPPWbfeOON4yilHclk8k+2bRe0AFmWxauqqsbZtt115MiRP1dWVhY8q8Xjcds0zQm2bZ9Ip9P7e893b8pkMqisrJwgy/LxlpaWP0ycOPGsBLCVDO5ly5Yphw8fHq8oSrvf7xd2WXPO5VgsNgXAiVAoVJKrmxAynjH2XDAYPAYI1SF0AbgFuTzac2zb/g6A/yCEPOEcJk85p+u8HjnswIEDNd/73vf2qar6ceS8XS9alvXniy66aIhjZvtzj615IOd8WCgUek3TtD6zbK1ZswZz585FLBYbzTkngUDg1aampoLbvnNwHEEpLVdVdbOmaejrcGqaJv70pz/BMIyxnHM7GAy+LpJNKxqN1lFKM6qqNmua1ufOGY1GMXToUNLS0jIWQEePdWMF1ApiGMZMAIdUVW3VNA2FJHw0GqU+n2+k2+3uXrJkyX6cRSqlbAh8Pl9ZJpP5MOe82e/3HxGxCABAMpmUOOcflSRpeygUOiqihvSYrH9jjDXPnDnzmGAt+CrHkP8HAD91vrvH0bezyFUqW4VcPpJViqLcd9ddd915zz33fINz/gVHlwvhnxWBjwIYDuB2R6UawTmXZ8+e/drjjz/ep745f/58zJ07F9FodBqAU8FgcG9Pv0FfVqBIJDKEEGKrqrr1pZdeKqjPqqqKhoaGD5w6depQKBTabxhGn8COxWJYtGgRMQxjMKX0uKqqeznnBfsPh8NobW29KJvN7lVV9Wih+Y/H41i1apXc3d09gnO+Kw/svvo3DCOf1HRme3v7gSVLluzPOwDfMXDnTVSSJJUxxiYHg8Hnw+FwFkDBGjU9kuNUlpWVTbcs6/mGhgZW7MCSX2TDMKoNw5ikquofKaV8x44dfS5cj3rsU5Er2pRALjKsJ3U6Rv67AHzn0ksvXXzllVfe9+lPfzre2dl5NXL1bx4EsAdvjjPhAH5CCPnjjTfe+HG3290diUSO/v73v0eholVer9cdCoXGEEKOhkKhtmLza1mWYhjGbMuyuoLB4GHDMPDhD3+4kAnWE4vFJv74xz+uXLZs2ay+Dqb5/gkhvurq6hnt7e0dwWDwhPNdIStYua7r02zb3jJgwADWl6k13z/nvLqrq2uyJElH82nrTgfsPC5SqVRZNBqdLcvy5hUrVrC8A/BsUlG1xElc6ZEkaWhXV9crIp3mTVuGYYywbbucc/6ay+UaIGrqi0ajIxRFKa+vr9/WC7yn10FcLoZcXMGnAUSRCwMtpOYon/vc56a2tLTsv++++4ijxvTpXbvjjjuU7u5u68iRI680NDSkAbjvvPPODPrwVpqm6eaczwCwpVheQ13XEYlEFEVRpm3YsGHT448/bhWzoMTj8SpCyIQpU6as37dvX3Wx/nVd9wKYFA6H10EgPDYWiw1Jp9ODL7vssg1lZWX4yEc+Uqz/agCjVVVdX6xvBxdDCCG1qqquO4MUd2cP3M7WktJ1fZckYK+7+eabEQwGEYvFptm2faK7u3t7RUWFS3RA8Xh8ajqdbqusrNwuoj8GAoHk6NGjr6GUHps1a1bjb3/7W/v555+npwE05syZwwDMp5R+6J577mlobW29m1J6hDE2w+12f7m1tbX9N7/5Dcm/UHPnzmXhcLiqs7NzVFVV1ZZLLrkku3r1akiS9B+Off5/H3vsMZqfloEDB/Lnn39+MGOshlK6/qMf/ai1Zs2aPndHp/8BiqKMqqio2PTVr37VvvLKK2lfu9Po0aP5s88+OyqbzcoVFRXrr7jiCsvlciU551I2m0Xv35o7dy7TNG2YJEkDy8vL1z/11FNIp9N99v+Vr3yFRSKRCdlsNjVixIjNHR0dpL29nZzu5hJjDFdffTUzDGMsY0yura3d8OSTT1Lb7vvexNy5c5lpmuMIIempU6duyudof7tI+FSq6/pY27bbotHoydPoZkSW5ZnhcHh9IBAgFRUVF2YymUORSOSQpmlQFMUtSVJlIBBo7UNS1BJCTgEYwzlPhkKh/fPnz8fq1avzO8GFXq/3tdN42RQAC2tqao4vW7bsz0ePHnX3lXqBMYbvf//73zp48OAXrr/++rnPPvvsFzo7O4c/8MAD8Tlz5rw4dOjQHy1evPiRTCZDnJeBWZY1nHOelGW5458Y4KSiosJqbIFFbVMAABmFSURBVGy8mDE287bbbns4n6GAcz6Ac14lSdJOznkxwcFs2x5BCOmilJ50VMSCUpUxNgSALEnSPs655JgL4Xa7WXd3t9Tr2TljbAwh5KTjwaVFdjNYljWWENIlSdIRznkxlVWyLGs0gHZFUY7lx1PozAZgCiGkxbbt1tbW1t2rVq2y8DbSWY0tCQQCks/nG8kY2xKJRJL5KggC6dcUy7LGeL3effX19akf//jHuPbaa4taJJErjHrf8ePHm2+66aaCzEuXLp3c0tJyNWNs+f/93/+9ilzc751XXnnlHgAnDh486F20aNGmvPdP1/UJlNJ9oVBoex9dbgFAli5dyp0XcHI2m91RLLH95z73OTzzzDNIJBIXZLPZTZFIpKXIToaGhgZomjZFUZRNoiXFDcO4LJ1O/zkWixUMlkokEshms0SSpGmMsZdFC23FYrF/27t37zM/+MEPrGJqEWNMBjDZ7XY/GwgE2vAO0VuOLdE0DRUVFZxSWuXz+epcLldLPu92kQNL3ltVRSmd4vP5DtTX16cAiAB7MnLxBPc6DggUWcDpzz//fLVt2wPxzziSGgBbe1hYTgCAoiiyYRhTAJwKhULbi+j73PEYfvh73/vebF3XDxc7H1x88cXuSCQyxbbtfcWADQCHDx+WY7HYh23bbskfNHtRhfOi518Gn2EY02zb3lQI2PksBADKJEn6UDab3R0KhY4Vs2SFw+Fy0zRnptPpV37wgx9Y+Rj83pT/XlGUAbIsz3K5XM0NDQ3vGLDfMrjzpp5UKlUDYEAoFNpaX1+fFWm3ZMkSNDY21nq93uG2ba+zLCtbxJyY/+ds5KL6VjmWjUK/oxiGMbOysnLfK6+80uXYZ9uQC7rJIhcNOAFALYDnfvvb30rJZHIaIeSgpmlH+3pB84cg0zQ9pmlekEgkXm5paeEAFhBCyvta6Hg87vZ4PJMURTmQt1j0Zb92zIZlgwYNmh4KhV4yDKOjj4PmZDhZXg3DKGeMTfR6vbsopR2F5t/JQlDBGBurqupfIpFId1+WrJ/+9Kdv2PddLtd4xtgmn8/XDZzeta9pGtauXQtVVStt2x4eCoX+oapqEu8wnTG4Gxoa8nlExnDOaxhje3pJhIIvRDQaHZlKpcpDodAWkXx6zqJ+FcAnkAuA6igEjkQi4ZUkaZJt29u+853vnHKcN2nkStoNAJAP0VQBvLp8+fKX//a3v00PBoPrVVU9Vcim6+QNr+KcT2CMrXfSsa0F8EvkkuT8y0JHIpEqznldJpPZqKpqZyFVTVVVmKY5sKKiYlwoFFpXzAUBgN1zzz1DAYxMJpMbFy9enC4SH4JIJDKkoqJirKqqm3q+sKczEFx55ZUwDKOOEFKVTqc3hEIh1ld8SH59NU0bUVFRMSYUCm0SwcW7Srqujw2Hw1UAsGXLFjix0xMXLFhQ5niaZhXQ/9yJRGJQD31tuqqq3h6fBzk1G/syD34AwLchmDQnkUiMicVikxsaGnq/vKrjaWwE8CSAOwH85pJLLpkVj8cvFJ2LRCJRF4lExjpmRJG5G5xIJC4S7d8wjKmGYYwQZB/38Y9//BPLly+fVoJFakY8Hh9aAv8kwzAmivJHo9FJ8Xh8FCBeEPWcsZYoikLC4fAlbW1t6+64445MT2tJX+CWJKmyubn5+JgxYy4khGyrra3t/OY3v/kGuCVJam9oaDidakJrampWDx069PnFixf/qaWlRSl04rdte6AkSXWc8xd7H5gHDRpk/eEPf6j5xS9+YXi93j0zZ85c++UvfznT3d3tsyxru8ABmwAYTggZzhh7pa+UYLIs8yeeeOKC9evXj/H7/esVRSGEkB2C/U8ghCic883F1ofkaJQsyyydTu8ttgM654FpjvWkRWD9JUrpKEppyrKsQyJ4kaT/v72rD47quu7n3vexu9IKrT4s8SWEhCQExMKYYFxaZ9rx1NMmaT46ydiOx6WeMniIh4ACSFrt7nv7vqU24BpPmbEm45SObcjgNHHjppOOPY7DtDY1JkGGCDDfGPElGa2MpN33cW//2H3Jhmjfe2TcxmCdfzSS7jt7333n3T3nd8/5HabDcZxzw8PDx3bt2mX/ITfkW0VL6He+852qqampeYSQX+zYsaNkQHHzdYSQmsbGxnqE0GAikQjan6YeALaMjY39fX9//6Xh4WHkkUBECCFLKKUZAPi36QaOjIzAF7/4xdF9+/Zh0zR3P/TQQ+jy5cuXKisrR/3K6BzHQSzLLkEInSOEHPKi6nUcB/bu3fvvjz/++GPbt29/bN68eel169Zdt6zSYYVt2xzLsm0MwwyapvmRVyJWAbYjAPBZhNAQpXTCj0zWsqwwy7KLAOAtSmnWa/6FpWAppR2U0kHHcab81ofjOD6bzbZOTU29WVZWZjU1Nf3BiYICGzellHAcVzc1NcVUVlYOZTKZX+dulyrtWr9+PQwMDADDMLWEkMqGhoZ3zp8/75sDXJDPQp5/QnccZ2zdunXT+ncNDQ1w+fJlzLJsu23bF0VRHPWhjwgDQMsjjzxSDQDv6ro+tnbtWti9e/e0Ef++ffvAMIwQIaSZEHI6mUxmAnyNh1588cXlqqr+4MEHH/zpiRMnbnR1deU89JdjjBuz2ez7XgiHe1/xeLy8srJyMcMwh7Zt21YG+cLcsVL6dV2vZFl2Ps/zQ11dXVN++tPpdIxl2cbJycl3dF333X1lWa5yHKfBsqyjqqrmApbbfjICyldffRU4jluAEMKpVOr46Oio45fZBgAwMDAAkiQ1IYSqMMan1q5dS4JcB/muBg8CQHK6h+YGjZIkwf79+xEAfO7w4cMnksnkaKkXx81M6+jomBWLxf7zypUrb7u6pzNswzBg37598MQTT7C2ba9MJpO/8jJsl3lLEIQoIaQDAI68/fbbUwBwDQBcg5oDBRqGb3zjG7Bv3z7o7u4udxynPZVKHcUYf+QVJBcCtVhZWdnSnp6eQ9u2bcsV0J7W6a4pbDx1CKEFL7/88tGJiYkpr/kXDLuG5/nWZDJ5+K677rJ9XBxQFKWB47i6wcHBQYzxJ8awAxt3oYDXMk3zUpBdV5ZloJSCYRiLeJ6npmme8PPXbNt2/7/BhcK9xicSCVBVNdbS0rISY/zWnj177PXr13uiO7quNzz88MOzr1+/3vnGG29MeGXQxeNx6O/vj7W1ta1wHOcAQgiee+45L/wXZFluYlm2HgB+IYqiPY3+KACkEEIVL730EqiqWh2Lxdpra2vf9VvXAoLSyrJshGXZd29Cc+h0UKVhGDWRSKQqHo+/98tf/rKkfrctuSzL7TzPk/Hx8XcMw4DOzs6SOzxCCCRJWoAQ4uLx+PGXX3759+kW94d3Sz7/+c/TwslboJeB4zhOVdUOy7JOy7J8XVVVz7Yh4XCYDgwMRCFPRnm8AKl5iqZpcwkhkVQqdbD4m8LDuCsIIbMKUJ+n7kKdZpNpmlZxeZQbAJcIuGsZhgklEoljHqrfB4C9lNK/6+np+YlpmhlZlt8NsqaKoiwwTdOWJOnSTf8ahZtIeaqrq3E6nW42TfO6KIrH/XTv3bsXSZLUSggZC3LQwrIsK0nSHELIlCiKn9jGWh97bkkikUDhcHjFmTNnBp9//nk7SG7J6tWrW44ePRqfO3fuP61bt+6kVzcJx3EoxnihaZoQCoV86QoQQoQQ0mJZlsMwzFmvQK0ghFLalsvlxiORyJUgbFeWZbURQm7wPD8c4GuZEEKWmaZ5ORKJfMjzPORyOeTx9Y9t215KKR3mef7DUu5B0efSXC63PBwOn6KUTvhNxjRNFiG0zDTNkxUVFRNea48QgsnJSRQKhZZns9khWZavBeVvvB3QkpI3TSkliqJUMgyzKJvNHn7++eedjRs3ujzeXpevPHDgwBdWrFjRv3HjxtErV66EPawCA8ASSulFjPGHlNKwjxUhx3GWAMBFnuevU0ojHsfjlBDCAcASjPGpUCg0QQiJ+CAWHKW0hWGYixzHjVFKw176HccJI4TaMcZDoVAod+rUqdgrr7zyt1/96le/19bWlis2coQQNU0zCgANLMueRwjdmO5+eZ6nhBCwbZsSQmZRSpt5nj9a6H8T9pg/JYRUcRxXizE+y3Gc7TiO13pS0zSreJ6fTykdwhhPbdiw4RNr2L/3Ic5NGDZs3boVdF1/SFGU1umA+5sPcYrkTwpBo69omsarqtra19dX6zXO/ex4PM5omnb/1q1bA3UdkGU5Istys6Io5V6nk64fumXLFl5RlJWxWAz5zNv1f6O6rjclEgn+JtrnZZDnl55d7C93dXWVKYoSpOi6FgAqNU2LSpLUXBzcTiduTk8qlarSNK3NDZ791lMQhDLDMJYC/IYy746R6YzbfRCqqq7QNO3PSl1bwrg3AcDfeH1mUWuSiGEYn+3o6EB+D6Ow+OWiKC40DAO7kFgpqLJggDW6ri//yle+EshN6+zsjKTT6fYXXnjht/SUmr+maQ2KorQkk8lS+iMF1MM1pEpZlttv1lNC7l2zZs3jO3bsuGvz5s1+frtrsPWyLC8GADh06JDv/fb19S3u7e2tfO6551DAUr/bBy3x8DUZWZYXWpZ1ybbtkYCXhQGgG/Ipo//itWNIkgSaplWxLLsgHo8fHBwcpIVd2WsHnj81NVV36tSpc4WDjmmTexRFgYGBAZBlud6yrKre3t7DP/rRj6jfi6ZpWlNNTU1tOp0+9v7775cMZPv7+2FwcBAVqpGYVCp1UlXVUvqnAOAUAMBTTz1177lz52oFQTjmfkuU+uqnlKIlS5bMv3DhwpVvf/vb1y5enJ7R7Otf/zpQSmFsbAxrmvYZx3FAEITjsizDvfeWzgpoa2tDkiS1E0LGdV3PHDlyhN5K04Db1udWVZUnhCxDCJ3kef6G4zh1PhE2hXyaaRcA/HPBuEtKT08P9PX1LbRtm0smk74R/+bNm1F1dXUTAFBRFM8MDAyAFzRYQETabNvOSpJ0MkBcgURRbM/lchlZlof3798PDzzwQMnxV69eRatWrWpHCGWSyeRZN7e9lKTTaYwxXrpnzx48NDT0BAD8w9NPP53x2IUZwzCWlpeXnx0aGvqw1Evs/v3LX/4yvu+++zpyudxxWZanXPjS6/kCwBJCyJne3t7xT3Lg+LEatyzLNZZlzQ6FQkd7e3vd3JKS4++55x7S2dm5AgDWQD5pyRdukiSpLZfLDYuiGKiVR11dXZtpmlfT6fR1L1eh6B4WZ7PZs5qm5YLg/DzPL81kMmf6+/snAcDTsFeuXIljsdiikZGR0zt27MjdhGb8juzcuZMZGxtbevz48aGhoSEb8hmLWyCfr/47Od/d3d0RSun8q1evHj948KAvt/iGDRv4+fPnL04kEoFqKJPJZNS27XrHcY4pipILcrZx20OBlmV9yDAMhMPhBfF4/Mh0UGAJ4/jTaDT6F4888si+5uZmL2OllFKWUtpCKb3EMMx4gB016jjOLITQRwihcZ97oo7jhDHGLQihk5RSvxM1SimtpJSWA8ANrxPEIsSlglK6kGGYowEMyaGU1jqOE2UYZgQhNAWQZ5IKhULk6tWrofLycrusrMzJZrMMzSfQVAPAPErpYYZhGJZlKUIITNOcLv3DcRynHgCiGOMPEEJ+R+nEcZw5ABBFCB3PZrOng7QDv+13bkopiUajs23bniw27ADyKCFk9vj4eI/XIUvBr60oPLj/kGXZt+jhscce45ctW7aqt7f3jYCuVDWldHYikfjXIOMTiUQsEolUJ5PJnwUZL4riXQzDhARB+EGQ8bquN1JK7WQy+ZrP0D8CgHc1TatxHOdm/Q2F53hmGqRmuW3b10VR/K8g8+nr65vrOM4NQRB+fieAIL4B5TPPPOP6zMsIIWYymfwgQE0kQP7UrBsAMpDPm55W3CR2SZLmcRzXODExcdLLsN3PVlV19t13373o4sWLb3lNwj0yl2W52bKsatM0jwWB7kRRnB+NRmcnk8m3A744bZTSyKJFi454QXFuYCrL8kLbtplEInEwgPrhmpqanS+99FKHIAjv3aS/vgAHuhAocBwHiqK0EEKuiaL4Ky/6YjfnRlGUZaZpEkEQDgXtUHfbG/emTZvgmWee4W3bPpvNZj8s+GR+l1VBvjDgdQD4SalBhmHAs88+C6IozmFZNpxMJo+4zJ7TyQsvvODSRizAGJf39PQMNTY2elIfP/nkk6DreqNt25YkSSez2awXFRhMTk4iTdMWhsNhHI/Hj/lBX48++ijSNK0dITQuy/L5wcFB6hWohcNhLElSi+M4piAIvvWfqVSKMQyjYs2aNfrRo0fnAADnpT8Wi7GJRGIpx3FjqVRqGACmpS92kZjR0VFG1/XPOI7zgSAIl2+hQ93tL+7OXeoQx/W5iypxaiB/KFEdBHM1DKNF1/WFAAAutOa146mqukKW5VkuXOj14nzrW99CkiR1uFUng4OlOWNcck1VVZd3d3f76u/r64Pt27djVVVXaZpW6beORd8g94iiGPXDrwu9iFhN0+5XFIX3iJGqoFAgLEkSK4piR6F3j+95wNe+9jVGUZQHXO7F/0uCnE+kz71p06Yg/jjMnTs3B3nWpy9A/tRxyicQ5FRVXWZZ1hlBEDIAAK2trSXHnzhxgtE0bRUh5JggCOMuXOiBwVfU19cvZFn2THd390cAAB0dJSvh4Ny5c5xhGJ+bmpp6y0VEvPRTSqtzuVwdpXQomUz6IjqnT58u0zTtc6Zpvu7yYHshELZt11VWVlZijA8nEgkXESkO7pIA8L0CsgLxeDzGcdyK8fHxn/X19QU5D6gBgDmhUOhgTU0N3K6IyMeGlpRKnFq9ejW+du3a5kwmU9ff3/+PFy5c4DxQCAoALMdxKwkhB508RZEnwsGybBnkCXvO+nVUo5RCJBJhKaX3UkoPuCQ7XpdgjKMMwzTYtn2CUkoCjI8wDHN3Lpf7H6+KlqLxsxiGmWea5hGEEOsHz3AcF8MYN2Wz2V9Mpx8hRDOZDPvd7363t7y8fM83v/nNqxjjOblc7j2MsR+7Fy2k5dYyDPP6hQsX7F27dt2RXsfHkTiFDxw4sBkAhkRRfPbSpUuoFOsaQojath1jWbaG47gfZ7NZTx5pjDE4jlOHECrHGO+3LIv4lFMRjHG9ZVkVCKEf+/FUY4ypZVnzWZZlCSFvIoQ88WhCCGEYZh6llM/lcq96YfsuFMcwTIvjOCYh5E0/NjpKKWVZdiHGOGdZ1k+99NfW1sL4+PhTHR0df/7aa6+xGzZsePPEiROeL1o2m6VlZWXNADCKEHqP53nnTjXsW5ISPncl5Fvd/bHXte7JWTqdnqsoyuJt27YhL/KX4t7tQaquXd9SUZRGd3xA/U3pdLrBzz998cUX3YBzkaZpzQAAr7zyiu+aybK8WFXVOX7+bIGqmFEUpU3TtPlBnkdhDT9jGEYtAMDSpd7F7z09Payu60tkWZ4D4H/I9Wk37rpC4Bio5F/X9XskSWoMErh8//vfh3Q63WwYxjKvccVcGKqq3pdOp2uL4TwPVATJsny3ruuLAUonVhXPU1XVVf39/RVB7jWdTmNd11e6bURK3W/R35GmafcIghAJ8hKnUilGkqT7+vv7o0EMVZZlRlGUpU8//XR0xpL9jft+AND8EBEXk9Y0rVVRlHo/w9Y0DXbu3Amqqv4aUQhCcK+qaruqqhV++mVZBsMwQrIsL5UkKVqM8/q9mIqihAB+u5XcdPM3DKNCluWWdDodCfoiKIqycuPGjWyQF98wjBpJkhb09/d7VjetXbvWveewqqqrt27diotRoU+D3FJnhYL/+lcFNyQB/p3BwpqmLSaEXEqlUldKReTuAx0ZGWFv3Lix3HGcIUmSbvhF/LquR2VZbslms+cKLLGe+imlIcuylpimedzNWZmOOcktONA0LSbLclMmkzlSaPH9676bN8OCAACO45QTQpp27dp1Kp1Ol0SLfvjDH7ov5V0MwyxIpVLvxuNx22/+iqJUEULqRFE8PzIykvOCWXfv3g09PT1VjuPMSyaTB770pS8RgEAtVz59aMmOHTvmSZL0eCaTeR/ybE2e0tnZGa6qqmqzLGtIURTfo3RBEGaFQqGFQQjMCztSFca4gef5I11dXb49LEVRjCGEGtLp9HsBd9S5GONZExMTx4I0EO3r65ttmmZMEIRjAeffTCkFURRPBxkvSdJClmXZRCJxMsh4RVEWEUKoq/92zOr7/zLuCMa4q6Gh4ejmzZvfmZycZMAjKQghNB8AQo7jnMEYexoeIYRijOcXkJGrCCErwJybMMambdsf+JHRFObZgBCilNJLhd+px3wAY9yKEMpQSoM2pmokhGQxxpf91pRSCizLtjuOcxlK0FbcJAyltBFjPF6Yj/8NU7oYAIYJIecFQRiDT7lUQ77caXYRRNgO+dNGgDxHX3NQZYqizAk6dv369YyiKNVBx+/duxcEQagPOl4URXbLli2BA6l4PM709fXVBR2fSqXYRCIRuHPE9u3bOUmSqoKOFwSBbW1tRbcwn7CiKNHpguJPq1QAwKsA8DP4DU3AzwHgL2/VfXEDwP/r5JsgD614zK2yjAYJZIv136oRBenoVqwzCHQ3Y8il5X7IH5k3AsCTkC8qYGaWZUbuFHkdAP4b8qVgMzIjd5R8D/LddefOLMWM3O6Ci36KkCdmzwDAX88szYzcCRKGfOfch9y4BAAOzyzLjNwJ8jDkW0mXF36vhzztQk/R32ZkRm47+V+oPh1foYtE0gAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TD6: Solving the spherical laplacian using spectral methods\n",
    "\n",
    "The goal of this TD is to solve a **laplacian problem** in **spherical geometry** with Dirichlet boundary conditions using **spectral methods**. The spatial domain $\\Omega$ of this problem is the spherical shell bounded by two concentric spheres $\\mathcal{S}(0, R_1)$ (sphere of radius $R_1$ centered at the origin) and $\\mathcal{S}(0, R_2)$ with $R_1 < R_2$. Thus we want to solve:\n",
    "$$\n",
    "   \\begin{cases}\n",
    "        \\Delta f (r, \\theta, \\varphi) = 0 \\text{ with } (r, \\theta, \\varphi) \\in \\Omega\\\\\n",
    "        f (R_1, \\theta, \\varphi) = g_1 (\\theta, \\varphi) \\\\\n",
    "        f (R_2, \\theta, \\varphi) = g_2 (\\theta, \\varphi)\n",
    "   \\end{cases}\n",
    "$$\n",
    "with $g_1$ and $g_2$ two boundary functions. To simplify this problem, we will assume **$\\varphi$ invariance** so we can forget about this variable and we will use the following boundary functions:\n",
    "$$\n",
    "   \\begin{cases}\n",
    "        g_1 (\\theta) = \\cos (\\theta)\\\\\n",
    "        g_2 (\\theta) = 0\n",
    "   \\end{cases}\n",
    "$$\n",
    "\n",
    "We will adopt the following convention to define spherical coordinates $(r, \\theta, \\varphi)$ in a cartesian system:\n",
    " <div class=\"row\">\n",
    "  <div class=\"column\" style=\"float:left; margin:0; width:50%;\">$$\n",
    "\\begin{align}&\\\\&\\\\ x&=r \\, \\sin\\theta \\, \\cos\\varphi \\\\ y&=r \\, \\sin\\theta \\, \\sin\\varphi \\\\ z&=r \\, \\cos\\theta\\end{align}\n",
    "$$ </div>\n",
    "  <div class=\"column\" style=\"float:left; margin:0; width:50%;\">![3D_Spherical.png](attachment:3D_Spherical.png)</div>\n",
    "</div> \n",
    "\n",
    "$(r, \\theta, \\varphi)$ variables should satisfy: $r \\geq 0$, $\\theta \\in [0, \\pi]$ and $\\varphi \\in [0, 2\\pi[$.\n",
    "\n",
    "\n",
    "We will need the usual Python packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To draw matplotlib plots within this notebook.\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Spherical harmonics\n",
    "\n",
    "[Spherical harmonics](https://en.wikipedia.org/wiki/Spherical_harmonics) are a set of angular functions which can be defined by:\n",
    "$$\n",
    "Y^m_l(\\theta,\\varphi) = \\sqrt{\\frac{2l+1}{4\\pi} \\frac{(l-m)!}{(l+m)!}}\n",
    "  e^{i m \\varphi} P^m_l(\\cos(\\theta))\n",
    "$$\n",
    "with $P_l^m$ an [associated Legendre polynomial](https://en.wikipedia.org/wiki/Associated_Legendre_polynomials). These functions are eigenfunctions of the angular laplacian operator $\\Delta_{\\text{ang}}$:\n",
    "$$\n",
    "    \\Delta_{\\text{ang}}Y^m_l = \\left(\\frac{\\partial^2}{\\partial \\theta^2} + \\frac1{\\tan \\theta}\\frac{\\partial}{\\partial \\theta}+\\frac1{\\sin^2\\theta}\\frac{\\partial^2}{\\partial\\varphi^2}\\right)Y_l^m = -l(l+1)Y_l^m\n",
    "$$\n",
    "\n",
    "In Python, we can manipulate these functions using [scipy.special.sph_harm](https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.special.sph_harm.html) routine. We will be very careful to the fact that the documentation page of this routine uses another convention for the variables: $\\theta$ and $\\varphi$ variables are inverted.\n",
    "\n",
    "**1)** Plot the real component of a few spherical harmonics in 2D using [pcolor](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.pcolor.html) routine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import sph_harm\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For information, we can also have a 3D view of these functions using the following routine:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from matplotlib import cm\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "\n",
    "def plot_3d_sph(l, m, res=30):\n",
    "    theta = np.linspace(0, np.pi, res)\n",
    "    phi = np.linspace(0, 2 * np.pi, res)\n",
    "    theta, phi = np.meshgrid(theta, phi)\n",
    "\n",
    "    # The Cartesian coordinates of the unit sphere\n",
    "    x = np.sin(theta) * np.cos(phi)\n",
    "    y = np.sin(theta) * np.sin(phi)\n",
    "    z = np.cos(theta)\n",
    "\n",
    "    # Calculate the spherical harmonic Y(l,m) and normalize to [0,1]\n",
    "    fcolors = sph_harm(m, l, phi, theta).real\n",
    "    fmax, fmin = fcolors.max(), fcolors.min()\n",
    "    fcolors = (fcolors - fmin) / (fmax - fmin)\n",
    "\n",
    "    # Set the aspect ratio to 1 so our sphere looks spherical\n",
    "    fig = plt.figure(figsize=plt.figaspect(1.))\n",
    "    ax = fig.add_subplot(111, projection='3d')\n",
    "    m = cm.ScalarMappable(cmap=cm.seismic)\n",
    "    m.set_array(fcolors)\n",
    "    ax.plot_surface(x, y, z, rstride=1, cstride=1, facecolors=cm.seismic(fcolors))\n",
    "    ax.set_axis_off()\n",
    "    plt.colorbar(m)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any square-integrable function $f(\\theta, \\phi)$ can be expanded as a linear combination of spherical harmonics:\n",
    "$$\n",
    "    f(\\theta, \\varphi) = \\sum_{l=0}^\\infty \\sum_{m=-l}^l f_l^m \\, Y_l^m(\\theta, \\varphi)\n",
    "$$\n",
    "The coefficients can be computed through the following relation:\n",
    "$$\n",
    "f_l^m=\\int_{\\Omega} f(\\theta,\\varphi)\\, Y_l^{m*}(\\theta,\\varphi)\\,d\\Omega = \\int_0^{2\\pi}d\\varphi\\int_0^\\pi \\,d\\theta\\,\\sin\\theta~f(\\theta,\\varphi)~Y_l^{m*} (\\theta,\\varphi)\n",
    "$$\n",
    "\n",
    "**2)** Write a function which computes for a given function $f$ and given integers $l$ and $m$ the associated complex coefficient $f_l^m$. We can use the routine [scipy.integrate.dblquad](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.dblquad.html#scipy.integrate.dblquad) to compute this double integral. It is recommended to have a look at the [scipy.integrate](https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html) documentation page to see an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import dblquad\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: A spectral approach\n",
    "\n",
    "We want to solve the previous Laplacian equation using spectral methods. We will approximate the solution by an expansion over **Chebyshev polynomials for the radial variable $r$**, and over **spherical harmonics for the angular variable $\\theta$**. Our numerical solution should take the following form:\n",
    "$$\n",
    "    f^{\\text{appr}}(r, \\theta) = \\sum_{l=0}^L\\sum_{n=0}^Nf_{l,n}T_n(u(r))Y_l^0(\\theta) \\text{ with } u(r) = \\frac{r - (R_1+R_2)/2}{(R_2 - R_1)/2}\n",
    "$$\n",
    "$u$ function plays the role of a new variable satisfying $u \\in [-1, 1]$.\n",
    "\n",
    "To solve our problem, we will work in **spectral space**. Therefore, the goal is to **find $\\{f_{l,n}\\}_{0 \\leq l \\leq L, 0 \\leq n \\leq N}$ coefficients** so that $f^{\\text{appr}}$ can  satisfy the problem equation as well as the boundary conditions in an ideal case. But as there is no reason for the exact solution of this problem to be of the prescribed form, we will eventually need to make reasonable compromises so that $f^{\\text{appr}}$ verifies the boundary conditions.\n",
    "\n",
    "Let us define $L$ and $N$ variables with Python with some intial values. We are strongly encouraged to change these values during this TD."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 20\n",
    "L = 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3)** We recall the expression of the **scalar laplacian operator** in spherical coordinates:\n",
    "$$\n",
    "    \\Delta = \\frac{\\partial^2}{\\partial r^2} + \\frac{2}{r}\\frac{\\partial}{\\partial r} + \\frac1{r^2}\\Delta_{\\text{ang}}\n",
    "$$\n",
    "Show that in $\\Omega$:\n",
    "$$\n",
    "\\begin{align}\n",
    "    &\\Delta f^{\\text{appr}}(u, \\theta)=0\\\\ \\iff &\\sum_{l=0}^L\\sum_{n=0}^Nf_{l,n}\\left( u^2T_n^{\\prime\\prime}(u) + 2\\frac{R_1 + R_2}{R_2 - R_1}uT_n^{\\prime\\prime}(u) + \\left(\\frac{R_1 + R_2}{R_2 - R_1}\\right)^2T_n^{\\prime\\prime}(u) + 2uT_n^\\prime(u) + 2\\frac{R_1+R_2}{R_2-R_1}T_n^\\prime(u) - l(l+1)T_n(u)\\right)Y_l^0(\\theta)=0\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we should be able to rewrite the previous equation as:\n",
    "$$\n",
    "    \\sum_{l=0}^L\\sum_{n=0}^Ng_{l,n}T_n(u)Y_l^0(\\theta)=0\n",
    "$$\n",
    "introducing some new coefficients $\\{g_{l,n}\\}$ which must be linearly related to $\\{f_{l, n}\\}$. Actually, for a given $l$, there exists a matrix $A_l$ of shape $(N + 1)\\times(N+1)$ so that:\n",
    "$$\n",
    "    G_l = A_lF_l \\text{ with } F_l = (f_{l,n})_{0\\leq n \\leq N} \\text{ and } G_l = (g_{l,n})_{0\\leq n\\leq N}\n",
    "$$\n",
    "Provided we can numerically build these matrix $A_l$, we have the equivalence:\n",
    "$$\n",
    "    \\Delta f^{\\text{appr}}(u, \\theta) = 0 \\iff \\forall l \\in \\{0,...,L\\}, G_l = A_lF_l = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**4) Building Chebyshev differentiation matrix:** We first want to build a matrix of differentiation $D$ of shape $(N+1)\\times(N+1)$ in Chebyshev space. With $X_j=(\\delta_{i,j})_{0\\leq i\\leq N}$ the basis vector which represents $T_j$ polynomial in Chebyshev space, build $D$ so that $DX_j$ is the vector which represents $T_j^\\prime$ polynomial coordinates on the basis of Chebyshev polynomials $T_0, ..., T_N$.\n",
    "\n",
    "We will need the following result: for $h(x) = \\sum_{n=0}^{+\\infty}h_nT_n(x)$ an expansion over Chebyshev polynomials, with $h^\\prime(x)=\\sum_{n=0}^{+\\infty}h_n^\\prime T_n(x)$ we have:\n",
    "$$\n",
    "    h_n^\\prime=\\frac{2}{c_n}\\sum_{\\substack{p=n+1\\\\p + n \\text { odd}}}^{+\\infty}ph_p \\text{ where } c_0 = 2 \\text{ and }c_n = 1 \\text{ for } n\\geq 1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5) Building Chebyshev $u$ multiplication matrix**: Now we want to build a matrix $M$ of shape $(N+1)\\times(N+1)$ so that (using the same notations as the previous question) $MX_j$ is the vector that represents $uT_j(u)$ polynomial coordinates on the basis of Chebyshev polynomials $T_0, ..., T_N$.\n",
    "\n",
    "We will need the following result for this question:\n",
    "$$\n",
    "    \\forall m, n\\in\\mathbb{N}, T_m(u) T_n(u) = \\frac1{2}\\left(T_{m+n}(u) + T_{|m-n|}(u)\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**6)** What is the expression of $A_l$ as a function of $M$ and $D$ matrix, and as a function of $I_{N+1}$ matrix (which is the $(N+1)\\times(N+1)$ identity matrix)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: Enforcing boundary conditions\n",
    "\n",
    "Solving our laplacian problem demands that $f^{\\text{appr}}$ satisfies our boundary conditions:\n",
    "$$\n",
    "\\begin{cases}\n",
    "    f^{\\text{appr}} (R_1, \\theta) = g_1 (\\theta) \\\\\n",
    "    f^{\\text{appr}} (R_2, \\theta) = g_2 (\\theta)\n",
    "\\end{cases}\n",
    "$$\n",
    "We will enforce boundary conditions for $f^{\\text{appr}}$ using **Tau method**. We have:\n",
    "$$\n",
    " \\begin{cases}\n",
    "     f^{\\text{appr}}(R_1, \\theta) = \\sum_{l=0}^L\\sum_{n=0}^Nf_{l,n}(-1)^nY_l^0(\\theta)\\\\\n",
    "     f^{\\text{appr}}(R_2, \\theta) =\n",
    "     \\sum_{l=0}^L\\sum_{n=0}^Nf_{l,n}Y_l^0(\\theta)\n",
    " \\end{cases}\n",
    "$$\n",
    "\n",
    "**7)** With $g_1 (\\theta) = \\sum_{l=0}^Lg_{1,l}Y_l^0(\\theta)$ and $g_2 (\\theta) = \\sum_{l=0}^Lg_{2,l}Y_l^0(\\theta)$, what are the linear relations that should satisfy $\\{f_{l,n}\\}$ coefficients to enforce the boundary conditions?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**8)** Change $A_l$ matrix last two rows (these rows involve relations between higher order Chebyshev coefficients) and define a right-hand member $B_l$ so that the previous relations can be satisfied as well as most of the relations of part 2 through the system:\n",
    "$$\n",
    "    \\forall l\\in\\{0, ...,L\\}, A_lF_l = B_l\n",
    "$$\n",
    "\n",
    "With $g_1 (\\theta) = \\cos(\\theta)$ and $g_2 (\\theta)=0$, to get $\\{g_{1,l}\\}$ and $\\{g_{2,l}\\}$ coefficients we can either use the function we should have programmed in question 2) or guess the spherical harmonics decomposition of these simple functions..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**9)** Invert these linear systems using for example [scipy.linalg.lu_factor](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu_factor.html) and [scipy.linalg.lu_solve](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu_solve.html) routines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.linalg as la\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**10)** Now plot the numerical solution $f^{\\text{appr}}$ using [pcolor](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.pcolor.html) routine. We can use [numpy.polynomial.chebyshev.chebval](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.polynomial.chebyshev.chebval.html#numpy.polynomial.chebyshev.chebval) function to compute the values of a linear combination of Chebyshev polynomials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy.polynomial.chebyshev as cheb\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**11)** Compare the previous numerical solution to the exact solution of this problem:\n",
    "$$\n",
    "    f(r, \\theta) = (r/R_1)\\frac{\\left(R_1/r\\right)^3 - a^3}{1-a^3}\\cos (\\theta) \\text{ with } a = \\frac{R_1}{R_2}\n",
    "$$\n",
    "We can try several values of $N$ or $L$ and observe the difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**BONUS: 12)** We now want to check that our code is still working for boundary conditions involving multiple spherical harmonics. Numerically solve the same laplacian problem with:\n",
    "$$\n",
    "   \\begin{cases}\n",
    "        g_1 (\\theta) = \\exp\\left(\\cos\\left(\\theta\\right)\\right)\\\\\n",
    "        g_2 (\\theta) = 0\n",
    "   \\end{cases}\n",
    "$$\n",
    "We will then compare this numerical solution to the exact solution:\n",
    "$$\n",
    "    f(r,\\theta) = \\sum_{l=0}^{+\\infty}(r/R_1)^l\\frac{(R_1/r)^{2l+1}-a^{2l+1}}{1 - a^{2l+1}}g_lY_l^0(\\theta)\n",
    "$$\n",
    "with $\\{g_l\\}$ the spherical harmonics coefficients of $g_1$, meaning $g_1(\\theta)=\\sum_{l=0}^{+\\infty}g_lY_l^0(\\theta)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}