testing_cylinder_fdm.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cylinder_fdm_3d import CylinderFDM\n",
    "from icecream import ic\n",
    "import numpy as np\n",
    "from scipy.sparse.linalg import eigsh\n",
    "from time import time\n",
    "from matplotlib import pyplot as plt\n",
    "from fft_tdse.simulator import LaserPulse\n",
    "from erfgau import ErfgauPotential\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ic| self.m_i: array([0.])\n",
      "ic| cylinder_fdm_3d.py:169 in __init__() at 22:04:00.890\n",
      "ic| self.n_r: 300, self.n_z: 300, self.n_m: 1\n",
      "ic| self.r_max: 20, self.z_max: 20\n",
      "ic| self.n_dof: 90000\n",
      "ic| tt.shape: (1, 300, 300)\n",
      "    rr.shape: (1, 300, 300)\n",
      "    zz.shape: (1, 300, 300)\n",
      "ic| m: 0.0, N: 41.56142818209105\n",
      "ic| V_FFT[i, :, :]: array([[-0.05033417+0.j, -0.05067311+0.j, -0.05101666+0.j, ...,\n",
      "                            -0.05101666+0.j, -0.05067311+0.j, -0.05033417+0.j],\n",
      "                           [-0.05033332+0.j, -0.05067225+0.j, -0.05101578+0.j, ...,\n",
      "                            -0.05101578+0.j, -0.05067225+0.j, -0.05033332+0.j],\n",
      "                           [-0.05033191+0.j, -0.05067082+0.j, -0.05101431+0.j, ...,\n",
      "                            -0.05101431+0.j, -0.05067082+0.j, -0.05033191+0.j],\n",
      "                           ...,\n",
      "                           [-0.0356514 +0.j, -0.03577123+0.j, -0.03589147+0.j, ...,\n",
      "                            -0.03589147+0.j, -0.03577123+0.j, -0.0356514 +0.j],\n",
      "                           [-0.03559183+0.j, -0.03571106+0.j, -0.03583069+0.j, ...,\n",
      "                            -0.03583069+0.j, -0.03571106+0.j, -0.03559183+0.j],\n",
      "                           [-0.03553236+0.j, -0.035651  +0.j, -0.03577002+0.j, ...,\n",
      "                            -0.03577002+0.j, -0.035651  +0.j, -0.03553236+0.j]])\n",
      "ic| m_max: 0\n",
      "ic| m: 0.0\n",
      "ic| m_list[i]: 0.0, np.linalg.norm(V_m[i]): 41.56142818209105\n"
     ]
    }
   ],
   "source": [
    "\n",
    "n_r = 300\n",
    "n_z = 300\n",
    "n_m = 1\n",
    "solver = CylinderFDM(r_max = 20, z_max = 20, n_r = n_r, n_z = n_z , n_m = n_m)\n",
    "\n",
    "tt, rr, zz = solver.get_trz_meshgrid()\n",
    "ic(tt.shape, rr.shape, zz.shape)\n",
    "xx = rr*np.cos(tt)\n",
    "yy = rr*np.sin(tt)\n",
    "\n",
    "\n",
    "erfgau = ErfgauPotential(mu=1.0)\n",
    "V = erfgau.potential_radial((rr*rr + zz*zz)**.5)\n",
    "D = zz\n",
    "\n",
    "solver.set_realspace_potential(V, rotation_symmetric=False)\n",
    "\n",
    "#laser = LaserPulse(t0 = 0, omega = 0.057, T = 100, E0 = 0.05)\n",
    "#solver.set_td_potential(D, rotation_symmetric=False)\n",
    "#solver.set_td_potential_modulator(laser)\n",
    "\n",
    "#ic(solver.V_m[0][0,0], solver.V[0,0,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ic| 'Initial guess energy:'\n",
      "ic| E: -0.3960477716838697\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t = 0.00, E = -0.41387932011853495, delta = 0.021117285849218622\n",
      "t = 0.05, E = -0.4255404703765809, delta = 0.01707313307752508\n",
      "t = 0.10, E = -0.4338244918522796, delta = 0.014387806817378434\n",
      "t = 0.15, E = -0.44010177720359794, delta = 0.012523407136044058\n",
      "t = 0.20, E = -0.4451081333097773, delta = 0.011183804464786343\n",
      "t = 0.25, E = -0.4492643720955761, delta = 0.010190356507179296\n",
      "t = 0.30, E = -0.45282309010381017, delta = 0.009429984010004186\n",
      "t = 0.35, E = -0.45594237257818543, delta = 0.008829208036364163\n",
      "t = 0.40, E = -0.45872480902282625, delta = 0.008339529851609968\n",
      "t = 0.45, E = -0.46123939480840637, delta = 0.00792857691228756\n",
      "t = 0.50, E = -0.4635340673478557, delta = 0.007574523760984364\n",
      "t = 0.55, E = -0.4656432981122005, delta = 0.007262504678452566\n",
      "t = 0.60, E = -0.4675926487451243, delta = 0.006982281926418636\n",
      "t = 0.65, E = -0.46940170901650785, delta = 0.006726719145995892\n",
      "t = 0.70, E = -0.4710859092441204, delta = 0.006490775498352808\n",
      "t = 0.75, E = -0.47265776286518285, delta = 0.006270838504982709\n",
      "t = 0.80, E = -0.47412764228074233, delta = 0.006064278349024609\n",
      "t = 0.85, E = -0.475504348004221, delta = 0.005869147826531243\n",
      "t = 0.90, E = -0.4767954634794311, delta = 0.005683978718138397\n",
      "t = 0.95, E = -0.47800763688787723, delta = 0.005507642450030835\n",
      "t = 1.00, E = -0.4791467545256783, delta = 0.0053392539465392935\n",
      "t = 1.05, E = -0.4802180919022131, delta = 0.0051781047275124776\n",
      "t = 1.10, E = -0.48122640475511835, delta = 0.005023615963511617\n",
      "t = 1.15, E = -0.48217601691091505, delta = 0.004875305257898811\n",
      "t = 1.20, E = -0.48307087201363874, delta = 0.004732762943087751\n",
      "t = 1.25, E = -0.48391458880381455, delta = 0.0045956350208572385\n",
      "t = 1.30, E = -0.4847104931143149, delta = 0.004463610776574126\n",
      "t = 1.35, E = -0.48546165521799356, delta = 0.004336413704970076\n",
      "t = 1.40, E = -0.4861709113242153, delta = 0.00421379479868132\n",
      "t = 1.45, E = -0.48684089035963596, delta = 0.004095527534235139\n",
      "t = 1.50, E = -0.48747402949832225, delta = 0.0039814040858555284\n",
      "t = 1.55, E = -0.4880725942790209, delta = 0.003871232433479392\n",
      "t = 1.60, E = -0.4886386904883471, delta = 0.0037648341265482292\n",
      "t = 1.65, E = -0.4891742798032984, delta = 0.00366204253215106\n",
      "t = 1.70, E = -0.48968118927182325, delta = 0.0035627014435651518\n",
      "t = 1.75, E = -0.490161123782762, delta = 0.003466663959066767\n",
      "t = 1.80, E = -0.49061567384897536, delta = 0.0033737915651130056\n",
      "t = 1.85, E = -0.49104632572697804, delta = 0.003283953375452692\n",
      "t = 1.90, E = -0.4914544679286043, delta = 0.0031970254903550305\n",
      "t = 1.95, E = -0.4918413995394224, delta = 0.003112890449335474\n",
      "t = 2.00, E = -0.49220833573434164, delta = 0.0030314367574754138\n",
      "t = 2.05, E = -0.49255641468040595, delta = 0.002952558470360805\n",
      "t = 2.10, E = -0.49288670224678954, delta = 0.0028761548263029177\n",
      "t = 2.15, E = -0.4932001977740995, delta = 0.0028021299171950963\n",
      "t = 2.20, E = -0.4934978381184561, delta = 0.002730392391365001\n",
      "t = 2.25, E = -0.49378050250090505, delta = 0.0026608551832772334\n",
      "t = 2.30, E = -0.49404901599292955, delta = 0.002593435266065872\n",
      "t = 2.35, E = -0.4943041536117748, delta = 0.002528053423723838\n",
      "t = 2.40, E = -0.494546643332903, delta = 0.002464634040418423\n",
      "t = 2.45, E = -0.49477716956240464, delta = 0.002403104904892712\n",
      "t = 2.50, E = -0.49499637574637956, delta = 0.0023433970282874476\n",
      "t = 2.55, E = -0.49520486732547725, delta = 0.0022854444740080552\n",
      "t = 2.60, E = -0.49540321399939863, delta = 0.0022291841984866563\n",
      "t = 2.65, E = -0.4955919522488764, delta = 0.0021745559018654017\n",
      "t = 2.70, E = -0.4957715873040544, delta = 0.0021215018877664316\n",
      "t = 2.75, E = -0.4959425953036458, delta = 0.0020699669314253374\n",
      "t = 2.80, E = -0.49610542500853017, delta = 0.002019898155553845\n",
      "t = 2.85, E = -0.49626049965543, delta = 0.001971244913371138\n",
      "t = 2.90, E = -0.4964082184509139, delta = 0.001923958678302395\n",
      "t = 2.95, E = -0.49654895816699385, delta = 0.0018779929398934274\n",
      "t = 3.00, E = -0.4966830744455048, delta = 0.0018333031055325815\n",
      "t = 3.05, E = -0.4968109031750749, delta = 0.001789846407606701\n",
      "t = 3.10, E = -0.49693276163154615, delta = 0.0017475818157494123\n",
      "t = 3.15, E = -0.4970489496691202, delta = 0.0017064699538659963\n",
      "t = 3.20, E = -0.4971597507188463, delta = 0.0016664730216445496\n",
      "t = 3.25, E = -0.49726543282142743, delta = 0.0016275547202826953\n",
      "t = 3.30, E = -0.4973662495026454, delta = 0.0015896801821793588\n",
      "t = 3.35, E = -0.4974624406709125, delta = 0.0015528159043570952\n",
      "t = 3.40, E = -0.4975542333858355, delta = 0.0015169296853975171\n",
      "t = 3.45, E = -0.49764184263995404, delta = 0.0014819905656853298\n",
      "t = 3.50, E = -0.49772547203435313, delta = 0.0014479687707697792\n",
      "t = 3.55, E = -0.4978053144610034, delta = 0.0014148356576656365\n",
      "t = 3.60, E = -0.49788155269752826, delta = 0.0013825636639251014\n",
      "t = 3.65, E = -0.4979543600039268, delta = 0.0013511262593239892\n",
      "t = 3.70, E = -0.4980239006468286, delta = 0.0013204979000141224\n",
      "t = 3.75, E = -0.4980903304224455, delta = 0.0012906539850038494\n",
      "t = 3.80, E = -0.4981537971193091, delta = 0.0012615708148357699\n",
      "t = 3.85, E = -0.4982144409773736, delta = 0.0012332255523398385\n",
      "t = 3.90, E = -0.49827239509700566, delta = 0.0012055961853462034\n",
      "t = 3.95, E = -0.4983277858428128, delta = 0.0011786614912500846\n",
      "t = 4.00, E = -0.49838073320552545, delta = 0.001152401003326356\n",
      "t = 4.05, E = -0.49843135115789916, delta = 0.001126794978698386\n",
      "t = 4.10, E = -0.4984797479753089, delta = 0.0011018243678708592\n",
      "t = 4.15, E = -0.4985260265499175, delta = 0.0010774707857415305\n",
      "t = 4.20, E = -0.4985702846750834, delta = 0.0010537164840123397\n",
      "t = 4.25, E = -0.49861261532309714, delta = 0.0010305443249241984\n",
      "t = 4.30, E = -0.49865310689786563, delta = 0.001007937756245155\n",
      "t = 4.35, E = -0.4986918434807854, delta = 0.0009858807874444562\n",
      "t = 4.40, E = -0.4987289050554965, delta = 0.0009643579669901214\n",
      "t = 4.45, E = -0.4987643677259566, delta = 0.0009433543607105609\n",
      "t = 4.50, E = -0.49879830391647306, delta = 0.0009228555311641777\n",
      "t = 4.55, E = -0.49883078256549185, delta = 0.0009028475179644586\n",
      "t = 4.60, E = -0.49886186930384196, delta = 0.000883316819010774\n",
      "t = 4.65, E = -0.4988916266270921, delta = 0.0008642503725779371\n",
      "t = 4.70, E = -0.49892011405467124, delta = 0.0008456355402202136\n",
      "t = 4.75, E = -0.49894738828330515, delta = 0.0008274600904484116\n",
      "t = 4.80, E = -0.49897350332925156, delta = 0.0008097121831399727\n",
      "t = 4.85, E = -0.4989985106650381, delta = 0.0007923803546455355\n",
      "t = 4.90, E = -0.49902245934672074, delta = 0.0007754535035564681\n",
      "t = 4.95, E = -0.4990453961360702, delta = 0.0007589208771001368\n",
      "t = 5.00, E = -0.49906736561433646, delta = 0.0007427720581318116\n",
      "t = 5.05, E = -0.4990884102915611, delta = 0.0007269969526932188\n",
      "t = 5.10, E = -0.49910857070859865, delta = 0.0007115857781101884\n",
      "t = 5.15, E = -0.4991278855348659, delta = 0.0006965290516023498\n",
      "t = 5.20, E = -0.49914639165989505, delta = 0.0006818175793803239\n",
      "t = 5.25, E = -0.49916412428093926, delta = 0.0006674424462063826\n",
      "t = 5.30, E = -0.49918111698516265, delta = 0.0006533950053963984\n",
      "t = 5.35, E = -0.4991974018282615, delta = 0.0006396668692414623\n",
      "t = 5.40, E = -0.49921300940828034, delta = 0.0006262498998297458\n",
      "t = 5.45, E = -0.49922796893625915, delta = 0.0006131362002487052\n",
      "t = 5.50, E = -0.4992423083026641, delta = 0.0006003181061505548\n",
      "t = 5.55, E = -0.4992560541408289, delta = 0.0005877881776630135\n",
      "t = 5.60, E = -0.4992692318867551, delta = 0.0005755391916298942\n",
      "t = 5.65, E = -0.4992818658363119, delta = 0.0005635641341657213\n",
      "t = 5.70, E = -0.49929397919903545, delta = 0.0005518561935100069\n",
      "t = 5.75, E = -0.4993055941496644, delta = 0.0005404087531674573\n",
      "t = 5.80, E = -0.4993167318766544, delta = 0.0005292153853207634\n",
      "t = 5.85, E = -0.49932741262864966, delta = 0.0005182698445039362\n",
      "t = 5.90, E = -0.4993376557583084, delta = 0.000507566061523982\n",
      "t = 5.95, E = -0.49934747976415833, delta = 0.0004970981376200954\n",
      "t = 6.00, E = -0.4993569023302443, delta = 0.0004868603388494425\n",
      "t = 6.05, E = -0.49936594036392384, delta = 0.00047684709068971845\n",
      "t = 6.10, E = -0.499374610031671, delta = 0.00046705297284843535\n",
      "t = 6.15, E = -0.49938292679327245, delta = 0.0004574727142706335\n",
      "t = 6.20, E = -0.4993909054342257, delta = 0.0004481011883352119\n",
      "t = 6.25, E = -0.4993985600966157, delta = 0.0004389334082327851\n",
      "t = 6.30, E = -0.49940590430851217, delta = 0.00042996452251638136\n",
      "t = 6.35, E = -0.4994129510119037, delta = 0.0004211898108179078\n",
      "t = 6.40, E = -0.49941971258927886, delta = 0.0004126046797233354\n",
      "t = 6.45, E = -0.4994262008890021, delta = 0.00040420465879950596\n",
      "t = 6.50, E = -0.49943242724936665, delta = 0.000395985396766549\n",
      "t = 6.55, E = -0.4994384025216159, delta = 0.0003879426578094553\n",
      "t = 6.60, E = -0.4994441370917594, delta = 0.0003800723180231015\n",
      "t = 6.65, E = -0.49944964090147276, delta = 0.00037237036198523667\n",
      "t = 6.70, E = -0.4994549234678639, delta = 0.00036483287945188495\n",
      "t = 6.75, E = -0.49945999390247064, delta = 0.00035745606217045604\n",
      "t = 6.80, E = -0.49946486092923603, delta = 0.000350236200805482\n",
      "t = 6.85, E = -0.4994695329017071, delta = 0.0003431696819723648\n",
      "t = 6.90, E = -0.49947401781941814, delta = 0.0003362529853750121\n",
      "t = 6.95, E = -0.49947832334351777, delta = 0.0003294826810429509\n",
      "t = 7.00, E = -0.49948245681163506, delta = 0.0003228554266639816\n",
      "t = 7.05, E = -0.4994864252521428, delta = 0.0003163679650084712\n",
      "t = 7.10, E = -0.4994902353976348, delta = 0.0003100171214419658\n",
      "t = 7.15, E = -0.4994938936979372, delta = 0.0003037998015222338\n",
      "t = 7.20, E = -0.4994974063323627, delta = 0.00029771298867762946\n",
      "t = 7.25, E = -0.4995007792215586, delta = 0.0002917537419637092\n",
      "t = 7.30, E = -0.49950401803866273, delta = 0.0002859191938946614\n",
      "t = 7.35, E = -0.4995071282201356, delta = 0.00028020654834730085\n",
      "t = 7.40, E = -0.4995101149758643, delta = 0.00027461307853408694\n",
      "t = 7.45, E = -0.49951298329905813, delta = 0.0002691361250432146\n",
      "t = 7.50, E = -0.4995157379755175, delta = 0.0002637730939427301\n",
      "t = 7.55, E = -0.49951838359256323, delta = 0.00025852145494643344\n",
      "t = 7.60, E = -0.49952092454756436, delta = 0.00025337873963943023\n",
      "t = 7.65, E = -0.4995233650560493, delta = 0.0002483425397607188\n",
      "t = 7.70, E = -0.49952570915950195, delta = 0.0002434105055410644\n",
      "t = 7.75, E = -0.4995279607327557, delta = 0.0002385803440936686\n",
      "t = 7.80, E = -0.4995301234910961, delta = 0.0002338498178563315\n",
      "t = 7.85, E = -0.4995322009970813, delta = 0.00022921674308250343\n",
      "t = 7.90, E = -0.4995341966669648, delta = 0.0002246789883800159\n",
      "t = 7.95, E = -0.4995361137769134, delta = 0.00022023447329541296\n",
      "t = 8.00, E = -0.4995379554689303, delta = 0.00021588116694237227\n",
      "t = 8.05, E = -0.4995397247565201, delta = 0.0002116170866726395\n",
      "t = 8.10, E = -0.49954142453009664, delta = 0.00020744029678795888\n",
      "t = 8.15, E = -0.4995430575621586, delta = 0.00020334890729141507\n",
      "t = 8.20, E = -0.49954462651222925, delta = 0.0001993410726770585\n",
      "t = 8.25, E = -0.4995461339316162, delta = 0.00019541499075615093\n",
      "t = 8.30, E = -0.49954758226792206, delta = 0.0001915689015189662\n",
      "t = 8.35, E = -0.4995489738693818, delta = 0.00018780108603086337\n",
      "t = 8.40, E = -0.49955031098901076, delta = 0.00018410986536137152\n",
      "t = 8.45, E = -0.4995515957885647, delta = 0.000180493599545315\n",
      "t = 8.50, E = -0.4995528303423134, delta = 0.0001769506865745911\n",
      "t = 8.55, E = -0.4995540166407601, delta = 0.00017347956141984807\n",
      "t = 8.60, E = -0.4995551565940021, delta = 0.00017007869508093735\n",
      "t = 8.65, E = -0.4995562520351262, delta = 0.0001667465936651521\n",
      "t = 8.70, E = -0.4995573047233805, delta = 0.00016348179749212812\n",
      "t = 8.75, E = -0.4995583163472348, delta = 0.00016028288022514154\n",
      "t = 8.80, E = -0.4995592885272353, delta = 0.00015714844802698367\n",
      "t = 8.85, E = -0.4995602228189142, delta = 0.00015407713874049565\n",
      "t = 8.90, E = -0.49956112071534214, delta = 0.00015106762109235914\n",
      "t = 8.95, E = -0.4995619836497889, delta = 0.00014811859391964253\n",
      "t = 9.00, E = -0.4995628129981257, delta = 0.00014522878541824968\n",
      "t = 9.05, E = -0.49956361008118627, delta = 0.00014239695241272674\n",
      "t = 9.10, E = -0.49956437616701665, delta = 0.00013962187964648132\n",
      "t = 9.15, E = -0.49956511247302643, delta = 0.00013690237909186022\n",
      "t = 9.20, E = -0.499565820168076, delta = 0.00013423728927970162\n",
      "t = 9.25, E = -0.49956650037443545, delta = 0.0001316254746472689\n",
      "t = 9.30, E = -0.49956715416969644, delta = 0.00012906582490428113\n",
      "t = 9.35, E = -0.49956778258858653, delta = 0.00012655725441658788\n",
      "t = 9.40, E = -0.4995683866246548, delta = 0.00012409870160650133\n",
      "t = 9.45, E = -0.4995689672320282, delta = 0.00012168912836973985\n",
      "t = 9.50, E = -0.49956952532698323, delta = 0.00011932751950812038\n",
      "t = 9.55, E = -0.499570061789428, delta = 0.00011701288217763363\n",
      "t = 9.60, E = -0.499570577464443, delta = 0.00011474424535151138\n",
      "t = 9.65, E = -0.4995710731636199, delta = 0.00011252065929768624\n",
      "t = 9.70, E = -0.4995715496664836, delta = 0.00011034119507027872\n",
      "t = 9.75, E = -0.49957200772171695, delta = 0.00010820494401464745\n",
      "t = 9.80, E = -0.4995724480484432, delta = 0.00010611101728559804\n",
      "t = 9.85, E = -0.4995728713374146, delta = 0.00010405854537834887\n",
      "t = 9.90, E = -0.4995732782521129, delta = 0.00010204667767199222\n",
      "t = 9.95, E = -0.499573669429905, delta = 0.00010007458198482062\n",
      "t = 10.00, E = -0.49957404548303963, delta = 9.814144414135426e-05\n",
      "t = 10.05, E = -0.49957440699967576, delta = 9.624646755073531e-05\n",
      "t = 10.10, E = -0.49957475454482536, delta = 9.438887279601023e-05\n",
      "t = 10.15, E = -0.49957508866132905, delta = 9.256789723411461e-05\n",
      "t = 10.20, E = -0.49957540987067367, delta = 9.078279460622369e-05\n",
      "t = 10.25, E = -0.4995757186738966, delta = 8.903283465800812e-05\n",
      "t = 10.30, E = -0.49957601555238695, delta = 8.731730276973821e-05\n",
      "t = 10.35, E = -0.49957630096865036, delta = 8.563549959589147e-05\n",
      "t = 10.40, E = -0.49957657536708167, delta = 8.39867407138316e-05\n",
      "t = 10.45, E = -0.49957683917467677, delta = 8.237035628142655e-05\n",
      "t = 10.50, E = -0.4995770928017237, delta = 8.078569070341228e-05\n",
      "t = 10.55, E = -0.49957733664248166, delta = 7.923210230612774e-05\n",
      "t = 10.60, E = -0.4995775710757777, delta = 7.770896302040662e-05\n",
      "t = 10.65, E = -0.4995777964656521, delta = 7.621565807250041e-05\n",
      "t = 10.70, E = -0.4995780131619415, delta = 7.475158568251044e-05\n",
      "t = 10.75, E = -0.49957822150082404, delta = 7.331615677072886e-05\n",
      "t = 10.80, E = -0.4995784218053941, delta = 7.190879467081699e-05\n",
      "t = 10.85, E = -0.4995786143861193, delta = 7.052893485037346e-05\n",
      "t = 10.90, E = -0.4995787995413986, delta = 6.917602463837236e-05\n",
      "t = 10.95, E = -0.4995789775580074, delta = 6.784952295915618e-05\n",
      "t = 11.00, E = -0.4995791487115537, delta = 6.65489000731371e-05\n",
      "t = 11.05, E = -0.4995793132669465, delta = 6.527363732367751e-05\n",
      "t = 11.10, E = -0.4995794714787712, delta = 6.402322689027968e-05\n",
      "t = 11.15, E = -0.4995796235917735, delta = 6.27971715476871e-05\n",
      "t = 11.20, E = -0.4995797698411124, delta = 6.159498443083694e-05\n",
      "t = 11.25, E = -0.49957991045289263, delta = 6.0416188805529065e-05\n",
      "t = 11.30, E = -0.49958004564439623, delta = 5.926031784459476e-05\n",
      "t = 11.35, E = -0.4995801756245395, delta = 5.8126914409409464e-05\n",
      "t = 11.40, E = -0.49958030059406716, delta = 5.7015530836792135e-05\n",
      "t = 11.45, E = -0.4995804207459967, delta = 5.592572873076155e-05\n",
      "t = 11.50, E = -0.49958053626584653, delta = 5.485707875942027e-05\n",
      "t = 11.55, E = -0.4995806473319673, delta = 5.3809160456676326e-05\n",
      "t = 11.60, E = -0.49958075411578673, delta = 5.2781562028583465e-05\n",
      "t = 11.65, E = -0.4995808567821344, delta = 5.177388016425778e-05\n",
      "t = 11.70, E = -0.49958095548942927, delta = 5.07857198513446e-05\n",
      "t = 11.75, E = -0.49958105039001177, delta = 4.981669419575115e-05\n",
      "t = 11.80, E = -0.4995811416302832, delta = 4.88664242456378e-05\n",
      "t = 11.85, E = -0.49958122935102034, delta = 4.793453881946067e-05\n",
      "t = 11.90, E = -0.49958131368757386, delta = 4.702067433817562e-05\n",
      "t = 11.95, E = -0.4995813947700377, delta = 4.6124474661181814e-05\n",
      "t = 12.00, E = -0.49958147272352943, delta = 4.524559092617475e-05\n",
      "t = 12.05, E = -0.4995815476683085, delta = 4.4383681392749806e-05\n",
      "t = 12.10, E = -0.49958161971999554, delta = 4.353841128943977e-05\n",
      "t = 12.15, E = -0.49958168898979055, delta = 4.2709452664484666e-05\n",
      "t = 12.20, E = -0.4995817555845533, delta = 4.18964842399234e-05\n",
      "t = 12.25, E = -0.49958181960707815, delta = 4.109919126907152e-05\n",
      "t = 12.30, E = -0.49958188115617536, delta = 4.031726539719662e-05\n",
      "t = 12.35, E = -0.49958194032687103, delta = 3.955040452552742e-05\n",
      "t = 12.40, E = -0.49958199721049124, delta = 3.879831267815841e-05\n",
      "t = 12.45, E = -0.4995820518948901, delta = 3.806069987218006e-05\n",
      "t = 12.50, E = -0.49958210446455587, delta = 3.7337281990580014e-05\n",
      "t = 12.55, E = -0.49958215500066433, delta = 3.662778065815263e-05\n",
      "t = 12.60, E = -0.49958220358132643, delta = 3.5931923120228335e-05\n",
      "t = 12.65, E = -0.49958225028163455, delta = 3.524944212392086e-05\n",
      "t = 12.70, E = -0.4995822951737772, delta = 3.458007580232288e-05\n",
      "t = 12.75, E = -0.4995823383272058, delta = 3.392356756114743e-05\n",
      "t = 12.80, E = -0.49958237980867676, delta = 3.327966596796325e-05\n",
      "t = 12.85, E = -0.4995824196823835, delta = 3.264812464381456e-05\n",
      "t = 12.90, E = -0.49958245801007106, delta = 3.202870215739044e-05\n",
      "t = 12.95, E = -0.49958249485110845, delta = 3.142116192148483e-05\n",
      "t = 13.00, E = -0.4995825302625858, delta = 3.082527209174369e-05\n",
      "t = 13.05, E = -0.4995825642993856, delta = 3.0240805467697276e-05\n",
      "t = 13.10, E = -0.49958259701431784, delta = 2.966753939592355e-05\n",
      "t = 13.15, E = -0.49958262845814605, delta = 2.9105255675469267e-05\n",
      "t = 13.20, E = -0.49958265867967644, delta = 2.855374046515664e-05\n",
      "t = 13.25, E = -0.49958268772588005, delta = 2.8012784193218034e-05\n",
      "t = 13.30, E = -0.4995827156419069, delta = 2.7482181468587233e-05\n",
      "t = 13.35, E = -0.4995827424711556, delta = 2.6961730994411257e-05\n",
      "t = 13.40, E = -0.49958276825540116, delta = 2.6451235483279077e-05\n",
      "t = 13.45, E = -0.49958279303477643, delta = 2.595050157442433e-05\n",
      "t = 13.50, E = -0.49958281684791667, delta = 2.545933975264558e-05\n",
      "t = 13.55, E = -0.499582839731966, delta = 2.4977564269026667e-05\n",
      "t = 13.60, E = -0.49958286172263333, delta = 2.4504993063383883e-05\n",
      "t = 13.65, E = -0.49958288285424174, delta = 2.4041447688399772e-05\n",
      "t = 13.70, E = -0.4995829031598644, delta = 2.358675323538671e-05\n",
      "t = 13.75, E = -0.4995829226712333, delta = 2.3140738261622952e-05\n",
      "t = 13.80, E = -0.49958294141892484, delta = 2.2703234719407782e-05\n",
      "t = 13.85, E = -0.4995829594323108, delta = 2.2274077886407234e-05\n",
      "t = 13.90, E = -0.4995829767396493, delta = 2.185310629772332e-05\n",
      "t = 13.95, E = -0.49958299336811274, delta = 2.1440161679236424e-05\n",
      "t = 14.00, E = -0.4995830093438417, delta = 2.1035088882561648e-05\n",
      "t = 14.05, E = -0.49958302469197097, delta = 2.0637735821249484e-05\n",
      "t = 14.10, E = -0.49958303943665117, delta = 2.024795340842924e-05\n",
      "t = 14.15, E = -0.49958305360116584, delta = 1.986559549575837e-05\n",
      "t = 14.20, E = -0.4995830672078587, delta = 1.9490518813619027e-05\n",
      "t = 14.25, E = -0.4995830802782247, delta = 1.9122582912784e-05\n",
      "t = 14.30, E = -0.49958309283294944, delta = 1.8761650107053154e-05\n",
      "t = 14.35, E = -0.49958310489193136, delta = 1.8407585417280978e-05\n",
      "t = 14.40, E = -0.49958311647427694, delta = 1.8060256516600668e-05\n",
      "t = 14.45, E = -0.49958312759839707, delta = 1.7719533676645537e-05\n",
      "t = 14.50, E = -0.499583138281969, delta = 1.7385289715168716e-05\n",
      "t = 14.55, E = -0.49958314854203184, delta = 1.7057399944451024e-05\n",
      "t = 14.60, E = -0.4995831583949415, delta = 1.6735742121111764e-05\n",
      "t = 14.65, E = -0.4995831678564133, delta = 1.6420196396719473e-05\n",
      "t = 14.70, E = -0.4995831769416017, delta = 1.6110645269630066e-05\n",
      "t = 14.75, E = -0.4995831856650576, delta = 1.580697353769433e-05\n",
      "t = 14.80, E = -0.4995831940407799, delta = 1.5509068252064363e-05\n",
      "t = 14.85, E = -0.4995832020822278, delta = 1.5216818671896846e-05\n",
      "t = 14.90, E = -0.49958320980236076, delta = 1.493011621999667e-05\n",
      "t = 14.95, E = -0.4995832172136076, delta = 1.464885443951195e-05\n",
      "t = 15.00, E = -0.4995832243279603, delta = 1.4372928951349914e-05\n",
      "t = 15.05, E = -0.49958323115694225, delta = 1.4102237412609437e-05\n",
      "t = 15.10, E = -0.49958323771161317, delta = 1.3836679475831319e-05\n",
      "t = 15.15, E = -0.4995832440026206, delta = 1.3576156749081222e-05\n",
      "t = 15.20, E = -0.49958325004022847, delta = 1.3320572756925144e-05\n",
      "t = 15.25, E = -0.49958325583427077, delta = 1.3069832902173065e-05\n",
      "t = 15.30, E = -0.4995832613942104, delta = 1.2823844428329953e-05\n",
      "t = 15.35, E = -0.4995832667291803, delta = 1.2582516383024171e-05\n",
      "t = 15.40, E = -0.49958327184793827, delta = 1.2345759581989407e-05\n",
      "t = 15.45, E = -0.49958327675889763, delta = 1.2113486573910458e-05\n",
      "t = 15.50, E = -0.49958328147016684, delta = 1.1885611605946991e-05\n",
      "t = 15.55, E = -0.49958328598953944, delta = 1.1662050590005575e-05\n",
      "t = 15.60, E = -0.499583290324492, delta = 1.1442721069623893e-05\n",
      "t = 15.65, E = -0.49958329448225064, delta = 1.1227542187643938e-05\n",
      "t = 15.70, E = -0.49958329846973776, delta = 1.101643465448593e-05\n",
      "t = 15.75, E = -0.49958330229360914, delta = 1.0809320717021023e-05\n",
      "t = 15.80, E = -0.49958330596027223, delta = 1.0606124128272473e-05\n",
      "t = 15.85, E = -0.4995833094759018, delta = 1.0406770117442698e-05\n",
      "t = 15.90, E = -0.49958331284639246, delta = 1.0211185360856412e-05\n",
      "t = 15.95, E = -0.4995833160774682, delta = 1.0019297953297677e-05\n",
      "t = 16.00, E = -0.49958331917459686, delta = 9.831037379958307e-06\n",
      "t = 16.05, E = -0.49958332214304374, delta = 9.646334489089024e-06\n",
      "t = 16.10, E = -0.4995833249878684, delta = 9.465121464985512e-06\n",
      "t = 16.15, E = -0.4995833277139494, delta = 9.287331801794713e-06\n",
      "t = 16.20, E = -0.4995833303259669, delta = 9.11290027755276e-06\n",
      "t = 16.25, E = -0.4995833328284262, delta = 8.941762929007734e-06\n",
      "t = 16.30, E = -0.4995833352256391, delta = 8.773857026822098e-06\n",
      "t = 16.35, E = -0.49958333752178047, delta = 8.609121051325778e-06\n",
      "t = 16.40, E = -0.49958333972085445, delta = 8.447494668687072e-06\n",
      "t = 16.45, E = -0.49958334182668956, delta = 8.28891870766658e-06\n",
      "t = 16.50, E = -0.49958334384300884, delta = 8.133335136755288e-06\n",
      "t = 16.55, E = -0.49958334577334357, delta = 7.980687041856166e-06\n",
      "t = 16.60, E = -0.49958334762114565, delta = 7.830918604328183e-06\n",
      "t = 16.65, E = -0.49958334938965787, delta = 7.683975079576651e-06\n",
      "t = 16.70, E = -0.4995833510820776, delta = 7.539802775909297e-06\n",
      "t = 16.75, E = -0.4995833527014187, delta = 7.39834903403784e-06\n",
      "t = 16.80, E = -0.4995833542506143, delta = 7.259562206810659e-06\n",
      "t = 16.85, E = -0.49958335573247753, delta = 7.1233916394384625e-06\n",
      "t = 16.90, E = -0.4995833571496733, delta = 6.989787650081615e-06\n",
      "t = 16.95, E = -0.499583358504826, delta = 6.858701510827631e-06\n",
      "t = 17.00, E = -0.49958335980041313, delta = 6.730085429168813e-06\n",
      "t = 17.05, E = -0.4995833610388622, delta = 6.603892529622538e-06\n",
      "t = 17.10, E = -0.49958336222241884, delta = 6.480076835851439e-06\n",
      "t = 17.15, E = -0.49958336335335685, delta = 6.358593253240042e-06\n",
      "t = 17.20, E = -0.4995833644337638, delta = 6.2393975515957835e-06\n",
      "t = 17.25, E = -0.4995833654657044, delta = 6.122446348395461e-06\n",
      "t = 17.30, E = -0.49958336645113943, delta = 6.007697092213425e-06\n",
      "t = 17.35, E = -0.49958336739194514, delta = 5.8951080465875955e-06\n",
      "t = 17.40, E = -0.4995833682899593, delta = 5.784638274189898e-06\n",
      "t = 17.45, E = -0.4995833691469189, delta = 5.676247621270423e-06\n",
      "t = 17.50, E = -0.499583369964491, delta = 5.56989670237199e-06\n",
      "t = 17.55, E = -0.49958337074429743, delta = 5.465546885485922e-06\n",
      "t = 17.60, E = -0.49958337148788345, delta = 5.363160277455025e-06\n",
      "t = 17.65, E = -0.4995833721967417, delta = 5.262699709484765e-06\n",
      "t = 17.70, E = -0.4995833728722978, delta = 5.1641287232535275e-06\n",
      "t = 17.75, E = -0.4995833735159138, delta = 5.067411557064257e-06\n",
      "t = 17.80, E = -0.499583374128921, delta = 4.972513132346801e-06\n",
      "t = 17.85, E = -0.49958337471259073, delta = 4.879399040406547e-06\n",
      "t = 17.90, E = -0.49958337526813257, delta = 4.788035529509749e-06\n",
      "t = 17.95, E = -0.4995833757967243, delta = 4.698389492123851e-06\n",
      "t = 18.00, E = -0.49958337629947974, delta = 4.61042845245276e-06\n",
      "t = 18.05, E = -0.49958337677749154, delta = 4.524120554255205e-06\n",
      "t = 18.10, E = -0.4995833772317929, delta = 4.439434548774755e-06\n",
      "t = 18.15, E = -0.4995833776633675, delta = 4.3563397831903835e-06\n",
      "t = 18.20, E = -0.49958337807319864, delta = 4.27480618886814e-06\n",
      "t = 18.25, E = -0.4995833784621674, delta = 4.19480427024057e-06\n",
      "t = 18.30, E = -0.4995833788312018, delta = 4.116305093725957e-06\n",
      "t = 18.35, E = -0.4995833791811294, delta = 4.0392802767953475e-06\n",
      "t = 18.40, E = -0.4995833795127466, delta = 3.9637019774718276e-06\n",
      "t = 18.45, E = -0.49958337982686535, delta = 3.889542883816049e-06\n",
      "t = 18.50, E = -0.4995833801242088, delta = 3.816776203732156e-06\n",
      "t = 18.55, E = -0.499583380405513, delta = 3.74537565498456e-06\n",
      "t = 18.60, E = -0.4995833806714753, delta = 3.67531545533949e-06\n",
      "t = 18.65, E = -0.4995833809227566, delta = 3.60657031296608e-06\n",
      "t = 18.70, E = -0.49958338115998036, delta = 3.5391154169564875e-06\n",
      "t = 18.75, E = -0.49958338138377256, delta = 3.4729264281546047e-06\n",
      "t = 18.80, E = -0.4995833815947349, delta = 3.407979469969589e-06\n",
      "t = 18.85, E = -0.4995833817933935, delta = 3.3442511196104276e-06\n",
      "t = 18.90, E = -0.4995833819803255, delta = 3.2817183992745244e-06\n",
      "t = 18.95, E = -0.4995833821560296, delta = 3.2203587675907133e-06\n",
      "t = 19.00, E = -0.4995833823210178, delta = 3.160150111302432e-06\n",
      "t = 19.05, E = -0.4995833824757576, delta = 3.1010707370263528e-06\n",
      "t = 19.10, E = -0.4995833826207177, delta = 3.043099363175595e-06\n",
      "t = 19.15, E = -0.49958338275633607, delta = 2.986215112035833e-06\n",
      "t = 19.20, E = -0.49958338288303417, delta = 2.9303975021056762e-06\n",
      "t = 19.25, E = -0.4995833830012242, delta = 2.8756264403687014e-06\n",
      "t = 19.30, E = -0.4995833831112826, delta = 2.8218822149901244e-06\n",
      "t = 19.35, E = -0.49958338321359114, delta = 2.769145487896495e-06\n",
      "t = 19.40, E = -0.49958338330850505, delta = 2.7173972876961514e-06\n",
      "t = 19.45, E = -0.49958338339636765, delta = 2.666619002579077e-06\n",
      "t = 19.50, E = -0.49958338347751186, delta = 2.6167923735025026e-06\n",
      "t = 19.55, E = -0.4995833835522557, delta = 2.567899487432638e-06\n",
      "t = 19.60, E = -0.49958338362088067, delta = 2.5199227706213914e-06\n",
      "t = 19.65, E = -0.49958338368369326, delta = 2.472844982306396e-06\n",
      "t = 19.70, E = -0.49958338374099, delta = 2.4266492081314853e-06\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ic| 'Time per time step:'\n",
      "ic| delta_t/len(t_range): 0.012095324993133545\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t = 19.75, E = -0.49958338379300954, delta = 2.3813188540880406e-06\n",
      "t = 19.80, E = -0.49958338384001627, delta = 2.3368376402466894e-06\n",
      "t = 19.85, E = -0.49958338388225104, delta = 2.2931895948679446e-06\n",
      "t = 19.90, E = -0.49958338391996143, delta = 2.2503590484403303e-06\n",
      "t = 19.95, E = -0.4995833839533625, delta = 2.20833062796068e-06\n"
     ]
    }
   ],
   "source": [
    "# Set up an initial ground state guess, here a gaussian.\n",
    "# Note that we are working with a reduced wavefunction.\n",
    "psi = rr**.5 *  np.exp(-0.5*(rr**2 + zz**2))\n",
    "psi /= np.linalg.norm(psi)\n",
    "\n",
    "get_energy = lambda psi: np.sum(psi.conj() * solver.apply_hamiltonian(psi)).real\n",
    "E = get_energy(psi)\n",
    "ic('Initial guess energy:')\n",
    "ic(E)\n",
    "\n",
    "# Abuse the propagator to calculate the ground state using imaginary time propagation.\n",
    "dt = -0.05j\n",
    "\n",
    "# Set up the splitting scheme for the Crank-Nicolson propagator.\n",
    "solver.setup_splitting_scheme(dt)\n",
    "\n",
    "t_range = np.arange(0, 20, np.abs(dt))\n",
    "start = time()\n",
    "for k,t in enumerate(t_range):\n",
    "    psi_new = solver.propagate_crank_nicolson(psi, t)\n",
    "\n",
    "    psi_new /= np.linalg.norm(psi_new)\n",
    "    delta = np.linalg.norm(psi - psi_new)\n",
    "    \n",
    "    psi = psi_new\n",
    "    E = get_energy(psi)\n",
    "    print(f't = {t:.2f}, E = {E}, delta = {delta}')\n",
    "    \n",
    "delta_t = time()-start\n",
    "ic('Time per time step:')\n",
    "ic(delta_t/len(t_range))\n",
    "psi_ground_state = psi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ic| self.m_i: array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7., -8., -7., -6., -5., -4.,\n",
      "                     -3., -2., -1.])\n",
      "ic| cylinder_fdm_3d.py:169 in __init__() at 22:04:09.302\n",
      "ic| self.n_r: 300, self.n_z: 300, self.n_m: 16\n",
      "ic| self.r_max: 20, self.z_max: 20\n",
      "ic| self.n_dof: 1440000\n",
      "ic| m: 0.0, N: 41.56142818209105\n",
      "ic| V_FFT[i, :, :]: array([[-0.05033417+0.j, -0.05067311+0.j, -0.05101666+0.j, ...,\n",
      "                            -0.05101666+0.j, -0.05067311+0.j, -0.05033417+0.j],\n",
      "                           [-0.05033332+0.j, -0.05067225+0.j, -0.05101578+0.j, ...,\n",
      "                            -0.05101578+0.j, -0.05067225+0.j, -0.05033332+0.j],\n",
      "                           [-0.05033191+0.j, -0.05067082+0.j, -0.05101431+0.j, ...,\n",
      "                            -0.05101431+0.j, -0.05067082+0.j, -0.05033191+0.j],\n",
      "                           ...,\n",
      "                           [-0.0356514 +0.j, -0.03577123+0.j, -0.03589147+0.j, ...,\n",
      "                            -0.03589147+0.j, -0.03577123+0.j, -0.0356514 +0.j],\n",
      "                           [-0.03559183+0.j, -0.03571106+0.j, -0.03583069+0.j, ...,\n",
      "                            -0.03583069+0.j, -0.03571106+0.j, -0.03559183+0.j],\n",
      "                           [-0.03553236+0.j, -0.035651  +0.j, -0.03577002+0.j, ...,\n",
      "                            -0.03577002+0.j, -0.035651  +0.j, -0.03553236+0.j]])\n",
      "ic| m: 1.0, N: 0.0\n",
      "ic| V_FFT[i, :, :]: array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           ...,\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])\n",
      "ic| m: 2.0, N: 0.0\n",
      "ic| V_FFT[i, :, :]: array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           ...,\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])\n",
      "ic| m: 3.0, N: 0.0\n",
      "ic| V_FFT[i, :, :]: array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           ...,\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])\n",
      "ic| m: 4.0, N: 0.0\n",
      "ic| V_FFT[i, :, :]: array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           ...,\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])\n",
      "ic| m: 5.0, N: 0.0\n",
      "ic| V_FFT[i, :, :]: array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           ...,\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])\n",
      "ic| m: 6.0, N: 0.0\n",
      "ic| V_FFT[i, :, :]: array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           ...,\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])\n",
      "ic| m: 7.0, N: 0.0\n",
      "ic| V_FFT[i, :, :]: array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           ...,\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])\n",
      "ic| m: -8.0, N: 0.0\n",
      "ic| V_FFT[i, :, :]: array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           ...,\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])\n",
      "ic| m: -7.0, N: 0.0\n",
      "ic| V_FFT[i, :, :]: array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           ...,\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])\n",
      "ic| m: -6.0, N: 0.0\n",
      "ic| V_FFT[i, :, :]: array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           ...,\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])\n",
      "ic| m: -5.0, N: 0.0\n",
      "ic| V_FFT[i, :, :]: array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           ...,\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])\n",
      "ic| m: -4.0, N: 0.0\n",
      "ic| V_FFT[i, :, :]: array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           ...,\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])\n",
      "ic| m: -3.0, N: 0.0\n",
      "ic| V_FFT[i, :, :]: array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           ...,\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])\n",
      "ic| m: -2.0, N: 0.0\n",
      "ic| V_FFT[i, :, :]: array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           ...,\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])\n",
      "ic| m: -1.0, N: 0.0\n",
      "ic| V_FFT[i, :, :]: array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           ...,\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],\n",
      "                           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])\n",
      "ic| m_max: 0\n",
      "ic| m: 0.0\n",
      "ic| m_list[i]: 0.0, np.linalg.norm(V_m[i]): 41.56142818209105\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAHFCAYAAAAaD0bAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABaAklEQVR4nO3daXiTZd428DNLk3Tfd7qytVDWFpAdRYuAIiM+4or7TMdxFBjmQcEZFV+tOj6Owygwo+AyjMooqOhUoexbZSktSymlQFdo6b63SZPc74c0wdoC3a+kOX/H0Q/cvZr+cwPJmWuVSZIkgYiIiMiOyEUXQERERNTXGICIiIjI7jAAERERkd1hACIiIiK7wwBEREREdocBiIiIiOwOAxARERHZHQYgIiIisjsMQERERGR3GICIqMs+/vhjyGQyHDt2THQpVunRRx9FeHi46DKIqB0MQERERGR3GICIyC5IkoTGxkbRZRCRlWAAIqJe1dTUhD/84Q8YPXo03N3d4eXlhYkTJ+Lbb79t0/bLL7/EhAkT4O7uDicnJ0RGRuLxxx9v1aampgbLli1DREQEVCoVgoODsXjxYtTX17dqJ5PJ8Mwzz2DdunWIjo6GWq3GJ598cs06w8PDcccdd+Drr7/GyJEjodFoEBkZidWrV7dqZx72y83NbXV9z549kMlk2LNnz3XvR08+RyLqOqXoAoiof9NqtaioqMCyZcsQHBwMnU6HHTt24O6778ZHH32ERYsWAQBSUlKwcOFCLFy4EC+//DI0Gg3y8vKwa9cuy2M1NDRg+vTpKCwsxIoVKzBy5EhkZGTgz3/+M06dOoUdO3ZAJpNZ2n/zzTfYv38//vznPyMgIAB+fn7XrTU9PR2LFy/Gyy+/jICAAPz73//Gc889B51Oh2XLlnX7XvTGcySiLpKIiLroo48+kgBIR48e7fDP6PV6qbm5WXriiSekMWPGWK6//fbbEgCpqqrqmj+bmJgoyeXyNr/vq6++kgBISUlJlmsAJHd3d6mioqJDdYWFhUkymUxKT09vdf22226T3NzcpPr6ekmSrj7nnJycVu12794tAZB2795tufbII49IYWFhvfYciajrOARGRL3uyy+/xOTJk+Hi4gKlUgkHBwesX78emZmZljbjxo0DANx77734z3/+g0uXLrV5nO+//x4xMTEYPXo09Hq95WvWrFntDj/dcsst8PT07HCdw4cPx6hRo1pde+CBB1BTU4Pjx4934hm3rzeeIxF1DQMQEfWqLVu24N5770VwcDA2btyIlJQUHD16FI8//jiampos7aZNm4ZvvvkGer0eixYtwoABAxATE4PPP//c0ubKlSs4efIkHBwcWn25urpCkiSUlZW1+t2BgYGdqjUgIOCa18rLyzv1WO3pjedIRF3DOUBE1Ks2btyIiIgIbNq0qdXcFa1W26btXXfdhbvuugtarRY//fQTEhMT8cADDyA8PBwTJ06Ej48PHB0dsWHDhnZ/l4+PT6s/d3auTHFx8TWveXt7AwA0Gk279Xc0mPT0cySirmEAIqJeJZPJoFKpWoWR4uLidleBmanVakyfPh0eHh7Ytm0b0tLSMHHiRNxxxx14/fXX4e3tjYiIiB6vNSMjAydOnGg1DPbZZ5/B1dUVY8eOBQDLxoYnT57E0KFDLe22bt3aqd8l6jkSkQkDEBF1265du9osCweAOXPm4I477sCWLVvw9NNP45577kFBQQFeffVVBAYGIjs729L2z3/+MwoLCzFz5kwMGDAAVVVV+Nvf/gYHBwdMnz4dALB48WJs3rwZ06ZNw5IlSzBy5EgYjUbk5+dj+/bt+MMf/oAJEyZ0+XkEBQVh3rx5ePnllxEYGIiNGzciOTkZb775JpycnACY5vEMHToUy5Ytg16vh6enJ77++mscOHDgho9vDc+RiEwYgIio25YvX97u9ZycHDz22GMoKSnBunXrsGHDBkRGRuL5559HYWEhXnnlFUvbCRMm4NixY1i+fDlKS0vh4eGBuLg47Nq1C8OHDwcAODs7Y//+/XjjjTfwz3/+Ezk5OXB0dERoaChuvfXWbh87MXr0aDz22GN46aWXkJ2djaCgILzzzjtYsmSJpY1CocB3332HZ555BgkJCVCr1bjvvvvw3nvvYe7cudd9fGt4jkRkIpMkSRJdBBGRaOHh4YiJicH3338vuhQi6gNcBUZERER2hwGIiIiI7A6HwIiIiMjusAeIiIiI7A4DEBEREdkdBiAiIiKyO9wHqB1GoxGXL1+Gq6trp7fSJyIiIjEkSUJtbS2CgoIgl1+/j4cBqB2XL19GSEiI6DKIiIioCwoKCjBgwIDrtmEAaoerqysA0w10c3MTXA0RERF1RE1NDUJCQizv49fDANQO87CXm5sbAxAREZGN6cj0FU6CJiIiIrvDAERERER2hwGIiIiI7A4DEBEREdkdBiAiIiKyOwxAREREZHcYgIiIiMjuMAARERGR3WEAIiIiIrvDAERERER2R3gAWrNmDSIiIqDRaBAbG4v9+/dft/3evXsRGxsLjUaDyMhIrFu3rtX3P/74Y8hksjZfTU1Nvfk0iIiIyIYIDUCbNm3C4sWLsXLlSqSlpWHq1KmYPXs28vPz222fk5ODOXPmYOrUqUhLS8OKFSvw7LPPYvPmza3aubm5oaioqNWXRqPpi6dERERENkAmSZIk6pdPmDABY8eOxdq1ay3XoqOjMX/+fCQmJrZpv3z5cmzduhWZmZmWawkJCThx4gRSUlIAmHqAFi9ejKqqqi7XVVNTA3d3d1RXV/MwVKJ+xGiUUNmgg8EowUmthIua50ET9Sedef8W9r9fp9MhNTUVzz//fKvr8fHxOHToULs/k5KSgvj4+FbXZs2ahfXr16O5uRkODg4AgLq6OoSFhcFgMGD06NF49dVXMWbMmGvWotVqodVqLX+uqanp6tMiIitiNEr4KaccP5wqxtHcClworUOz4epnPm9nFYYFuWH6EF/MHhGIYA9HgdUSUV8SFoDKyspgMBjg7+/f6rq/vz+Ki4vb/Zni4uJ22+v1epSVlSEwMBBRUVH4+OOPMWLECNTU1OBvf/sbJk+ejBMnTmDw4MHtPm5iYiJeeeWVnnliRCScTm/El6kFWLvnAgorG9t8XyYDJAkor9dhf3YZ9meX4fWkTNweE4Df3TwIw4PcBVRNRH1JeP+vTCZr9WdJktpcu1H7n1+/6aabcNNNN1m+P3nyZIwdOxZ///vfsXr16nYf84UXXsDSpUstf66pqUFISEjnnggRWYVDF8qwYssp5JY3AABc1UrMHRmIaUN8MSLYHf5uGjgoZGjQGXChtA5HcyuxLaMYR3IqkHSqGD+cLsbCuBA8PzsKHk4qwc+GiHqLsADk4+MDhULRprenpKSkTS+PWUBAQLvtlUolvL292/0ZuVyOcePGITs7+5q1qNVqqNXqTj4DIrImWr0Br/03E5+m5AEAfFzUeHrGQDwwIRQaB0Wb9s5qJUYO8MDIAR54YkoEsopr8fdd2fj+ZBG+OFqAvedK8c69ozFxYPuvLURk24StAlOpVIiNjUVycnKr68nJyZg0aVK7PzNx4sQ27bdv3464uDjL/J9fkiQJ6enpCAwM7JnCicjqlNQ04f5//mQJPw9OCMXuZdPx+JSIdsNPe4YGuOK9B8biq4SJiPBxRlF1Ex748CesP5ADgWtFiKiXCF0Gv3TpUnz44YfYsGEDMjMzsWTJEuTn5yMhIQGAaWhq0aJFlvYJCQnIy8vD0qVLkZmZiQ0bNmD9+vVYtmyZpc0rr7yCbdu24eLFi0hPT8cTTzyB9PR0y2MSUf9ysbQO898/iOP5VXDTKPHRY+Pw2q9GwFXT/oeiG4kL98L3v5+Ce2IHQJKAV78/g5e2ZsBgZAgi6k+EzgFauHAhysvLsWrVKhQVFSEmJgZJSUkICwsDABQVFbXaEygiIgJJSUlYsmQJ3n//fQQFBWH16tVYsGCBpU1VVRV+/etfo7i4GO7u7hgzZgz27duH8ePH9/nzI6LelVVciwc/PIyyOi0ifZ2x/pFxiPBx7vbjOquV+Ms9IzHE3wWJP5zFpyl5aNQZ8OaCkZDLrz1HkYhsh9B9gKwV9wEisn45ZfW4Z+0hlNfrEB3ohn89MR4+Lj0/l2/rictY/EUajBKwMC4EbywYcd2FGkQkTmfev4UfhUFE1FmltVos2nAY5fU6DA9yw+dPTeiV8AMA80YF4d37xkAuAzYdK8A7yed65fcQUd9iACIim9KoM+Cxj4+goKIRoV5O+Pix8b2+XH3eqCC8sWAkAODvu87jy2MFvfr7iKj3MQARkc2QJAkrvj6F05dq4OWswiePj4eva99sYXFvXAieuXkQAOCFLaeQmlfRJ7+XiHoHAxAR2YyNP+Xh67RLUMhleP+BsT0y4bkzlt42BHNHBEJvlPDMZ2moqNf16e8nop7DAERENuFEQRVWfX8GAPD87VFCNiiUy2V4856RiPQ17RO0eFM6jFweT2STGICIyOo16gxYsikdzQYJs2MC8OTUCGG1uKiVWPPgWGgc5Nh3rhQfHrgorBYi6joGICKyem/8kImLZfXwd1Mj8W7xy9CjAtzw8p3DAQBvbz+H7Cu1Qushos5jACIiq7Y/uxSftBxx8Zd7RlnNAaULx4Xg5qG+0OmNWPblCegNRtElEVEnMAARkdVq0Onx/OZTAIBFE8MwbYiv4IqukslkSLx7JNw0SpworMY/9nEojMiWMAARkdVavfM8LlU1ItjDEc/PjhJdThsB7hq8PM80FLZ6ZzYKKhoEV0REHcUARERWKau4Fh/uN/WqvDJvOJxUQo8uvKZfjQnGpIHe0OqNeGlrBk+OJ7IRDEBEZHWMRgkvfnMKeqOE24b549Zh/qJLuiaZTIZVd8XAQSHDrrMlSD5zRXRJRNQBDEBEZHW2nriMo7mVcHRQWIaYrNkgPxc8OTUSAPDKd2fQqDMIroiIboQBiIisSlOzAX/ZlgUAeOaWQQj2cBRcUcf8vqXWS1WN2HAwR3Q5RHQDDEBEZFU+OpiLS1WNCHTX4Ikp4jY87CwnlRJ/nDUUALB2zwWU12kFV0RE18MARERWo6JehzW7zwMA/jhrKDQOCsEVdc68UUGICXZDnVaPv+86L7ocIroOBiAishqrd2ajVqvH8CA3zB8dLLqcTpPLZVgxOxqA6eDWnLJ6wRUR0bUwABGRVSioaMC/D5t2fF45JxpyudjjLrpq0iAfzBjqC71Rwtstc5mIyPowABGRVXh/93k0GyRMHeyDSYN8RJfTLctvN23a+N9TRcgq5jlhRNaIAYiIhCuoaMBXqYUAgMW3DhZcTfdFB7phzogAAKZhPSKyPgxARCTc33dlQ2809f7EhnmJLqdHPDvTFOT+e6oIZ4trBFdDRL/EAEREQuWV12Pz8UsAgCW3DRFcTc+JCnDD3BGBANgLRGSNGICISKj3dp2HwShh+hBfjA31FF1OjzL3AiWdKkZmEXuBiKwJAxARCVNU3Yiv00y9P8/1g7k/vzQ0wNXSC7Ru7wXB1RDRzzEAEZEwGw7kQG+UMCHCq9/1/pj9dsZAAMD3J4tQWNkguBoiMmMAIiIhqhub8dnhfABAwvSBgqvpPTHB7pgyyAcGo4T1B3hGGJG1YAAiIiH+fTgP9ToDhvq7YsZQX9Hl9KpfTzOdFP/FkQJUNegEV0NEAAMQEQnQ1GzARwdzAZjCgUxmm7s+d9TUwT6IDnRDY7MBG3/KE10OEYEBiIgE+CbtEkprtQh01+DOUUGiy+l1MpkMCdNNvUAfH8pFU7NBcEVExABERH1Kkq7OhXl8cgRUSvt4GZozIhDBHo4oq9Phm5aVb0Qkjn288hCR1Th0oRzZJXVwUimwcHyI6HL6jINCjkcmhQEw9QJJkiS4IiL7xgBERH3q40O5AIAFYwfATeMgtpg+dm9cCDQOcpwtrsWRnArR5RDZNQYgIuozBRUN2Jl5BQAsvSH2xMNJhV+NCQYAfJKSK7YYIjvHAEREfWbj4TwYJWDKIB8M8nMVXY4QiyaGAwC2ZVxBUXWj2GKI7BgDEBH1iaZmAzYdLQAALJpof70/ZtGBbhgf4QWDUcK/f8oXXQ6R3WIAIqI+8W36JVQ1NGOApyNmRvuLLkeoRyeFAwA+P5LPJfFEgjAAEVGvkyQJnxwybQD48E1hUMj798aHNxI/zB+B7hqU1+uQdKpIdDlEdokBiIh63cnCapwpqoFKKce9cfaz9P1alAo57h8fCsB0PAYR9T0GICLqdV8cNc11mRMTAE9nleBqrMO9cSGQy4AjuRU4X1Inuhwiu8MARES9ql6rx9b0ywCA+1p6PQgIcNfglig/AMCmo5wMTdTXGICIqFd9f/Iy6nUGRPg4Y0KEl+hyrMp940yB8KvUQmj1nAxN1JcYgIioV33eMsdl4biQfn/qe2fNGOoLfzc1KhuasT3jiuhyiOwKAxAR9ZqzxTVIL6iCUi7DgrEDRJdjdZSKq5PCv+AwGFGfYgAiol5jXuF02zB/+LqqBVdjne6NC4FMBhw8X4688nrR5RDZDQYgIuoVTc0GfJ12CQAnP19PiJcTpg72BQDLTtlE1PsYgIioV2zLKEZ1YzOCPRwxdZCP6HKs2n3jTMNgX6ddgsEoCa6GyD4wABFRr9h83NT7syB2AOR2vvPzjcyM9oO7owOKqpuQcqFcdDlEdoEBiIh63JWaJhzILgUALBgbLLga66dWKnDnqEAAwJbjhYKrIbIPDEBE1OO+Tb8EowTEhXkizNtZdDk24e6WVXI/nC5GnVYvuBqi/o8BiIh6lCRJ2JxqGv66m0vfO2xMiAcifZzR2GzAj6eLRZdD1O8xABFRjzpTVIOsK7VQKeWYOyJQdDk2QyaT4e6W4cLNqRwGI+ptwgPQmjVrEBERAY1Gg9jYWOzfv/+67ffu3YvY2FhoNBpERkZi3bp112z7xRdfQCaTYf78+T1cNRFdy5aWyc+3RfvD3clBcDW25VctPWYpF8tRWNkguBqi/k1oANq0aRMWL16MlStXIi0tDVOnTsXs2bORn9/+jqg5OTmYM2cOpk6dirS0NKxYsQLPPvssNm/e3KZtXl4eli1bhqlTp/b20yCiFnqDEd+mm4e/OPm5s4I9HDEx0hsA8E3LHkpE1DuEBqB33nkHTzzxBJ588klER0fj3XffRUhICNauXdtu+3Xr1iE0NBTvvvsuoqOj8eSTT+Lxxx/H22+/3aqdwWDAgw8+iFdeeQWRkZF98VSICMD+7DKU1eng7azCtCG+osuxSQtiTb1Am49fgiRxTyCi3iIsAOl0OqSmpiI+Pr7V9fj4eBw6dKjdn0lJSWnTftasWTh27Biam5st11atWgVfX1888cQTHapFq9Wipqam1RcRdd7mliXc80YHwUEhfITdJs2OCYCjgwI5ZfVIL6gSXQ5RvyXsFaqsrAwGgwH+/v6trvv7+6O4uP0VEMXFxe221+v1KCsrAwAcPHgQ69evxwcffNDhWhITE+Hu7m75CgkJ6eSzIaLapmYknzGdaH73GK7+6ipntRLxw02vc1tPXBZcDVH/JfwjmkzWeodYSZLaXLtRe/P12tpaPPTQQ/jggw/g49PxrfdfeOEFVFdXW74KCngeD1Fnbc+4Aq3eiIG+zogJdhNdjk27c2QQAOD7k0U8GoOolyhF/WIfHx8oFIo2vT0lJSVtennMAgIC2m2vVCrh7e2NjIwM5Obm4s4777R832g0AgCUSiWysrIwcODANo+rVquhVvOkaqLu+O6kqbdi3qjg636IoRubNsQX7o4OKK3V4vDFckziWWpEPU5YD5BKpUJsbCySk5NbXU9OTsakSZPa/ZmJEye2ab99+3bExcXBwcEBUVFROHXqFNLT0y1f8+bNw80334z09HQObRH1kop6HQ5km4ah7xjFvX+6S6WUY3ZMAAAOgxH1FqFDYEuXLsWHH36IDRs2IDMzE0uWLEF+fj4SEhIAmIamFi1aZGmfkJCAvLw8LF26FJmZmdiwYQPWr1+PZcuWAQA0Gg1iYmJafXl4eMDV1RUxMTFQqVRCnidRf/fj6WLojRKGB7lhoK+L6HL6hXmjTMNgP5wuhk5vFFwNUf8jbAgMABYuXIjy8nKsWrUKRUVFiImJQVJSEsLCwgAARUVFrfYEioiIQFJSEpYsWYL3338fQUFBWL16NRYsWCDqKRARgO9aeinubHnTpu6bEOkNP1c1Smq12J9dipnR7U8NIKKukUncaKKNmpoauLu7o7q6Gm5unMxJdD0lNU2YkLgTkgQcWH4zBng6iS6p33jluwx8dDAXd40Owt/uGyO6HCKr15n3b+GrwIjItn1/sgiSBMSGeTL89DDzMFjymSto1BkEV0PUvzAAEVG3mFd/3TmSk5972ugQD4R6OaFBZ8COzCuiyyHqVxiAiKjLCioakJZfBbkMmMMA1ONkMhnubFlVx9VgRD2LAYiIuuz7k0UAgJsiveHnqhFcTf80b5TpUNm9WaWoaWq+QWsi6igGICLqMvPqr3lc/dVrhga4YpCfC3QGI3ZyGIyoxzAAEVGXXCytw5miGijlMtzesmkf9Y45Lfc36VT75yQSUecxABFRl/xw2vRmPHmQDzycuMlobzLPr9p7rhR1Wr3gaoj6BwYgIuqSH06b5v/MZu9Prxvq74pIH2fo9EbsOlsiuhyifoEBiIg6raCiAacv1UAuA24bxh2Ke5tMJsPsEaag+cOpIsHVEPUPDEBE1Gk/tgx/TYjwhreLWnA19mF2jGkYbHdWCRp0HAYj6i4GICLqtKSW4a85Izj81VeGB7kh1MsJTc1G7MkqFV0Okc1jACKiTimqbkRafhVkMmDWcAagvvLzYbAkDoMRdRsDEBF1yraW4a/YUE/4uXHzw740p2UYbNfZEjQ182wwou5gACKiTklqCUDc+6fvjRzgjmAPRzToDBwGI+omBiAi6rDSWi2O5lYAYAASQSaTWbYdMG9DQERdwwBERB22/UwxJMnUEzHA00l0OXZp9gjTMNjOTA6DEXUHAxARddiPHP4SbkyIBwLcNKjT6nHwfJnocohsFgMQEXVIZb0Ohy6UA7i6Jw31Pblchvjhps0nk8/wcFSirmIAIqIOSc68AoNRQlSAKyJ8nEWXY9fih5l64Ha0/J0QUecxABFRh5iXv7P3R7wJkV5w1ShRVqdDWn6l6HKIbBIDEBHdUINOjwMt801mxfDsL9EcFHLMjPIDAGznMBhRlzAAEdEN7TtXBq3eiBAvRwz1dxVdDgGIb9mFe1tGMSSJw2BEncUAREQ3ZJ5se1t0AGQymeBqCACmDfGFSilHXnkDskvqRJdDZHMYgIjouvQGI3adbQlAwzj8ZS1c1EpMGeQDANieUSy4GiLbwwBERNeVmleJyoZmeDg5YFy4p+hy6GfiWwIp5wERdR4DEBFdl3n465ahflAq+JJhTWZG+0MmA04WVuNyVaPocohsCl/NiOiaJElCciaHv6yVr6sasaGmXrkdmewFIuoMBiAiuqbskjrklTdApZRj2hBf0eVQO8y7Qm/PYAAi6gwGICK6JvPw1+SB3nBWKwVXQ+25rWVX6J8ulqO6oVlwNUS2gwGIiK7JPLnW/CZL1ifCxxlD/F2gN0rYnVUiuhwim8EARETtulLThBMFVQCAW6P9xBZD12U+G2z7GS6HJ+ooBiAiapd5Uu3oEA/4uWkEV0PXY54HtCerFFq9QXA1RLaBAYiI2mXZ/Zmrv6xeTJA7/FzVaNAZ8NPFCtHlENkEBiAiaqNOq8eh8+UArm62R9ZLLpdhZssw5U4uhyfqEAYgImpj37lS6AxGhHs7YZCfi+hyqANmRpmC6s7MEh6OStQBDEBE1MbPh794+KltmDzIB2qlHJeqGpF1pVZ0OURWjwGIiFoxHX5qWk7N5e+2w1GlwOSWw1F3ZnI5PNGNMAARUStpBVWobjQdfhobxsNPbQnnARF1HAMQEbWyu6X3Z9pgXyjkHP6yJeZ5QGkFVSir0wquhsi6MQARUSu7s0oBALdEcfNDWxPgrkFMsBsk6WqQJaL2MQARkUVxdRMyi2ogk4GHn9qoW1p6gXYxABFdFwMQEVmYz5IaHeIBL2eV4GqoK8zHluw7x12hia6HAYiILMzDJjcP5fCXrTLvCl2vM+Awd4UmuiYGICICAGj1Bhw8XwaAAciWyeUyy/wtDoMRXRsDEBEBAI7lVqJeZ4CvqxrDg9xEl0PdMDPaNA9oR+YV7gpNdA0MQEQE4Orw14whvpBz+btNm9KyK3RhZSPOXakTXQ6RVWIAIiIAwK6WCdA3c/m7zWu1K/RZbopI1B4GICJCXnk9LpbWQymXYcpgH9HlUA8wzwPisRhE7WMAIiLsadn8MC7cE24aB8HVUE8wH4txPL8S5dwVmqgNBiAisuz/w9Vf/UeguyOGB7XsCt0ScInoKgYgIjvXqDMg5UI5AM7/6W9mtvx9mgMuEV0lPACtWbMGERER0Gg0iI2Nxf79+6/bfu/evYiNjYVGo0FkZCTWrVvX6vtbtmxBXFwcPDw84OzsjNGjR+Nf//pXbz4FIpuWcrEMWr0RwR6OGOznIroc6kHmQLvvXCn0BqPgaoisi9AAtGnTJixevBgrV65EWloapk6ditmzZyM/P7/d9jk5OZgzZw6mTp2KtLQ0rFixAs8++yw2b95saePl5YWVK1ciJSUFJ0+exGOPPYbHHnsM27Zt66unRWRTdp81DY/MGOoLmYzL3/uTkQNMR5rUNumRmlcpuhwiqyKTBO6SNWHCBIwdOxZr1661XIuOjsb8+fORmJjYpv3y5cuxdetWZGZmWq4lJCTgxIkTSElJuebvGTt2LObOnYtXX321Q3XV1NTA3d0d1dXVcHPjhnDUf0mShKlv7UZhZSPWPxJn2UCP+o8lm9LxddolJEwfiOdnR4kuh6hXdeb9W1gPkE6nQ2pqKuLj41tdj4+Px6FDh9r9mZSUlDbtZ82ahWPHjqG5ublNe0mSsHPnTmRlZWHatGnXrEWr1aKmpqbVF5E9uFBah8LKRqiUckwc6C26HOoFM4b6AgD2cB4QUSvCAlBZWRkMBgP8/Vt/4vT390dxcXG7P1NcXNxue71ej7KyMsu16upquLi4QKVSYe7cufj73/+O22677Zq1JCYmwt3d3fIVEhLSjWdGZDvMZ0XdFOkNJ5VScDXUG6YP8YVcBpwtrsXlqkbR5RBZDeGToH8550CSpOvOQ2iv/S+vu7q6Ij09HUePHsVrr72GpUuXYs+ePdd8zBdeeAHV1dWWr4KCgi48EyLbY57/c3NLLwH1Px5OKowN9QTA1WBEPyfsI5+Pjw8UCkWb3p6SkpI2vTxmAQEB7bZXKpXw9r7afS+XyzFo0CAAwOjRo5GZmYnExETMmDGj3cdVq9VQq9XdeDZEtqe2qRlHcysAcP+f/u7mKD8cy6vE7rOleHBCmOhyiKyCsB4glUqF2NhYJCcnt7qenJyMSZMmtfszEydObNN++/btiIuLg4PDtXevlSQJWi13QiX6uYPny6A3Soj0cUa4j7PocqgXmecBHTxfBq3eILgaIusgdAhs6dKl+PDDD7FhwwZkZmZiyZIlyM/PR0JCAgDT0NSiRYss7RMSEpCXl4elS5ciMzMTGzZswPr167Fs2TJLm8TERCQnJ+PixYs4e/Ys3nnnHXz66ad46KGH+vz5EVmzq8vf2fvT3w0LdIOfqxqNzQYcyakQXQ6RVRA663HhwoUoLy/HqlWrUFRUhJiYGCQlJSEszNRFW1RU1GpPoIiICCQlJWHJkiV4//33ERQUhNWrV2PBggWWNvX19Xj66adRWFgIR0dHREVFYePGjVi4cGGfPz8iayVJ0tXjL6I4/6e/k8lkuHmoHzYdK8Dus6WYOph/50RC9wGyVtwHiPq705eqccffD8BJpUDan2+DWqkQXRL1sh9PFyNhYyoifJyxe9kM0eUQ9Qqb2AeIiMQx7wkzeZAPw4+dmDzIGw4KGXLK6pFTVi+6HCLhGICI7JD5dHCu/rIfrhoHjAv3AsBNEYkABiAiu1NZr0NavulcqBnc/8eumAOvOQAT2TMGICI7sy+7FEYJiApwRZCHo+hyqA+ZJ7z/dLEcDTq94GqIxGIAIrIzu1uOv+Dyd/sz0NcFIV6O0OmNOHS+XHQ5REIxABHZEYNRwt5zpuGPW6IYgOyNeTk8wGMxiBiAiOzIicIqVDY0w1WjxNhQD9HlkADmALQnqxTcBYXsGQMQkR0xD39NG+ILpYL//e3RTZHeUCvluFTViOySOtHlEAnDV0AiO2LZ/Znzf+yWo0qBiQNNh0ebAzGRPWIAIrITJTVNOH2pBgCXv9s7cwDexQBEdowBiMhO7GmZ/DxqgDt8XNSCqyGRzAHoWF4lapqaBVdDJAYDEJGd4PJ3Mgv1dsJAX2cYjBIOZJeJLodICAYgIjvQbDBif8sb3c1c/k742a7QHAYjO8UARGQHjuVWok6rh7ezCiOD3UWXQ1bAHIR3Z5XCaORyeLI/DEBEdsB8+OX0ob6Qy2WCqyFrMC7cC84qBcrqtDhTVCO6HKI+xwBEZAe4/J1+SaWUY/IgHwAcBiP7xABE1M8VVjbg3JU6yGXAtMFc/k5XXR0GYwAi+8MARNTP7c4yLX+PDfOEu5OD4GrImpj3g0orqEJFvU5wNUR9iwGIqJ/b0zK8wdVf9EuB7o6ICnCFJAH7s0tFl0PUpxiAiPqxpmYDDl5oWf7O+T/UDsswGOcBkZ1hACLqxw7nVKCp2YgANw2iAlxFl0NWyByM954rhYHL4cmOMAAR9WO7LcNfvpDJuPyd2hob6gFXjRKVDc04UVgluhyiPsMARNRPSZJkOeySx1/QtSgVckwbYpoMvYfDYGRHGICI+qmcsnrkVzTAQSHDlJb9XojaYzkWI4sTocl+MAAR9VPmN7MJEd5wVisFV0PWbHpLD9CpS9UoqW0SXA1R32AAIuqnrp7+zs0P6fp8XdUYOcB0Rtxe9gKRnWAAIuqH6rV6HM4pB8D9f6hjzPPE9jAAkZ1gACLqhw6eL0OzQUKYtxMifZxFl0M24OaWnsJ92aVoNhgFV0PU+xiAiPoh8/yfm4f6cfk7dcjIAR7wclahtkmP43mVossh6nVdCkBZWVl4+eWXMXPmTAwcOBCBgYEYOXIkHnnkEXz22WfQarU9XScRdZAkSdiTxfk/1DkKucwyGZqrwcgedCoApaWl4bbbbsOoUaOwb98+jBs3DosXL8arr76Khx56CJIkYeXKlQgKCsKbb77JIEQkwNniWhRVN0HjIMdNkd6iyyEbYg7Me3g6PNmBTq2NnT9/Pv74xz9i06ZN8PLyuma7lJQU/PWvf8X//d//YcWKFd0ukog6bnfLm9ekgT7QOCgEV0O2ZNpgX8hkphB9uaoRQR6Ooksi6jWdCkDZ2dlQqVQ3bDdx4kRMnDgROp2uy4URUdfsOdsy/4erv6iTPJ1VGBPigeP5VdiTVYoHJoSKLomo13RqCOzn4efTTz9td4hLp9Ph008/bdOeiHpfdUMzUvNNE1hnDOH8H+q8q7tCcxiM+rcurwJ77LHHUF1d3eZ6bW0tHnvssW4VRURdsy/bdKL3YD8XhHg5iS6HbJC55/Dg+TJo9QbB1RD1ni4HIEmS2l1eW1hYCHd3924VRURdY/7UzuEv6qphgW7wdVWjQWfA0Rwuh6f+q9MHBI0ZMwYymQwymQwzZ86EUnn1IQwGA3JycnD77bf3aJFEdGNGo2Q5xuBmnv5OXSSXyzBjiC++TC3E7qwSTBnMg3Spf+p0AJo/fz4AID09HbNmzYKLi4vleyqVCuHh4ViwYEGPFUhEHXPqUjXK63VwUSsRF+4puhyyYTdH+VkC0J/uGCa6HKJe0ekA9NJLLwEAwsPDsXDhQmg0mh4viog6zzz8NXWwDxwU3OSdum7KYB8o5DJcLK1HfnkDQr05n4z6ny6/Sj7yyCM3DD+SJHX14Ymok8ynv3P4i7rLTeOAuDBTL+Kec1wNRv1TpwJQdHQ0Pvvssxvu75OdnY3f/va3ePPNN7tVHBF1TGmtFicKTasyefwF9QTzRHpzsCbqbzo1BPb+++9j+fLl+N3vfof4+HjExcUhKCgIGo0GlZWVOHPmDA4cOIAzZ87gmWeewdNPP91bdRPRz+w7Z5r8HBPsBj83DktT99081A9v/HAWhy6Uo6nZwF3Fqd/pVAC65ZZbcPToURw6dAibNm3CZ599htzcXDQ2NsLHxwdjxozBokWL8NBDD8HDw6OXSiaiX7Isf+fwF/WQIf4uCHLX4HJ1E1IulvPfFvU7nZ4EDQCTJk3CpEmTeroWIuoCvcFo6QGawTcp6iEymQwzovzw2eF87DlbwgBE/U6nJ0HPmTOn1Q7Qr732Gqqqqix/Li8vx7BhXDZJ1FeO51ehpkkPTycHjA7xEF0O9SNXj8Uo5aIW6nc6HYC2bdvW6gywN998ExUVFZY/6/V6ZGVl9Ux1RHRD5uGv6UN8oZC33Z2dqKsmDfSGSiFHfkUDLpbViy6HqEd1OgD98lMAPxUQiWVepcPhL+ppzmolJkR6AeBqMOp/uFsakQ27XNWIs8W1kMtMPUBEPc3872pPyzErRP1FpwOQ+RywX14jor5nHv4aE+oJT2eV4GqoPzLvB3Q4pxz1Wr3gaoh6TqdXgUmShEcffRRqtRoA0NTUhISEBDg7OwNAq/lBRNS7zMMSt/D0d+olkT7OCPVyQn5FAw6eL0P88ADRJRH1iE73AD3yyCPw8/ODu7s73N3d8dBDDyEoKMjyZz8/PyxatKjDj7dmzRpERERAo9EgNjYW+/fvv277vXv3IjY2FhqNBpGRkVi3bl2r73/wwQeYOnUqPD094enpiVtvvRVHjhzp7NMksnpNzQYcPF8OgPv/UO+RyWS4uWV38d0cBqN+pNM9QB999FGP/fJNmzZh8eLFWLNmDSZPnox//OMfmD17Ns6cOYPQ0NA27XNycjBnzhw89dRT2LhxIw4ePIinn34avr6+lhPo9+zZg/vvvx+TJk2CRqPBW2+9hfj4eGRkZCA4OLjHaicS7aeL5WhsNiDATYPoQFfR5VA/NiPKD5+k5GFPVgkkSeK0B+oXZJLAZVwTJkzA2LFjsXbtWsu16OhozJ8/H4mJiW3aL1++HFu3bkVmZqblWkJCAk6cOIGUlJR2f4fBYICnpyfee++9DvdM1dTUwN3dHdXV1XBzc+vksyLqGy9vzcDHh3Jx//gQJN49UnQ51I81NRsw6pXt0OqN+HHxVEQF8HWRrFNn3r+FrQLT6XRITU1FfHx8q+vx8fE4dOhQuz+TkpLSpv2sWbNw7NgxNDc3t/szDQ0NaG5uhpeXV88UTmQFJEnCLp7+Tn1E46DApIHeAIDdZzkMRv2DsABUVlYGg8EAf3//Vtf9/f1RXFzc7s8UFxe3216v16OsrKzdn3n++ecRHByMW2+99Zq1aLVa1NTUtPoismYXSuuRX9EAlUKOyYN8RJdDdsByOnwW9wOi/kH4PkC/HEu+0fhye+3buw4Ab731Fj7//HNs2bIFGs21T8hOTEy0TOJ2d3dHSEhIZ54CUZ8zr/6aEOkFZ3WXjvQj6pQZQ0wBKDWvEtWN7fe4E9kSYQHIx8cHCoWiTW9PSUlJm14es4CAgHbbK5VKeHt7t7r+9ttv4/XXX8f27dsxcuT150e88MILqK6utnwVFBR04RkR9R2e/k59LdTbCQN9nWEwSjiQ3X6PO5EtERaAVCoVYmNjkZyc3Op6cnLyNU+anzhxYpv227dvR1xcHBwcHCzX/vKXv+DVV1/Fjz/+iLi4uBvWolar4ebm1uqLyFrVNjXjSI7p/D3u/0N96erhqBwGI9sndAhs6dKl+PDDD7FhwwZkZmZiyZIlyM/PR0JCAgBTz8zPV24lJCQgLy8PS5cuRWZmJjZs2ID169dj2bJlljZvvfUWXnzxRWzYsAHh4eEoLi5GcXEx6urq+vz5EfWGA9ll0BslRPo4I9zHWXQ5ZEfM84D2ZJXCaOQ5kGTbhE4eWLhwIcrLy7Fq1SoUFRUhJiYGSUlJCAsLAwAUFRUhPz/f0j4iIgJJSUlYsmQJ3n//fQQFBWH16tWWPYAA08aKOp0O99xzT6vf9dJLL+Hll1/uk+dF1Jssq7/Y+0N9LC7cE84qBcrqtMi4XIMRA9xFl0TUZUL3AbJW3AeIrJXRKGH86ztRVqfFxicmYMpgrgCjvvXrT49h+5krWHrbEDw7c7DocohasYl9gIio8zIu16CsTgtnlQLjI7i3FfW9q8NgnAdEto0BiMiGmIe/pgz2gUrJ/77U92a0nAuWVlCFinqd4GqIuo6voEQ2ZFcWT38nsQLdHREV4ApJAvZnc1dosl0MQEQ2oqxOi5OFVQCAGdz/hwQy//szb8hJZIsYgIhsxN6sUkgSMDzIDf5u197ZnKi33dwyDLb3XCkMXA5PNooBiMhGcPiLrMXYME+4apSobGjGiZZeSSJbwwBEZAOaDUbsO2eab8H9f0g0B4Uc0wabeoH2cBiMbBQDEJENSM2rRG2THl7OKowa4CG6HCLLarDdWZwITbaJAYjIBpgnm04f4guFXCa4GiJgeksAOnWpGiU1TYKrIeo8BiAiG7CTx1+QlfFz1WBUy1EYOzkMRjaIAYjIyuWW1eN8SR2Ucpll2IHIGtwa7Q8A2Jl5RXAlRJ3HAERk5Xa0vLlMiPSCm8ZBcDVEV81sCUD7s8vQqDMIroaocxiAiKycOQDNjPIXXAlRa9GBrgj2cIRWb8TB82WiyyHqFAYgIitW3dCMo7mVAK4ONxBZC5lMhpnRpnlpOzgMRjaGAYjIiu05VwKDUcIQfxeEejuJLoeoDcs8oLMlMHJXaLIhDEBEVmxHpml1zUz2/pCVmhDpBWeVAqW1Wpy8VC26HKIOYwAislLNBiP2tBx/weEvslZqpcKyJxBXg5EtYQAislJHcypQ26SHt7MKo0M8RJdDdE3mCfrmHksiW8AARGSlzG8mt0T5cfdnsmo3R/lBLgMyi2pQWNkguhyiDmEAIrJCkiRdXf7O4S+ycl7OKsSGeQIAdnFXaLIRDEBEVuh8SR3yKxqgUsoxdbCP6HKIbsgc1JPPcB4Q2QYGICIrZB7+mjTQG85qpeBqiG7MPFH/p4vlqG1qFlwN0Y0xABFZIQ5/ka0Z6OuMcG8nNBsk7M/mrtBk/RiAiKxMeZ0Wx/PNuz/z9HeyDTKZzNILxF2hyRYwABFZmV1nSyBJwPAgNwS6O4ouh6jDzD2Wu8+adjAnsmYMQERWZmcmNz8k2xQX7gl3RwdUNjRbejGJrBUDEJEV0eoN2JddCoABiGyPg0KOGS27Qu/gajCycgxARFbk0PlyNOgM8HdTIybYTXQ5RJ1mDu7JnAdEVo4BiMiKbMsoBgDEDwuATMbdn8n2TB/qCweFDBdL63G+pE50OUTXxABEZCUMRsmyidys4QGCqyHqGjeNAyYONG3eaQ70RNaIAYjISqTmVaK8Xgd3RwdMiPQSXQ5Rl80abhoG284ARFaMAYjISpg/Lc+M8oODgv81yXbdNswfMhlworAaRdWNosshahdfZYmsgCRJV+f/cPiLbJyfqwZjQ02Ho27P4GRosk4MQERW4ExRDQorG6FxkGP6EF/R5RB12+0tQZ7zgMhaMQARWQHzp+Spg33hqFIIroao+8wT+Q/nVKCyXie4GqK2GICIrID5UzJXf1F/EerthKgAVxiMEs8GI6vEAEQkWH55A84W10Ihl/HwU+pXZlmGwRiAyPowABEJZu79mRDhBQ8nleBqiHqOOQDtzy5Fg04vuBqi1hiAiATj8Bf1V9GBrgj1coJWb8TerFLR5RC1wgBEJFBprRapLadmxw/n4afUv8hkMsumiFwNRtaGAYhIoOQzVyBJwKgB7gh0dxRdDlGPM/ds7jxbAp3eKLgaoqsYgIgE4uaH1N+NDfWEj4satU16/HSxXHQ5RBYMQESCVDc249CFMgBXz04i6m/kchluG2b69/0jh8HIijAAEQmy48wVNBskDPF3wSA/V9HlEPWa22NMPZzbM67AYJQEV0NkwgBEJEjSqSIAwJwRgYIrIepdEyO94aZRoqxOi2O5FaLLIQLAAEQkRE1TM/Znm4a/GICov1Mp5ZZ5bubgTyQaAxCRADvOXIHOYMQgPxcM8efwF/V/c1uCftLpYg6DkVVgACISgMNfZG8mD/KBm0aJ0loOg5F1YAAi6mM1Tc3Yd840/DWXAYjsxM+Hwf7LYTCyAgxARH1sZ6Zp+GugrzOG+LuILoeoz8wdaQr8P3AYjKwAAxBRH/vvSdNeKHNHBEImkwmuhqjvTB7IYTCyHsID0Jo1axAREQGNRoPY2Fjs37//uu337t2L2NhYaDQaREZGYt26da2+n5GRgQULFiA8PBwymQzvvvtuL1ZP1Dm1Tc3Yl206FHLOSA5/kX3hMBhZE6EBaNOmTVi8eDFWrlyJtLQ0TJ06FbNnz0Z+fn677XNycjBnzhxMnToVaWlpWLFiBZ599lls3rzZ0qahoQGRkZF44403EBDA4wXIuuxqOQ8p0tcZQ7n6i+wQh8HIWggNQO+88w6eeOIJPPnkk4iOjsa7776LkJAQrF27tt3269atQ2hoKN59911ER0fjySefxOOPP463337b0mbcuHH4y1/+gvvuuw9qtbqvngpRh/z3ZMvqrxgOf5F9+vkw2FEOg5FAwgKQTqdDamoq4uPjW12Pj4/HoUOH2v2ZlJSUNu1nzZqFY8eOobm5uddqJeoJdVo99pxrGf7i6i+yUyql3HJCPDdFJJGEBaCysjIYDAb4+7c+BNLf3x/Fxe0fmFdcXNxue71ej7Kysi7XotVqUVNT0+qLqKftOHMFOr0RET7OiA7k8BfZL/P8t6RTHAYjcYRPgv7lMIAkSdcdGmivfXvXOyMxMRHu7u6Wr5CQkC4/FtG1fJt+CQBw56ggDn+RXZs80Afujg4oq9PiSA6HwUgMYQHIx8cHCoWiTW9PSUlJm14es4CAgHbbK5VKeHt7d7mWF154AdXV1ZavgoKCLj8WUXsq6nWWs7/mjQoSXA2RWKZhMNPr/HcnLwuuhuyVsACkUqkQGxuL5OTkVteTk5MxadKkdn9m4sSJbdpv374dcXFxcHBw6HItarUabm5urb6IelLSqSLojRKGB7lhkB83PySaNyoYgOn/hk5vFFwN2SOhQ2BLly7Fhx9+iA0bNiAzMxNLlixBfn4+EhISAJh6ZhYtWmRpn5CQgLy8PCxduhSZmZnYsGED1q9fj2XLllna6HQ6pKenIz09HTqdDpcuXUJ6ejrOnz/f58+PyGxruulT7l2j2ftDBAATB3rD11WNqoZm7GtZHEDUl4QGoIULF+Ldd9/FqlWrMHr0aOzbtw9JSUkICwsDABQVFbXaEygiIgJJSUnYs2cPRo8ejVdffRWrV6/GggULLG0uX76MMWPGYMyYMSgqKsLbb7+NMWPG4Mknn+zz50cEAJerGnEktwIymWn+DxEBCrkMd440/X/49gSHwajvySTzLGKyqKmpgbu7O6qrqzkcRt32j70XkPjDWYyP8MJ/fjNRdDlEVuNkYRXmvXcQGgc5jr14G1zUStElkY3rzPu38FVgRP3d1hMc/iJqz4hgd0T4OKOp2YjkM+1vf0LUWxiAiHrR+ZJaZFyugVIuw5wYbn5I9HMymczyweCbNA6DUd9iACLqRebJz9OG+MLTWSW4GiLrc9do02qwA+fLUFanFVwN2RMGIKJeIkkSh7+IbiDCxxmjBrjDYJQsZ+UR9QUGIKJecqKwGrnlDdA4yHFrdPubexIRMK+lF8i8WzpRX2AAIuolW44XAgDihwXAmatbiK7pzpGBkMuA4/lVyC9vEF0O2QkGIKJeoNUbLMNfC2IHCK6GyLr5uWkwaaAPAPYCUd9hACLqBbvPlqCqoRn+bmpMGeQjuhwiqzd/jGkYbPPxQnB7OuoLDEBEveCrVNOn2PljgqGQ8+R3ohuZMyIAzioFcssbcCyvUnQ5ZAcYgIh6WHmdFnuySgAA94zl8BdRRziplJgzwrRX1lfHCgVXQ/aAAYioh32bfhl6o4SRA9wx2N9VdDlENuOelvly/z1VhAadXnA11N8xABH1sM0tq78WsPeHqFPGhXsh1MsJdVo9tmXwaAzqXQxARD3obHENMi7XwEEhwzye/E7UKXK5zPLB4UsOg1EvYwAi6kGbU00v2jOj/Hn0BVEXLIg1rQY7dKEchZXcE4h6DwMQUQ9pNhjxdRr3/iHqjgGeTpg00BsAsOU49wSi3sMARNRDdp0tQVmdFj4uKswY6iu6HCKbZZ4M/VUq9wSi3sMARNRDvjiSD8DU++Og4H8toq66PSYALmol8isacDinQnQ51E/xVZqoB1yqasSec6UAgPvGhQquhsi2OamUuGOkaU8g8wcLop7GAETUA/5ztACSBNwU6YUIH2fR5RDZvPvHmz5IJJ0uRmW9TnA11B8xABF1k8Eo4ctjBQCuvmgTUfeMHOCO4UFu0OmNlr21iHoSAxBRN+07V4rL1U3wcHLArOEBossh6hdkMhkemGD6QPHZkXxOhqYexwBE1E2ftcxRuHvMAGgcFIKrIeo/7hodDGeVAhdL6/HTRU6Gpp7FAETUDSU1Tdh11nTw6f3jQwRXQ9S/uKiVmDfatDHiZ5wMTT2MAYioG/5zrAAGo4TYME8efErUCx5omVf34+kilNdpBVdD/QkDEFEX6Q1G/Puw6VPpA5z8TNQrRgxwx4hgdzQbJHyVysnQ1HMYgIi6KPnMFRRVN8HbWYW5LXuWEFHPM0+G/vxIPoxGToamnsEARNRFn6TkAgDuGx/Cyc9EvWjeqCC4qpXILW/A3uxS0eVQP8EARNQFWcW1+OliBRRyGR6cECa6HKJ+zVmtxL3jTIsMPj6YK7YY6jcYgIi6wNz7Ez/MH0EejmKLIbIDj0wMh0wG7D1XivMldaLLoX6AAYiok6obm/H18UsAgEUTw8UWQ2QnQr2dMDPKHwDwacsHEKLuYAAi6qSvUgvR2GzAUH9X3BTpJbocIrvx2ORwAKb/g9WNzWKLIZvHAETUCQajZPn0uWhSGGQymdiCiOzIpIHeGOLvggadwXL+HlFXMQARdULymWLklTfA3dEBvxoTLLocIrsik8nw6KQIAMCnKXkwcEk8dQMDEFEHSZKEf+y7CAB4+KYwOKmUgisisj+/GhMMd0cH5Fc0YGfmFdHlkA1jACLqoNS8SqTlV0GlkOORSeGiyyGyS44qhWVjRPMHEqKuYAAi6qB/trzY3j02GL6uasHVENmvxyaFQ6WQIzWvEkdzeUo8dQ0DEFEHXCytQ3JLd/uTUyMEV0Nk3/zcNFgQa5qDt27PBcHVkK1iACLqgA8P5ECSgFuj/TDIj6e+E4n21NRIyGTAzrMlyCquFV0O2SAGIKIbKKltwuaWU6ifmhopuBoiAoBIXxfcPjwAAPCPfewFos5jACK6gQ/350CrN2JMqAfGR3DjQyJrkTB9IABga/plXKpqFFwN2RoGIKLrKK/T4l8peQCAZ28ZzI0PiazIqBAPTIz0ht4o4QOuCKNOYgAiuo71B3LQ2GzAiGB3zBjqK7ocIvqF3908CADw2ZF8XKlpElwN2RIGIKJrqGrQ4ZNDuQCA398yiL0/RFZo8iBvxIV5Qqc3Yi1XhFEnMAARXcOGg7mo1xkQHeiG24b5iy6HiNohk8mw5LYhAEy9QMXV7AWijmEAImpHdWMzPjqYA4C9P0TWbtJAb4wP94JOb8SaPedFl0M2ggGIqB3/2HsBtU16DPG/utSWiKyTTCbD4lsHAwC+OFKAomquCKMbYwAi+oWSmiZsaOn9+eOsKMjl7P0hsnYTB3pjfIQXdAYj3t/NXiC6MQYgol/4285sNDUbMTbUA7dG+4kuh4g6QCaTYWnLXKAvjhTgYmmd4IrI2jEAEf1MTlk9vjhaAABYfnsU5/4Q2ZCbIr1xS5Qf9EYJb/2YJbocsnIMQEQ/807yORiMEmYM9cWESG/R5RBRJz0/OwpyGfBjRjFS83hSPF0bAxBRi7T8Snx34jIA4H9nRQmuhoi6Yoi/K+6NCwEAvJ50FpIkCa6IrJXwALRmzRpERERAo9EgNjYW+/fvv277vXv3IjY2FhqNBpGRkVi3bl2bNps3b8awYcOgVqsxbNgwfP31171VPvUTRqOEV747AwC4e2wwhgW5Ca6IiLpqyW1D4OigQGpeJbZlFIsuh6yU0AC0adMmLF68GCtXrkRaWhqmTp2K2bNnIz8/v932OTk5mDNnDqZOnYq0tDSsWLECzz77LDZv3mxpk5KSgoULF+Lhhx/GiRMn8PDDD+Pee+/F4cOH++ppkQ36Jv0S0guq4KxS4Pnb2ftDZMv83TR4amoEAOC1pEw0NRsEV0TWSCYJ7B+cMGECxo4di7Vr11quRUdHY/78+UhMTGzTfvny5di6dSsyMzMt1xISEnDixAmkpKQAABYuXIiamhr88MMPlja33347PD098fnnn3eorpqaGri7u6O6uhpubuwJ6O/qtXrc/PYelNRqsfz2KPx2xkDRJRFRN9Vr9bj1nb0oqm7C4lsHY/GtQ0SXRH2gM+/fwnqAdDodUlNTER8f3+p6fHw8Dh061O7PpKSktGk/a9YsHDt2DM3Nzddtc63H7Eumk8Vz8eWxAnx/8jJ2nb2C05eqUVmv4zi1QO/vPo+SWi3CvJ3w+JRw0eUQUQ9wVivx4txhAIC1ey4gv7xBcEX2S6s3ILesHsdyK/Dj6SJ8lVqIzw7n48fTYocnlaJ+cVlZGQwGA/z9W5+x5O/vj+Li9m9KcXFxu+31ej3KysoQGBh4zTbXekwA0Gq10Gq1lj/X1NR09ul0SH5FA/70bUa733NSKRAV4IqYYHfEBLtjYqQ3QryceqUOuup8SS0+2H8RALBiTjTUSoXgioiop8wZEYDJg7xx8Hw5Vn2fgQ8fGSe6pH5PpzfiZGEVDudU4FRhNc6V1CKvvAEGY9sP+bFhnrg9RtxO+8ICkNkv91mRJOm6e6+01/6X1zv7mImJiXjllVc6XHNXuWqUmB0TgMZmAxp1BtTr9Ciu1qKsTosGnQHH86twPL/K0j7c2wlTBvvg9uGBuCnSC0qF8Dnr/YrRKOH5zafQbJBwS5Qf4nngKVG/IpPJ8Mq8GMz+2z7syCzBjjNXcCv/n/e40lotks9cwbaMYhzOKUdTs7FNG0cHBXxd1fBxUcFV4wAHhRyD/FwEVHuVsADk4+MDhULRpmempKSkTQ+OWUBAQLvtlUolvL29r9vmWo8JAC+88AKWLl1q+XNNTQ1CQkI69Xw6YpCfK9Y+FNvmelOzAQUVDThTVIPTl6pxPL8K6QVVyC1vQG55Pjb+lA9vZxVmjwjAgrEDMDrEgxv09YDPj+bjWF4lnFQKvDo/hveUqB8a5OeCx6dE4B97L+JP357G+EgvuGkcRJdl8xp0eiSdKsaXxwpwJLcCP5/F4eWswvhwL8SGeWJogCuGBrjCz1Vtda+xwgKQSqVCbGwskpOT8atf/cpyPTk5GXfddVe7PzNx4kR89913ra5t374dcXFxcHBwsLRJTk7GkiVLWrWZNGnSNWtRq9VQq9XdeTrdonFQYLC/Kwb7u+Ku0cEAgNqmZvx0sQK7zl7Bj6eLUV6vw8afTGFoeJAbHpwQhrtGB8FZLbwTzyZdqWnCG0lnAQDL4oci2MNRcEVE1FsWzxyCH08XI6+8Aa//NxNvLBgpuiSblVlUg40/5WFr+mXUavWW66MGuGNWTABmRvljiL+L1YWd9ghdBbZp0yY8/PDDWLduHSZOnIh//vOf+OCDD5CRkYGwsDC88MILuHTpEj799FMApmXwMTEx+M1vfoOnnnoKKSkpSEhIwOeff44FCxYAAA4dOoRp06bhtddew1133YVvv/0WL774Ig4cOIAJEyZ0qC5rWwXWbDDi0IVyfJt2Cd+fKoJOb+pedNUo8fBNYXh8SgR8XMQFOFsjSRKe+jQVOzKvYNQAd2x5ejIUPPCUqF87fLEcC//5EwDgX0+Mx9TBvoIrsh2SJCHlQjnW7buIfedKLddDvZywcFwI5o8JtpoPkZ15/xYagADTRohvvfUWioqKEBMTg7/+9a+YNm0aAODRRx9Fbm4u9uzZY2m/d+9eLFmyBBkZGQgKCsLy5cuRkJDQ6jG/+uorvPjii7h48SIGDhyI1157DXfffXeHa7K2APRzVQ06fJVaiH8fzkdOWT0AQK2UY+G4EDw1NZITpzvgiyP5eH7LKagUcnz7zGREB1rX3zER9Y6Xvj2NT1LyEOzhiG1LpsGFPejXJUkStmVcwfu7z+PUpWoAgFwGzI4JxIMTQnFTpDfkVvbh0aYCkDWy5gBkZjRK2JF5BWv2XEB6QRUAwEEhwwPjQ/HMLYPh68oeofbkltVjzur9aNAZ8MLsKPxmOvf8IbIX9Vo9bv/bPhRUNOKe2AF4+39GiS7JKkmShP3ZZXh7exZOFpqCj8ZBjnvjQvDklEiEelvvB20GoG6yhQBkJkkSUi6W4/3d53HwfDkA05L6xydH4NfTIznZ72f0BiP+5x8pSMuvwoQIL3z21E0c+iKyM0dyKnDfP1NglIC/LhyFX40ZILokq3IstwJ/2ZaFwzmmg2TN7yePT4mAl7NKcHU3xgDUTbYUgH7u0PkyvLktCydaeoQ8nByweOZgPHhTGBy4hB5v/ngWa/dcgKtaiR8WT8UAT+v9FENEvefdHefw7o5sOKsU+P7ZqYjwcRZdknCFlQ14PSkTSadMq6hVSjkevikMv50x0KbmmDIAdZOtBiDg6pjt29uzcL6kDgAw2M8Ff7pjGKYNsd9Jf9szivHrf6UCAP5+/xjcOSpIcEVEJIrBKOGBD37C4ZwKxAS74auESdA42OcmqI06A9btvYB1ey9AqzdCLgPujQvBszMHI8hKJjZ3BgNQN9lyADLTG4zYdKwAb2/LQmWD6ZiQW6P9sHLuMLv7tJNbVo873zuA2iY9HpscjpfuHC66JCISrKi6EXP+th+VDc24e0ww/u/eUTaxdLunSJKEpFPFeD0pE5eqGgEAN0V64aU7h9v0whAGoG7qDwHIrLqhGX/bmY1PU3KhN0pwUMjw2OQI/P6WQXC1g/lBtU3NuGdtCrKu1CIuzBOf//omDgcSEQDg4PkyLNpwBAajhJVzovHUtEjRJfWJs8U1eHlrBn66aJrnE+zhiJVzozE7JsDmQyADUDf1pwBkdr6kDq/99wx2Z5n2cPBxUeF/Z0XhntgBVreMsac0G4x4/OOj2J9dBl9XNb57ZgoC3DWiyyIiK/LxwRy8/N0ZyGXA+kfH4eahfqJL6jVVDTq8k3wOG3/Kg1EybaGSMH0gEqYPhKOqfwwBMgB1U38MQGa7z5bg1e/P4GLLHkKjBrjjpXnDMTbUU3BlPUuSJLyw5RS+OFoARwcF/vObiRgxwF10WURkZX7+WuGkUuCzp27C6BAP0WX1KINRwmdH8vF/27NQ1TIlYs6IAKyYE93vFoMwAHVTfw5AgOm03k8O5eJvO7NR17KV+d1jgrF8dhT83fpHD8nb27Lw3u7zkMuAfz4cxwMQieiadHojnvjE1Fvs4eSArxImYpCfq+iyesThi+V4+bszyCyqAQAM9XfFS/OGYdJAH8GV9Q4GoG7q7wHIrKS2CX/5MQtfphYCAJxVCvzulkF4YkoE1Erb7Q5dvTMb7ySfAwC8etdwPDwxXGxBRGT16rV6PPDhYZwoqEKAmwabfnMTwrxtd8HIpapGvJ6Uif+eLAIAuDs64A/xQ/DA+FAo+/E8SAagbrKXAGR2oqAKL3+XgbT8KgBAmLcT/jR3GGZG+9nUhDhJkvD+7vN4e7sp/Lw4NxpPTrWPSY1E1H2V9Trc+48UZJfUwc9VjX8/OQGD/W2rJ6hRZ8A/9pmWtTc1m5a13z8+FH+IH2oTGxl2FwNQN9lbAAJMR2t8nXYJb/x4FqW1WgDAtCG++PMdwzDIz0VwdTdmNEpY9f0ZfHwoFwDwv7cPxdMzBoktiohsTkltEx7+8AiyrtTCy1mFTx8fj5hg658/2N6y9gkRpmXtw4Ls430MYADqNnsMQGZ1Wj3e23Ue6w9cRLNBglIuwyOTwvHcrYOt9liNRp0By746Yenq/dMdw/DElAjBVRGRraqs12HRhiM4dakaTioF/rpwNGYNDxBd1jUdza3Amz+cxbG8SgCmZe0r5kRjzgjbX9beWQxA3WTPAcgsp6we/+/7M9h5tgQA4OWsQsL0SDx8U7hVLZfML2/AbzamIrOoBg4KGd7+n1G4a3Sw6LKIyMbVNDXj6Y3HceB8GQBgWfwQPD1jkFVtG5JVXIu3fjxreZ3WOJiWtf9mWv9Z1t5ZDEDdxAB01Z6sEqz6/gwulpqWzfu6qvH0jIG4f3yo8K3jk04V4fnNJ1HTpIe3swrvPTAWEwd6C62JiPoPvcGI//ffTMvQ+pRBPvi/e0cJXy17vqQWa3ZfwNfplyBJgEIuw71xIVh862DhtYnGANRNDECt6Q1GbDl+Cat3ZaOw0jS2HOCmwRNTIrBwfEifD42V12nx560ZliGvMaEeWPPgWAS62965NURk/TYdzcdLWzPQ1GyEh5MD/jR3GH41JrjPe4NOFVZjzZ7z+DGjGOZ37jkjAvCH+KEY6Gv9czX7AgNQNzEAtU+nN+LL1AK8t+s8iqqbAJiWzt87LgSPTYpAqHfvbqjV1GzAJ4dy8d7u86ht0kMhl+HpGQPx+1sGQ6Xsv8s6iUi88yV1WLwpDacvmfbTGRvqgZfuHI5RvbxpYrPBiO0ZV7DxpzykXCy3XI8f5o9nbhmEkQN69/fbGgagbmIAuj6t3oBv0i5h/YEcnLtSZ7k+IcILC2IHYM6IQLiolT32+2qamrHpSAE2HMyxBK/hQW54c8FIm1idQUT9g05vxPoDOfj7rmw06AwAgBlDffH0jEEYF+7ZYxOOJUnC2eJafHfiMr5MLbSszFXIZZg3Kgi/nTEQQ2xseX5fYQDqJgagjpEkCfuzy/DhgRzszy61dMlqHOSYNNAHN0f5YcYQXwzwdOz0C0NTswEpF8rx3YnL2JZRjPqWF5sANw3+OGuokO5nIiLAdJL8X37Mwjfpl2Bsed0b6OuMe2JDcNswfwz0de7Sa97xvEocvFCGH08X40LLvEsA8HFR475xIbh/QiiCPTjUfz0MQN3EANR5l6sa8XXaJWw+XmiZMG3m66rGqAEeGBbkhhBPRwR7OsLd0QEaBwUUMhnqtHrUNDYjv6IBOWX1SCuoQnpBFXR6o+Uxhvi74MkpkZg3Okj45GsiIgDIK6/Hur0X8E3aZTQ2GyzXgz0cERvmiaEBrhjo6wwfFzU8nFRQyGUwGCXUafUordWiqLoRZ4trcbaoBqcv17R6zVMp5ZgxxBd3jQ7GbcP8OczfQQxA3cQA1HWSJCHrSi12nS3B7rMlOJ5fBYOxa//E/FzVmDMiEHeMDERsWM91LxMR9aTapmYknSrC9yeLcDinolWQ6Qw/VzUmD/LB9CG+mBntB1cr3XvNmjEAdRMDUM9p1BlwpqgaaflVuFBaj8LKBlyqakRdkx5NzQYYjBJcNEq4qJUI9nRCpI8zhga4YkKEFyJ8Ot+NTEQkUqPOgKO5FTh9uRpZxbXIr2hARb0OlfU6SADkMhmcVQr4umng56rGYD8XRAW6YXiQGyL5mtdtDEDdxABERERkezrz/s1BRSIiIrI7DEBERERkdxiAiIiIyO4wABEREZHdYQAiIiIiu8MARERERHaHAYiIiIjsDgMQERER2R0GICIiIrI7DEBERERkdxiAiIiIyO4wABEREZHdYQAiIiIiu8MARERERHZHKboAayRJEgCgpqZGcCVERETUUeb3bfP7+PUwALWjtrYWABASEiK4EiIiIuqs2tpauLu7X7eNTOpITLIzRqMRly9fhqurK2QyWY8+dk1NDUJCQlBQUAA3N7cefWy6ive5b/A+9w3e577De903eus+S5KE2tpaBAUFQS6//iwf9gC1Qy6XY8CAAb36O9zc3Pifqw/wPvcN3ue+wfvcd3iv+0Zv3Ocb9fyYcRI0ERER2R0GICIiIrI7DEB9TK1W46WXXoJarRZdSr/G+9w3eJ/7Bu9z3+G97hvWcJ85CZqIiIjsDnuAiIiIyO4wABEREZHdYQAiIiIiu8MARERERHaHAagPrVmzBhEREdBoNIiNjcX+/ftFl2TTEhMTMW7cOLi6usLPzw/z589HVlZWqzaSJOHll19GUFAQHB0dMWPGDGRkZAiquH9ITEyETCbD4sWLLdd4n3vGpUuX8NBDD8Hb2xtOTk4YPXo0UlNTLd/nfe4Zer0eL774IiIiIuDo6IjIyEisWrUKRqPR0ob3uvP27duHO++8E0FBQZDJZPjmm29afb8j91Sr1eL3v/89fHx84OzsjHnz5qGwsLB3CpaoT3zxxReSg4OD9MEHH0hnzpyRnnvuOcnZ2VnKy8sTXZrNmjVrlvTRRx9Jp0+fltLT06W5c+dKoaGhUl1dnaXNG2+8Ibm6ukqbN2+WTp06JS1cuFAKDAyUampqBFZuu44cOSKFh4dLI0eOlJ577jnLdd7n7quoqJDCwsKkRx99VDp8+LCUk5Mj7dixQzp//rylDe9zz/h//+//Sd7e3tL3338v5eTkSF9++aXk4uIivfvuu5Y2vNedl5SUJK1cuVLavHmzBED6+uuvW32/I/c0ISFBCg4OlpKTk6Xjx49LN998szRq1ChJr9f3eL0MQH1k/PjxUkJCQqtrUVFR0vPPPy+oov6npKREAiDt3btXkiRJMhqNUkBAgPTGG29Y2jQ1NUnu7u7SunXrRJVps2pra6XBgwdLycnJ0vTp0y0BiPe5ZyxfvlyaMmXKNb/P+9xz5s6dKz3++OOtrt19993SQw89JEkS73VP+GUA6sg9raqqkhwcHKQvvvjC0ubSpUuSXC6Xfvzxxx6vkUNgfUCn0yE1NRXx8fGtrsfHx+PQoUOCqup/qqurAQBeXl4AgJycHBQXF7e672q1GtOnT+d974Lf/e53mDt3Lm699dZW13mfe8bWrVsRFxeH//mf/4Gfnx/GjBmDDz74wPJ93ueeM2XKFOzcuRPnzp0DAJw4cQIHDhzAnDlzAPBe94aO3NPU1FQ0Nze3ahMUFISYmJheue88DLUPlJWVwWAwwN/fv9V1f39/FBcXC6qqf5EkCUuXLsWUKVMQExMDAJZ72959z8vL6/MabdkXX3yB48eP4+jRo22+x/vcMy5evIi1a9di6dKlWLFiBY4cOYJnn30WarUaixYt4n3uQcuXL0d1dTWioqKgUChgMBjw2muv4f777wfAf9O9oSP3tLi4GCqVCp6enm3a9MZ7JQNQH5LJZK3+LElSm2vUNc888wxOnjyJAwcOtPke73v3FBQU4LnnnsP27duh0Wiu2Y73uXuMRiPi4uLw+uuvAwDGjBmDjIwMrF27FosWLbK0433uvk2bNmHjxo347LPPMHz4cKSnp2Px4sUICgrCI488YmnHe93zunJPe+u+cwisD/j4+EChULRJsCUlJW3SMHXe73//e2zduhW7d+/GgAEDLNcDAgIAgPe9m1JTU1FSUoLY2FgolUoolUrs3bsXq1evhlKptNxL3ufuCQwMxLBhw1pdi46ORn5+PgD+e+5Jf/zjH/H888/jvvvuw4gRI/Dwww9jyZIlSExMBMB73Rs6ck8DAgKg0+lQWVl5zTY9iQGoD6hUKsTGxiI5ObnV9eTkZEyaNElQVbZPkiQ888wz2LJlC3bt2oWIiIhW34+IiEBAQECr+67T6bB3717e906YOXMmTp06hfT0dMtXXFwcHnzwQaSnpyMyMpL3uQdMnjy5zTYO586dQ1hYGAD+e+5JDQ0NkMtbv/0pFArLMnje657XkXsaGxsLBweHVm2Kiopw+vTp3rnvPT6tmtplXga/fv166cyZM9LixYslZ2dnKTc3V3RpNuu3v/2t5O7uLu3Zs0cqKiqyfDU0NFjavPHGG5K7u7u0ZcsW6dSpU9L999/Ppaw94OerwCSJ97knHDlyRFIqldJrr70mZWdnS//+978lJycnaePGjZY2vM8945FHHpGCg4Mty+C3bNki+fj4SP/7v/9racN73Xm1tbVSWlqalJaWJgGQ3nnnHSktLc2y3UtH7mlCQoI0YMAAaceOHdLx48elW265hcvg+4P3339fCgsLk1QqlTR27FjLcm3qGgDtfn300UeWNkajUXrppZekgIAASa1WS9OmTZNOnTolruh+4pcBiPe5Z3z33XdSTEyMpFarpaioKOmf//xnq+/zPveMmpoa6bnnnpNCQ0MljUYjRUZGSitXrpS0Wq2lDe915+3evbvd1+RHHnlEkqSO3dPGxkbpmWeekby8vCRHR0fpjjvukPLz83ulXpkkSVLP9ysRERERWS/OASIiIiK7wwBEREREdocBiIiIiOwOAxARERHZHQYgIiIisjsMQERERGR3GICIiIjI7jAAERERkd1hACIiuzJjxgwsXrxYdBlEJBgDEBEREdkdHoVBRHbj0UcfxSeffNLqWk5ODsLDw8UURETCMAARkd2orq7G7NmzERMTg1WrVgEAfH19oVAoBFdGRH1NKboAIqK+4u7uDpVKBScnJwQEBIguh4gE4hwgIiIisjsMQERERGR3GICIyK6oVCoYDAbRZRCRYAxARGRXwsPDcfjwYeTm5qKsrAxGo1F0SUQkAAMQEdmVZcuWQaFQYNiwYfD19UV+fr7okohIAC6DJyIiIrvDHiAiIiKyOwxAREREZHcYgIiIiMjuMAARERGR3WEAIiIiIrvDAERERER2hwGIiIiI7A4DEBEREdkdBiAiIiKyOwxAREREZHcYgIiIiMjuMAARERGR3fn/Xv7uP01xqEwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set up a new solver, now with more angular momenta and\n",
    "# a time dependent laser pulse.\n",
    "n_m = 16\n",
    "solver = CylinderFDM(r_max = 20, z_max = 20, n_r = n_r, n_z = n_z , n_m = n_m)\n",
    "tt, rr, zz = solver.get_trz_meshgrid()\n",
    "xx = rr*np.cos(tt)\n",
    "yy = rr*np.sin(tt)\n",
    "V = erfgau.potential_radial((rr*rr + zz*zz)**.5)\n",
    "D = xx # yes we could put the potential along the z axis to conserve symmetry, but this is just an example.\n",
    "\n",
    "solver.set_realspace_potential(V)\n",
    "solver.set_td_potential(D)\n",
    "laser = LaserPulse(t0 = 0, omega = 0.057, T = 100, E0 = 0.05)\n",
    "solver.set_td_potential_modulator(laser)\n",
    "\n",
    "dt = 0.01\n",
    "solver.setup_splitting_scheme(dt)\n",
    "t_range = np.arange(0, 100+dt, dt)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(t_range, laser(t_range))\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('E(t)')\n",
    "plt.title('Laser pulse')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ic| solver.m_i: array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7., -8., -7., -6., -5., -4.,\n",
      "                       -3., -2., -1.])\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t = 0.00, mu = -6.673423354066664e-17, norm = 1.0000000000000233\n",
      "t = 0.10, mu = 3.9014422629988084e-10, norm = 1.000000000000024\n",
      "t = 0.20, mu = 6.260728699532239e-09, norm = 1.0000000000000242\n",
      "t = 0.30, mu = 3.165384065510127e-08, norm = 1.0000000000000213\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[8], line 14\u001b[0m\n\u001b[1;32m     11\u001b[0m mu_hist[k] \u001b[38;5;241m=\u001b[39m mu\n\u001b[1;32m     13\u001b[0m \u001b[38;5;66;03m# propagate\u001b[39;00m\n\u001b[0;32m---> 14\u001b[0m psi \u001b[38;5;241m=\u001b[39m \u001b[43msolver\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpropagate_crank_nicolson\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpsi\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mt\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     16\u001b[0m \u001b[38;5;66;03m# check norm conservation\u001b[39;00m\n\u001b[1;32m     17\u001b[0m norm \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mlinalg\u001b[38;5;241m.\u001b[39mnorm(psi)\n",
      "File \u001b[0;32m~/Code/fdm_tdse/cylinder_fdm_3d.py:517\u001b[0m, in \u001b[0;36mCylinderFDM.propagate_crank_nicolson\u001b[0;34m(self, psi, t)\u001b[0m\n\u001b[1;32m    515\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m i_m \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mn_m):\n\u001b[1;32m    516\u001b[0m     temp \u001b[38;5;241m=\u001b[39m psi[i_m, \u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m] \u001b[38;5;241m-\u001b[39m \u001b[38;5;241m0.5\u001b[39mj\u001b[38;5;241m*\u001b[39mdt\u001b[38;5;241m*\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mT_m[i_m] \u001b[38;5;241m@\u001b[39m psi[i_m, \u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m]\n\u001b[0;32m--> 517\u001b[0m     psi[i_m,\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mT_rho_lu\u001b[49m\u001b[43m[\u001b[49m\u001b[43mi_m\u001b[49m\u001b[43m]\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msolve\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtemp\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    519\u001b[0m \u001b[38;5;66;03m# Kinetic step, z direction\u001b[39;00m\n\u001b[1;32m    520\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m i_m \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mn_m):\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# since the ground state is independent of theta, all the components are identical.\n",
    "psi = np.zeros(solver.shape, dtype=complex)\n",
    "psi[0, ...] = psi_ground_state[0, ...]\n",
    "ic(solver.m_i)\n",
    "\n",
    "mu_hist = np.zeros(len(t_range), dtype=float)\n",
    "\n",
    "for k,t in enumerate(t_range):\n",
    "    # compute dipole moment\n",
    "    mu = np.sum(psi.conj() * np.fft.fft(D * np.fft.ifft(psi, axis=0), axis=0)).real\n",
    "    mu_hist[k] = mu\n",
    "    \n",
    "    # propagate\n",
    "    psi = solver.propagate_crank_nicolson(psi, t)\n",
    "\n",
    "    # check norm conservation\n",
    "    norm = np.linalg.norm(psi)\n",
    "    \n",
    "    if k % 10 == 0:\n",
    "        print(f't = {t:.2f}, mu = {mu}, norm = {norm}')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ic| 'Computing sparse matrix'\n",
      "ic| V_m_max: 0\n",
      "ic| m: 0, m_ind: 0\n",
      "ic| diagonals.shape: (3,)\n",
      "ic| data.shape: (3, 90000)\n",
      "    diagonals.shape: (3,)\n",
      "    diagonals: array([ 90000,      0, -90000])\n",
      "ic| 'Setting total Hamiltonian as self.H_tot'\n",
      "ic| hasattr(solver, 'H_kin'): True\n",
      "ic| hasattr(solver, 'H_pot'): True\n",
      "ic| hasattr(solver, 'H_pot_td'): False\n",
      "ic| time_taken_sparse: 0.10306525230407715\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.10306525230407715"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "ic('Computing sparse matrix')\n",
    "\n",
    "start = time()\n",
    "H_mat_sparse = solver.get_sparse_matrix_fast()\n",
    "ic(hasattr(solver, 'H_kin'))\n",
    "ic(hasattr(solver, 'H_pot'))\n",
    "ic(hasattr(solver, 'H_pot_td'))\n",
    "\n",
    "time_taken_sparse = time() - start\n",
    "ic(time_taken_sparse)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ic| np.sort(E): array([-0.49958419, -0.12482951, -0.12479072, -0.05453834, -0.05379616,\n",
      "                       -0.05220115, -0.02305298, -0.01257923, -0.00961363, -0.00445813])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-0.49958419, -0.12482951, -0.12479072, -0.05453834, -0.05379616,\n",
       "       -0.05220115, -0.02305298, -0.01257923, -0.00961363, -0.00445813])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E, U = eigsh(H_mat_sparse, k=10, sigma=-0.6)\n",
    "ic(np.sort(E))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}