Generalized_Mie_Near_Field_Plots.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generalized Mie Near Field Plots\n",
    "\n",
    "Chad Heaps\n",
    "February 25, 2017\n",
    "\n",
    "This was originally a separate utility, but is being added to mstm-scan to aid anyone running the calculations.  The saved scattered coefficient file should be the correct format to be read into the notebook.\n",
    "\n",
    "This jupyter notebook plots the scattered near field from a Generalized Mie calculation.  The vector components are calculated using Fortran code from Daniel Mackowski's MSTM code (http://www.eng.auburn.edu/~dmckwski/scatcodes/).  The Fortran code is compiled using fortranmagic.  You could easily compile it with f2py separately, too.  There is a Fortran module func_nf and a python class NearField that take care of most of the computation.  They are mostly stitched together from other scripts/codes.\n",
    "\n",
    "Only Fields in the Cartesian planes are currently supported.  \n",
    "\n",
    "The code requires a text file with vector spherical harmonic coefficients and some header information.  Right now it will only calculate fields external to the particles.\n",
    "\n",
    "The field for each particle is calculated individually and summed at the end.  The field is  expressed as\n",
    "\n",
    "\\begin{equation}\n",
    "{\\bf E}_\\text{cart} = \\sum_{l=1}^N \\sum_{m=-n}^n a_{lm} {\\bf N}_{lm}(r,\\theta,\\phi) + b_{lm} {\\bf M}_{lm}(r,\\theta,\\phi)\n",
    "\\end{equation}\n",
    "\n",
    "where the explicit forms for ${\\bf N}_{lm}$ and ${\\bf M}_{lm}$ may be found in the MSTM manual or other Mie theory resources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-26T17:33:11.533495Z",
     "start_time": "2018-03-26T17:33:11.158908Z"
    },
    "hideCode": false,
    "hidePrompt": false
   },
   "outputs": [],
   "source": [
    "import time\n",
    "import re\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-26T17:33:11.739879Z",
     "start_time": "2018-03-26T17:33:11.535716Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "$.getScript(\"https://raw.github.com/marijnh/CodeMirror/master/mode/fortran/fortran.js\", function () {\n",
       "IPython.config.cell_magic_highlight['magic_fortran'] = {'reg':[/^%%fortran/]};});\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "try: \n",
    "    %load_ext fortranmagic\n",
    "except:\n",
    "    %install_ext https://raw.github.com/mgaitan/fortran_magic/master/fortranmagic.py\n",
    "    %load_ext fortranmagic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-26T17:33:15.098620Z",
     "start_time": "2018-03-26T17:33:11.743697Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Ok. The following fortran objects are ready to use: func_nf\n"
     ]
    }
   ],
   "source": [
    "%%fortran -v\n",
    "module func_nf\n",
    "    implicit none\n",
    "    contains\n",
    "    subroutine calcconstants (nbc, fnr, bcof) \n",
    "        integer, intent(in) :: nbc\n",
    "        real*8, intent(inout), dimension(0:2*nbc) :: fnr\n",
    "        real*8, intent(inout), dimension(0:nbc, 0:nbc) :: bcof\n",
    "        \n",
    "        integer l, n\n",
    "        do n=1,2*nbc\n",
    "            fnr(n)=dsqrt(dble(n))\n",
    "        enddo\n",
    "        bcof(0,0)=1.d0\n",
    "        do n=0,nbc-1\n",
    "            do l=n+1,nbc\n",
    "                bcof(n,l)=fnr(n+l)*bcof(n,l-1)/fnr(l)\n",
    "                bcof(l,n)=bcof(n,l)\n",
    "            enddo\n",
    "            bcof(n+1,n+1)=fnr(n+n+2)*fnr(n+n+1)*bcof(n,n)/fnr(n+1)/fnr(n+1)\n",
    "        enddo\n",
    "    end subroutine calcconstants\n",
    "\n",
    "!\n",
    "!  tau are the vector spherical harmonic functions, normalized\n",
    "!\n",
    "!\n",
    "!  last revised: 15 January 2011\n",
    "\n",
    "    subroutine taufunc(cb, phip,nmax,pvec)\n",
    "         implicit none \n",
    "         integer, intent(in) :: nmax\n",
    "         real(8), intent(in) :: cb, phip\n",
    "         complex(8),intent(out) :: pvec(nmax*(nmax+2),2)\n",
    "         \n",
    "         real(8) :: drot(-1:1,0:nmax*(nmax+2)),tau(0:nmax+1,nmax,2),fnm\n",
    "         integer :: n,m,nn1,mn\n",
    "         complex(8) :: ci\n",
    "\n",
    "         ci = dcmplx(0.0, 1.0)\n",
    "\n",
    "         call rotcoef(cb,1,nmax,drot)\n",
    "\n",
    "         do n=1,nmax\n",
    "            nn1=n*(n+1)\n",
    "            fnm=sqrt(dble(n+n+1)/2.d0)/4.d0\n",
    "            do m=-n,-1\n",
    "               mn=nn1+m\n",
    "               tau(n+1,-m,1)=-fnm*(-drot(-1,mn)+drot(1,mn))\n",
    "               tau(n+1,-m,2)=-fnm*(drot(-1,mn)+drot(1,mn))\n",
    "            enddo\n",
    "            do m=0,n\n",
    "               mn=nn1+m\n",
    "               tau(m,n,1)=-fnm*(-drot(-1,mn)+drot(1,mn))\n",
    "               tau(m,n,2)=-fnm*(drot(-1,mn)+drot(1,mn))\n",
    "            enddo\n",
    "         enddo\n",
    "         !!Reshaping to py preferred indexing\n",
    "         do n=1,nmax\n",
    "            nn1=n*(n+1)\n",
    "               do m=-n,-1\n",
    "                  mn=nn1+m\n",
    "                  pvec(mn,1) = 2.0*tau(n+1,-m,1)*exp(ci*m*phip)\n",
    "                  pvec(mn,2) = 2.0*tau(n+1,-m,2)*exp(ci*m*phip)\n",
    "               enddo\n",
    "               do m=0,n\n",
    "                  mn=nn1+m\n",
    "                  pvec(mn,1) = 2.0*tau(m,n,1)*exp(ci*m*phip)\n",
    "                  pvec(mn,2) = 2.0*tau(m,n,2)*exp(ci*m*phip)\n",
    "               enddo\n",
    "         enddo\n",
    "         \n",
    "    end subroutine taufunc\n",
    "\n",
    "    subroutine evalf(cthp, rdpw, phip, norder, ibess, nmn)\n",
    "        complex*16, intent(out), dimension(2, norder*(norder+2),3) :: nmn\n",
    "        complex*16, intent(in) :: rdpw\n",
    "        real*8, intent(in) :: cthp, phip\n",
    "        integer, intent(in) :: norder, ibess\n",
    "\n",
    "        integer n, mn, nn1,m\n",
    "        real*8 fnm\n",
    "        !integer nbc\n",
    "        !real*8, dimension(0:2*(2*(2*norder)+4)) :: fnr\n",
    "        !real*8, dimension(0:2*(2*norder)+4, 0:2*(2*norder)+4) :: bcof\n",
    "        real(8) :: legen(0:0,0:norder*(norder+2))\n",
    "\n",
    "        complex*16 rhfp\n",
    "        complex*16 ci\n",
    "        complex*16 rhf(0:norder)\n",
    "        complex(8) :: tau(norder*(norder+2),2)\n",
    "\n",
    "        ci = dcmplx(0.0,1.0d0)\n",
    "        \n",
    "        !nbc = 2*(2*norder)+4 \n",
    "        !call calcconstants (nbc, fnr, bcof) \n",
    "\n",
    "        call taufunc(cthp, phip, norder, tau)\n",
    "\n",
    "        if(ibess.eq.1) then\n",
    "           call cricbessel(norder, rdpw, rhf)\n",
    "        elseif(ibess.eq.3) then\n",
    "           call crichankel(norder, rdpw, rhf)\n",
    "        endif \n",
    "        !Ricatti-Bessel (Hankel) function, so divide by kr\n",
    "        rhf(:) = rhf(:)/rdpw\n",
    "        \n",
    "        !!For legendre function for nmn(1,1,:)\n",
    "        call rotcoef(cthp,0,norder,legen)\n",
    "\n",
    "        !Evaluate VSH components\n",
    "        do n=1,norder\n",
    "           rhfp=rhf(n-1)-n*rhf(n)/rdpw\n",
    "           nn1=n*(n+1)\n",
    "           !Normalization that is only used for the N_r component\n",
    "           fnm = sqrt(dble((2*n+1))/(2.0*(n*(n+1))))\n",
    "           do m=-n,n\n",
    "              mn=nn1+m\n",
    "              !CWH 04-05-2017\n",
    "              !These are the correct VSH terms\n",
    "              !!N(r,theta,phi)\n",
    "              nmn(1, mn, 1) = fnm*(n*(n+1))*(rhf(n)/rdpw)*legen(0,mn)*exp(ci*m*phip)\n",
    "              nmn(1, mn, 2) = rhfp*tau(mn,1)\n",
    "              nmn(1, mn, 3) = ci*rhfp*tau(mn,2)\n",
    "               \n",
    "              !M(r,theta,phi)\n",
    "              nmn(2, mn, 1)= cmplx(0.0,0.0)\n",
    "              nmn(2, mn, 2)= ci*rhf(n)*tau(mn,2)\n",
    "              nmn(2, mn, 3)= -rhf(n)*tau(mn,1)\n",
    "\n",
    "           enddo\n",
    "        enddo \n",
    "        nmn = nmn*sqrt(2.0)\n",
    "    end subroutine evalf\n",
    "    \n",
    "    \n",
    "    subroutine loop_points(norder, ibess, tM, tE, &\n",
    "                           pt, pkr, pp, pEsphere, npoints,ndim)\n",
    "    \n",
    "    \n",
    "        integer, intent(in) :: ndim, npoints, norder\n",
    "        integer, dimension(npoints), intent(in) :: ibess\n",
    "        real*8, intent(in), dimension(npoints) :: pt\n",
    "        complex*16, intent(in), dimension(npoints) :: pkr\n",
    "        real*8, intent(in), dimension(npoints) :: pp\n",
    "        complex*16, intent(in), dimension(ndim) :: tM\n",
    "        complex*16, intent(in), dimension(ndim) :: tE\n",
    "        complex*16, intent(out), dimension(ndim, npoints, 3) :: pEsphere\n",
    "        \n",
    "        integer :: i, j, k\n",
    "        complex*16 :: nmn3(2, ndim, 3)\n",
    "                \n",
    "        pEsphere(:,:,:) = dcmplx(0.0, 0.0)\n",
    "        do i=1,npoints\n",
    "            nmn3(:,:,:) = 0.0\n",
    "            if(ibess(i).eq.0) then\n",
    "                nmn3(:,:,:) = 0.0\n",
    "            else\n",
    "                call evalf(cos(pt(i)), pkr(i), pp(i), norder, ibess(i), nmn3)\n",
    "            endif\n",
    "            !subroutine evalf(cthp, rdpw, phip, norder, ibess, nmn)\n",
    "\n",
    "            do j=1,ndim\n",
    "                do k=1,3\n",
    "                    pEsphere(j, i, k) = tM(j) * nmn3(1, j, k) + & \n",
    "                                        tE(j) * nmn3(2, j, k)\n",
    "                enddo\n",
    "            enddo\n",
    "        enddo\n",
    " \n",
    "\n",
    "    end subroutine loop_points\n",
    "    \n",
    "    \n",
    "    \n",
    "    subroutine rotcoef(cbe,kmax,nmax,dc)\n",
    "         implicit none\n",
    "         \n",
    "         integer :: kmax,nmax,k,m,in,n,knmax,nn1,kn,im,m1, nbc\n",
    "         real(8) :: cbe,sbe,dc(-kmax:kmax,0:nmax*(nmax+2)),cbe2,sbe2,dk0(-nmax-1:nmax+1),&\n",
    "                    dk01(-nmax-1:nmax+1),sben,dkt,fmn,dkm0,dkm1,dkn1\n",
    "        real*8, dimension(0:2*(2*(2*nmax)+4)) :: fnr\n",
    "        real*8, dimension(0:2*(2*nmax)+4, 0:2*(2*nmax)+4) :: bcof\n",
    "\n",
    "        nbc = 2*(2*nmax)+4 \n",
    "\n",
    "\n",
    "        call calcconstants (nbc, fnr, bcof) \n",
    "\n",
    "\n",
    "         sbe=dsqrt((1.d0+cbe)*(1.d0-cbe))\n",
    "         cbe2=.5d0*(1.d0+cbe)\n",
    "         sbe2=.5d0*(1.d0-cbe)\n",
    "         in=1\n",
    "         dk0(0)=1.d0\n",
    "         sben=1.d0\n",
    "         dc(0,0)=1.d0\n",
    "         dk01(0)=0.d0\n",
    "         do n=1,nmax\n",
    "            knmax=min(n,kmax)\n",
    "            nn1=n*(n+1)\n",
    "            in=-in\n",
    "            sben=sben*sbe/2.d0\n",
    "            dk0(n)=in*sben*bcof(n,n)\n",
    "            dk0(-n)=in*dk0(n)\n",
    "            dk01(n)=0.d0\n",
    "            dk01(-n)=0.d0\n",
    "            dc(0,nn1+n)=dk0(n)\n",
    "            dc(0,nn1-n)=dk0(-n)\n",
    "            do k=-n+1,n-1\n",
    "               kn=nn1+k\n",
    "               dkt=dk01(k)\n",
    "               dk01(k)=dk0(k)\n",
    "               dk0(k)=(cbe*dble(n+n-1)*dk01(k)-fnr(n-k-1)*fnr(n+k-1)*dkt)&\n",
    "                     /(fnr(n+k)*fnr(n-k))\n",
    "               dc(0,kn)=dk0(k)\n",
    "            enddo\n",
    "            im=1\n",
    "            do m=1,knmax\n",
    "               im=-im\n",
    "               fmn=1.d0/fnr(n-m+1)/fnr(n+m)\n",
    "               m1=m-1\n",
    "               dkm0=0.d0\n",
    "               do k=-n,n\n",
    "                  kn=nn1+k\n",
    "                  dkm1=dkm0\n",
    "                  dkm0=dc(m1,kn)\n",
    "                  if(k.eq.n) then\n",
    "                     dkn1=0.d0\n",
    "                  else\n",
    "                     dkn1=dc(m1,kn+1)\n",
    "                  endif\n",
    "                  dc(m,kn)=(fnr(n+k)*fnr(n-k+1)*cbe2*dkm1 &\n",
    "                          -fnr(n-k)*fnr(n+k+1)*sbe2*dkn1  &\n",
    "                          -dble(k)*sbe*dc(m1,kn))*fmn\n",
    "                  dc(-m,nn1-k)=dc(m,kn)*(-1)**(k)*im\n",
    "               enddo\n",
    "            enddo\n",
    "         enddo\n",
    "         end subroutine rotcoef\n",
    "\n",
    "\n",
    "    real(8) function envj(n,x)\n",
    "        implicit none\n",
    "        integer :: n\n",
    "        real(8) :: x\n",
    "        n=max(1,abs(n))\n",
    "        envj=0.5d0*dlog10(6.28d0*n)-n*dlog10(1.36d0*x/n)\n",
    "    end function envj\n",
    "!\n",
    "!     ===================================================\n",
    "!     Purpose: Determine the starting point for backward\n",
    "!              recurrence such that the magnitude of\n",
    "!              Jn(x) at that point is about 10^(-MP)\n",
    "!     Input :  x     --- Argument of Jn(x)\n",
    "!              MP    --- Value of magnitude\n",
    "!     Output:  MSTA1 --- Starting point\n",
    "!     ===================================================\n",
    "!\n",
    "!\n",
    "!  last revised: 15 January 2011\n",
    "!\n",
    "         integer function msta1(x,mp)\n",
    "         implicit none\n",
    "         integer :: mp,n0,n1,it,nn\n",
    "         real(8) :: x, a0,f1,f,f0\n",
    "         a0=dabs(x)\n",
    "         n0=int(1.1*a0)+1\n",
    "         f0=envj(n0,a0)-mp\n",
    "         n1=n0+5\n",
    "         f1=envj(n1,a0)-mp\n",
    "         do it=1,20\n",
    "            nn=n1-(n1-n0)/(1.0d0-f0/f1)\n",
    "            f=envj(nn,a0)-mp\n",
    "            if(abs(nn-n1).lt.1) exit\n",
    "            n0=n1\n",
    "            f0=f1\n",
    "            n1=nn\n",
    "            f1=f\n",
    "         enddo\n",
    "         msta1=nn\n",
    "    end function msta1\n",
    "!\n",
    "!     ===================================================\n",
    "!     Purpose: Determine the starting point for backward\n",
    "!              recurrence such that all Jn(x) has MP\n",
    "!              significant digits\n",
    "!     Input :  x  --- Argument of Jn(x)\n",
    "!              n  --- Order of Jn(x)\n",
    "!              MP --- Significant digit\n",
    "!     Output:  MSTA2 --- Starting point\n",
    "!     ===================================================\n",
    "!\n",
    "!\n",
    "!  last revised: 15 January 2011\n",
    "!\n",
    "         integer function msta2(x,n,mp)\n",
    "         implicit none\n",
    "         integer :: n,mp,n0,n1,it,nn\n",
    "         real(8) :: x,a0,hmp,ejn,obj,f0,f1,f\n",
    "         a0=dabs(x)\n",
    "         hmp=0.5d0*dble(mp)\n",
    "         ejn=envj(n,a0)\n",
    "         if (ejn.le.hmp) then\n",
    "            obj=mp\n",
    "            n0=int(1.1*a0)\n",
    "         else\n",
    "            obj=hmp+ejn\n",
    "            n0=n\n",
    "         endif\n",
    "         f0=envj(n0,a0)-obj\n",
    "         n1=n0+5\n",
    "         f1=envj(n1,a0)-obj\n",
    "         do it=1,20\n",
    "            nn=n1-(n1-n0)/(1.0d0-f0/f1)\n",
    "            f=envj(nn,a0)-obj\n",
    "            if (abs(nn-n1).lt.1) exit\n",
    "            n0=n1\n",
    "            f0=f1\n",
    "            n1=nn\n",
    "            f1=f\n",
    "         enddo\n",
    "         msta2=nn+10\n",
    "        end function msta2\n",
    "        \n",
    "subroutine cspherebessel(n,z,csj,csy)\n",
    "         integer :: n,nm,k,m\n",
    "         real(8) :: a0\n",
    "         complex(8) :: z,csj(0:n),csy(0:n),csa,csb,cs,cf0,cf1,cf\n",
    "         a0=cdabs(z)\n",
    "         nm=n\n",
    "         if (a0.lt.1.0d-60) then\n",
    "            csj=(0.d0,0.d0)\n",
    "            csy=(-1.d300,0.d0)\n",
    "            csy(0)=(1.d0,0.d0)\n",
    "            return\n",
    "         endif\n",
    "         csj=(0.d0,0.d0)\n",
    "         csj(0)=cdsin(z)/z\n",
    "         csj(1)=(csj(0)-cdcos(z))/z\n",
    "         if (n.ge.2) then\n",
    "            csa=csj(0)\n",
    "            csb=csj(1)\n",
    "            m=msta1(a0,200)\n",
    "            if (m.lt.n) then\n",
    "               nm=m\n",
    "            else\n",
    "               m=msta2(a0,n,15)\n",
    "            endif\n",
    "            cf0=0.0d0\n",
    "            cf1=1.0d0-100\n",
    "            do k=m,0,-1\n",
    "               cf=(2.0d0*k+3.0d0)*cf1/z-cf0\n",
    "               if (k.le.nm) csj(k)=cf\n",
    "               cf0=cf1\n",
    "               cf1=cf\n",
    "            enddo\n",
    "            if (cdabs(csa).gt.cdabs(csb)) cs=csa/cf\n",
    "            if (cdabs(csa).le.cdabs(csb)) cs=csb/cf0\n",
    "            do k=0,min(nm,n)\n",
    "               csj(k)=cs*csj(k)\n",
    "            enddo\n",
    "         endif\n",
    "         csy=(1.d200,0.d0)\n",
    "         csy(0)=-cdcos(z)/z\n",
    "         csy(1)=(csy(0)-cdsin(z))/z\n",
    "         do k=2,min(nm,n)\n",
    "            if (cdabs(csj(k-1)).gt.cdabs(csj(k-2))) then\n",
    "               csy(k)=(csj(k)*csy(k-1)-1.0d0/(z*z))/csj(k-1)\n",
    "            else\n",
    "               csy(k)=(csj(k)*csy(k-2)-(2.0d0*k-1.0d0)/z**3)/csj(k-2)\n",
    "            endif\n",
    "         enddo\n",
    "         end subroutine cspherebessel\n",
    "            \n",
    "subroutine crichankel(n,ds,xi)\n",
    "         implicit none\n",
    "         integer :: n,i\n",
    "         complex(8) :: ds,psi(0:n),chi(0:n),xi(0:n),ci,&\n",
    "                       psi0\n",
    "         ci = dcmplx(0.0,1.0d0)\n",
    "         xi(0)=-ci*cdexp(ci*ds)\n",
    "         psi0=cdsin(ds)\n",
    "         if(cdabs(xi(0))/cdabs(psi0).lt.1.d-6) then\n",
    "            xi(1)=-cdexp(ci*ds)*(ci+ds)/ds\n",
    "            do i=2,n\n",
    "               xi(i)=dble(i+i+1)/ds*xi(i-1)-xi(i-2)\n",
    "            enddo\n",
    "         else\n",
    "            call cspherebessel(n,ds,psi,chi)\n",
    "            do i=1,n\n",
    "               xi(i)=(psi(i)+ci*chi(i))*ds\n",
    "            enddo\n",
    "         endif\n",
    "end subroutine crichankel \n",
    "!\n",
    "!  ricatti-bessel function psi(n), complex argument\n",
    "!\n",
    "!\n",
    "!  last revised: 15 January 2011\n",
    "!\n",
    "         subroutine cricbessel(n,ds,psi)\n",
    "         implicit none\n",
    "         integer :: n,i\n",
    "         complex(8) :: ds,psi(0:n),chi(0:n)\n",
    "         call cspherebessel(n,ds,psi,chi)\n",
    "         do i=0,n\n",
    "            psi(i)=psi(i)*ds\n",
    "         enddo\n",
    "         return\n",
    "end subroutine cricbessel\n",
    "            \n",
    "         subroutine ricbessel(n,ds,eps,nmax,psi)\n",
    "         implicit none\n",
    "         integer :: n,nmax,ns,i\n",
    "         real(8) :: ds,dns,sn,psi(0:n),psit,ds2,sum,eps,err\n",
    "         if(int(ds).lt.n) then\n",
    "            ns=nint(ds+4.*(ds**.3333d0)+17)\n",
    "            ns=max(n+10,ns)\n",
    "            dns=0.d0\n",
    "            do i=ns-1,n,-1\n",
    "               sn=dble(i+1)/ds\n",
    "               dns=sn-1.d0/(dns+sn)\n",
    "            enddo\n",
    "            psi(n)=dns\n",
    "            psi(n-1)=dble(n)/ds-1.d0/(dns+dble(n)/ds)\n",
    "            do i=n-2,1,-1\n",
    "               sn=dble(i+1)/ds\n",
    "               psi(i)=sn-1.d0/(psi(i+1)+sn)\n",
    "            enddo\n",
    "            psit=dsin(ds)\n",
    "            psi(0)=psit\n",
    "            ds2=ds*ds\n",
    "            sum=psit*psit/ds2\n",
    "            do i=1,n\n",
    "               psit=psit/(dble(i)/ds+psi(i))\n",
    "               sum=sum+dble(i+i+1)*psit*psit/ds2\n",
    "               err=dabs(1.d0-sum)\n",
    "               psi(i)=psit\n",
    "               if(err.lt.eps) then\n",
    "                  nmax=i\n",
    "                  return\n",
    "               endif\n",
    "            enddo\n",
    "            nmax=n\n",
    "         else\n",
    "            psi(0)=dsin(ds)\n",
    "            psi(1)=psi(0)/ds-dcos(ds)\n",
    "            do i=1,n-1\n",
    "               sn=dble(i+i+1)/ds\n",
    "               psi(i+1)=sn*psi(i)-psi(i-1)\n",
    "            enddo\n",
    "            nmax=n\n",
    "         endif\n",
    "         end subroutine ricbessel\n",
    "\n",
    "            \n",
    "            \n",
    "end module func_nf\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-26T17:33:15.106110Z",
     "start_time": "2018-03-26T17:33:15.102061Z"
    }
   },
   "outputs": [],
   "source": [
    "# print(\"{}\".format(func_nf.__doc__))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-26T17:55:19.320644Z",
     "start_time": "2018-03-26T17:55:17.644131Z"
    },
    "hideCode": false,
    "hidePrompt": false
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "############################################################\n",
    "#                        Chad Heaps                        #\n",
    "#                        03-10-2017                        #\n",
    "#                                                          #\n",
    "#       Program to evaluate the near-field                 # \n",
    "#       electric field for the particle-centered           #\n",
    "#       VSH expansion from Mackowski's code SCSMFO         #\n",
    "#                                                          #\n",
    "#                                                          #\n",
    "#                                                          #      \n",
    "#                                                          #\n",
    "############################################################\n",
    "\n",
    "class Cluster(object):\n",
    "    def __init__(self):\n",
    "        self.npart = 1\n",
    "        #self.maxorder = 1\n",
    "        #self.wavelen =  0.532\n",
    "        #self.refmed  = 0.532\n",
    "        #self.coef = np.zeros([npart, maxorder*(maxorder+2), 2], dtype=np.complex)\n",
    "        #self.norder = np.zeros([npart], dtype=np.int)\n",
    "        #self.coord = np.zeros([npart,4])\n",
    "\n",
    "\n",
    "\n",
    "class NearField(object):\n",
    "    def __init__(self):\n",
    "        self.coef_filename = 'pan.dat'\n",
    "        self.output_file = 'test.out'\n",
    "        #size of the grid\n",
    "        self.nx = 25\n",
    "        self.ny = 25\n",
    "        self.nz = 1\n",
    "        self.xlim = np.array([-50, 50])\n",
    "        self.ylim = np.array([-50, 50])\n",
    "        self.zlim = np.array([0, 0])\n",
    "        self.plane = 'xy'\n",
    "        self.save_int = 1\n",
    "        self.iprint = 0\n",
    "\n",
    "    def run(self):\n",
    "        \"\"\"\n",
    "        Perform all of the necessary steps for the calculation\n",
    "        \"\"\"\n",
    "\n",
    "        cluster = self.read_coef_file(self.coef_filename)\n",
    "        \n",
    "        #Convert to microns\n",
    "        #cluster.coord/=1000.0 \n",
    "        #cluster.wavelen/=1000.0\n",
    "        #Scalar wave number for medium\n",
    "        self.wavenum_med  = (2.*np.pi*cluster.refmed)/cluster.wavelen \n",
    "        if self.iprint > 0:\n",
    "            with open(self.output, 'w') as outfile:\n",
    "                outfile.write('Number of '\n",
    "                              'particles: {}\\n'.format(cluster.npart))\n",
    "                outfile.write('Wavelength:'\n",
    "                              '{}\\n'.format(1000.*cluster.wavelen))\n",
    "                outfile.write('Maximum order of multipole expansions:'\n",
    "                              '{}\\n'.format(cluster.maxorder))\n",
    "\n",
    "\n",
    "        #Grid for the image\n",
    "        if self.plane == 'xy':\n",
    "            npoints = self.nx*self.ny\n",
    "            xx, yy = np.meshgrid(np.linspace(self.xlim[0], \n",
    "                                             self.xlim[1], \n",
    "                                             self.nx), \n",
    "                                 np.linspace(self.ylim[0], \n",
    "                                             self.ylim[1], \n",
    "                                             self.ny))\n",
    "            grid_x = xx.reshape(npoints)\n",
    "            grid_y = yy.reshape(npoints)\n",
    "            grid_z = np.zeros([npoints])\n",
    "            grid_z.fill(self.zlim[0])\n",
    "#             Take inputs in nm and convert to micron\n",
    "            grid_x/=1000.\n",
    "            grid_y/=1000.\n",
    "            grid_z/=1000. \n",
    "            plane_grid = np.vstack((grid_x, grid_y))\n",
    "\n",
    "        if self.plane == 'xz':\n",
    "            npoints = self.nx*self.nz\n",
    "            xx, zz = np.meshgrid(np.linspace(self.xlim[0], \n",
    "                                             self.xlim[1], \n",
    "                                             self.nx), \n",
    "                                 np.linspace(self.zlim[0], \n",
    "                                             self.zlim[1], \n",
    "                                             self.nz))\n",
    "            grid_x = xx.reshape(npoints)\n",
    "            grid_z = zz.reshape(npoints)\n",
    "            grid_y = np.zeros([npoints])\n",
    "            grid_y.fill(self.ylim[0])\n",
    "            #Take inputs in nm and convert to micron\n",
    "            grid_x/=1000.\n",
    "            grid_y/=1000.\n",
    "            grid_z/=1000. \n",
    "            plane_grid = np.vstack((grid_x, grid_z))\n",
    "\n",
    "        if self.plane == 'yz':\n",
    "            npoints = self.ny*self.nz\n",
    "            yy, zz = np.meshgrid(np.linspace(self.ylim[0], \n",
    "                                             self.ylim[1], \n",
    "                                             self.ny), \n",
    "                                 np.linspace(self.zlim[0], \n",
    "                                             self.zlim[1], \n",
    "                                             self.nz))\n",
    "            grid_y = yy.reshape(npoints)\n",
    "            grid_z = zz.reshape(npoints)\n",
    "            grid_x = np.zeros([npoints])\n",
    "            grid_x.fill(self.xlim[0])\n",
    "            #Take inputs in nm and convert to micron\n",
    "            grid_x/=1000.\n",
    "            grid_y/=1000.\n",
    "            grid_z/=1000.\n",
    "            plane_grid = np.vstack((grid_y, grid_z))\n",
    "        \n",
    "\n",
    "        if self.iprint > 0:\n",
    "            with open(self.output, 'w') as outfile:\n",
    "                outfile.write(\"Number of points in the image: {}\\n\".format(npoints))\n",
    "        #'part' prefix indicates array with coordinates for all particles\n",
    "        part_x = grid_x[:, None] - cluster.coord[:,1]\n",
    "        part_y = grid_y[:, None] - cluster.coord[:,2]\n",
    "        part_z = grid_z[:, None] - cluster.coord[:,3]\n",
    "        \n",
    "        part_r = np.sqrt(part_x**2 + part_y**2 + part_z**2)\n",
    "        #offset from zero slightly\n",
    "        part_r[np.where(part_r==0)] = 1e-10\n",
    "         \n",
    "        #Calc spherical coordinates\n",
    "        #Normalized to r for correct angles\n",
    "        part_phi =  np.angle((part_x/part_r)+1.0j*(part_y/part_r))\n",
    "        part_theta = np.arccos(part_z/part_r)\n",
    "      \n",
    "        #Determine which bessel function to use\n",
    "        rless  = part_r < cluster.coord[:,0]\n",
    "        rgreat = part_r > cluster.coord[:,0]\n",
    "        rlany = (1 - np.sum(rless,axis = -1))\n",
    "        \n",
    "        #the integer 1 and 3 correspond to ibess values\n",
    "        ibess = 0*rless + rlany[:,None]*3*rgreat\n",
    "        #ibess = self.intfield*1*rless + rlany[:,None]*self.extfield*3*rgreat\n",
    "        part_kr = 0.0*part_r*rless + self.wavenum_med*part_r*rgreat\n",
    "        \n",
    "        if self.iprint > 0:\n",
    "            with open(self.output, 'w') as outfile:\n",
    "                outfile.write('wavenum_med {}\\n'.format(self.wavenum_med))\n",
    "\n",
    "        maxdim = cluster.maxorder*(cluster.maxorder+2)\n",
    "        #Separated into TM and TE parts\n",
    "        empEspher = np.zeros([cluster.npart, 2, maxdim, npoints, 3], dtype=np.complex)\n",
    "        #Multiplied by coefficients removing one index\n",
    "        pEdat = np.zeros([cluster.npart, maxdim, npoints, 3], dtype=np.complex)\n",
    "        self.pEspher = np.zeros([cluster.npart, npoints, 3], dtype=np.complex)\n",
    "        self.pEcart  = np.zeros([cluster.npart, npoints, 3], dtype=np.complex)\n",
    "        t_startf = time.time() \n",
    "        for ipart in range(cluster.npart):\n",
    "            #Calculate field for each particle center in spherical coordinates\n",
    "            pkr =  part_kr[:,ipart]\n",
    "            pt =  part_theta[:,ipart]\n",
    "            pp =  part_phi[:,ipart]\n",
    "            nmax = cluster.norder[ipart]\n",
    "            ndim = nmax*(nmax+2)\n",
    "            tM = cluster.coef[ipart, :ndim, 0]\n",
    "            tE = cluster.coef[ipart, :ndim, 1]\n",
    "            pEdat[ipart, :ndim, :, :] = func_nf.loop_points(nmax, ibess[:,ipart],\n",
    "                                                              tM, tE, pt, pkr,\n",
    "                                                              pp)\n",
    "            #print(pEdat)\n",
    "            self.pEspher[ipart, :, :] = pEdat[ipart,:,:,:].sum(axis=0)\n",
    "            self.pEcart[ipart, :, :] = self.vec_spher_to_cart(self.pEspher[ipart, :, :], pt, pp)\n",
    "            if self.iprint > 0:\n",
    "                with open(self.output, 'w') as outfile:\n",
    "                    outfile.write('time to '\n",
    "                    'evaluate functions {}\\n'.format(time.time() - t_startf))\n",
    "        \n",
    "        Esphtot =  self.pEspher.sum(axis=0)  \n",
    "        Ecarttot = self.pEcart.sum(axis=0)\n",
    "        \n",
    "        self.Ecartint = (np.abs(Ecarttot)**2).sum(axis=-1)\n",
    "        self.Esphint = (np.abs(Esphtot)**2).sum(axis=-1)\n",
    "        self.Ecarttot = Ecarttot\n",
    "        self.grid_x = 1000.*grid_x\n",
    "        self.grid_y = 1000.*grid_y\n",
    "        self.grid_z = 1000.*grid_z\n",
    "        if self.save_int:\n",
    "            np.savetxt(self.name + '.int', np.vstack((plane_grid, self.Ecartint)).T) \n",
    "        if self.save_etot:\n",
    "            np.savetxt(self.name + '.etot', np.vstack((plane_grid, \n",
    "                                           self.Ecarttot[:, 0].real, \n",
    "                                           self.Ecarttot[:, 0].imag, \n",
    "                                           self.Ecarttot[:, 1].real, \n",
    "                                           self.Ecarttot[:,1].imag, \n",
    "                                           self.Ecarttot[:, 2].real, \n",
    "                                           self.Ecarttot[:,2].imag)).T)\n",
    "\n",
    "   \n",
    "        return\n",
    "\n",
    "    def vec_spher_to_cart(self, spv, theta, phi):\n",
    "        \"\"\" \n",
    "        Takes spherical vector components spv and expresses them in cartesian components\n",
    "        returns cartv\n",
    "        \"\"\"\n",
    "        if len(spv.shape) == 1:\n",
    "            spv = spv.reshape(-1,1)\n",
    "        cartv = np.zeros(spv.shape, dtype=np.complex)\n",
    "        st = np.sin(theta)\n",
    "        ct = np.cos(theta)\n",
    "        sph = np.sin(phi)\n",
    "        cph = np.cos(phi)\n",
    "        cartv[:,0] = ( spv[:,0] * st * cph \n",
    "                     + spv[:, 1] * ct * cph\n",
    "                     - spv[:, 2] * sph )\n",
    "        cartv[:,1] = ( spv[:,0] * st * sph \n",
    "                     + spv[:, 1] * ct * sph\n",
    "                     + spv[:, 2] * cph ) \n",
    "        cartv[:,2] = ( spv[:,0] * ct\n",
    "                     - spv[:, 1] * st )\n",
    "        return cartv\n",
    "    \n",
    "    def read_coef_file(self, file):\n",
    "        ##Extract information from coefficient file\n",
    "        find_float = re.compile('-?[0-9]*\\.?[0-9]+')\n",
    "        find_int = re.compile('-?[0-9]+\\s')\n",
    "        find_sci = re.compile('-?[0-9]*\\.[0-9]*E[+-]?[0-9]+')\n",
    "        coef_file = open(file, 'r')\n",
    "        coef_lines = coef_file.readlines()\n",
    "        coef_file.close()\n",
    "\n",
    "        cluster = Cluster() \n",
    "\n",
    "        prelim  = re.split('\\s+', coef_lines[0].strip())\n",
    "        cluster.npart = np.int(prelim[0])\n",
    "\n",
    "        cluster.maxorder = np.int(prelim[1])    \n",
    "        cluster.wavelen = np.float(prelim[2])\n",
    "        cluster.refmed  = np.float(prelim[3]) + 1.0j*np.float(prelim[4])\n",
    "        if cluster.refmed.imag != 0.0:\n",
    "            print('Warning, a complex refractive index medium was used in '\n",
    "                  'the scattering calculation.  This has not been tested '\n",
    "                  'for the near field code!')\n",
    "        \n",
    "        cluster.coef = np.zeros([cluster.npart, cluster.maxorder*(cluster.maxorder+2), 2], dtype=np.complex)\n",
    "        \n",
    "        cluster.norder = np.zeros([cluster.npart], dtype=np.int)\n",
    "        cluster.coord = np.zeros([cluster.npart,4])\n",
    "        i_line = 1\n",
    "        for i_part in range(cluster.npart):\n",
    "            part_prelim = re.split('\\s+', coef_lines[i_line].strip())\n",
    "            cluster.norder[i_part] = np.int(part_prelim[0])\n",
    "            cluster.coord[i_part,] = np.asarray(part_prelim[1:5])\n",
    "            i_line+=1\n",
    "            for i in range(cluster.norder[i_part]*(cluster.norder[i_part]+2)):\n",
    "                tc = find_sci.findall(coef_lines[i_line])\n",
    "                cluster.coef[i_part,i,0] = np.float(tc[0]) + 1.0j*np.float(tc[1])\n",
    "                cluster.coef[i_part,i,1] = np.float(tc[2]) + 1.0j*np.float(tc[3])\n",
    "                i_line+=1\n",
    "\n",
    "        return cluster\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Heterodimer near-field scattering\n",
    "\n",
    "Example 1 plots the near field for a dimer where one particle is silver and the other is gold.  The plot uses the scattered field at 360 nm.  The result is excitation at the silver nanoparticle but not the gold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-26T18:07:01.121148Z",
     "start_time": "2018-03-26T18:07:00.592035Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAF3CAYAAABTxMQNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xu0ZGdZ5/HfU+f0JfcmtAmddDAXAggoBEMAGVxIECMyCc6AE0VsIa6oSxAXOiaAMzIucQUdRVwqTito0AhkIpheiIQkEi+juXQukEsDuRCSJh3aSK6k093n1DN/1D5V796n3qq3bnvX3uf7Weus3rVr195vnbr0e97neZ/X3F0AAAAYrlV1AwAAAOqCjhMAAEAiOk4AAACJ6DgBAAAkouMEAACQiI4TAABAoko7Tma2ycwuM7Mvm9kuM3u5mR1tZlea2Z3Zv0+rso0AAAArqh5x+pCkz7n7cyW9UNIuSRdKutrdT5V0dXYbAACgclZVAUwzO1LSFyWd7EEjzOwrkl7l7nvMbIuka9z9OZU0EgAAIFDliNPJkv5d0p+b2c1m9mdmdpikY919jyRl/x5TYRsBAAC6Fiu+9oslvcPdrzOzD2mEsJyZnS/pfEla0ML3HmpHJV/YRmxoaUpv2Nz+JiLmdHmgips10eVZcmltssTPPu+PNelxPfyQu39Hmdf8oR84zP/jW8sTnePGL+2/wt3PmlKToqrsOO2WtNvdr8tuX6ZOx+mbZrYlCNXt7fdgd98uabskHdl6ur9s8YeSL2ytKXYYbIqDdtNsV4rUL895Ma9f4u0J2+XtCR8+wfUnvHbaJeb0dSvBVL9rpin1e6uE90eKyt9Dc/J7KMtVftnXy77mQ99a1nVXbJ3oHOu23L15Ss0ZqLJQnbs/KOl+M1vJXzpT0h2Sdkjalu3bJunyCpoHAACwSpUjTpL0DkmXmNl6SfdIeqs6nblLzew8SfdJetO4J18TI0t1GzWaRNXPNTbiNc7rGf4FnfLeGvAXb+x9nvRX+qBrT+mv7NTPYeWjChOofGRpmt9Pk5x3RiMz4e+3kvdJ+HtYY6NP5XEt1+R3W2nHyd1vkXR6n7vOLLstAACgGi6pXXXCaKKqR5wAAADUVj1GnKougAkAAFAbjRtxmjjXYJJcgYmvPb08CWvRJ542b0/411CYIxV7r8TyN8bIRZoo92nQNUvIY8lfbj6H70vPa5pVHtM0lfCeId+pmVyu5XmdOV3QuI4TAACoH3KcAAAAErikZTpO5TFNMGw+zvB3ySG5icNuVU/jb4hJX4ekUF/KJYrhiRHDI00J4a2+/Gy+dJsSkpv0eUz0+51R2Yvic5rXsC7S1GXEqQZBcwAAgPnQiBEnAABQXy6RHD7XygjPJYbHRg7/TBp2Y7bd9CXOtkt5raPhvJQZeVJaRfLEEN7IYY+KZxtVXr17VFMMyZXx3Ee9xlih3wnfN6XPuGOG3VTV5Te4NjtOAABgbric5HAAAIAkLi3Xo99EcjgAAECqZo841S2XaZz8pbJzpHLXrllOyThS8iQWFnrbqcmNkVym2PskuWp57O0wRu7TRGULZjT9vPYmyGuqdFWEooTXcFB7o++hkktdTBX5ThPpLPJbD83uOAEAgBowLasef4wTqgMAAJVydQbGJ/kZxsw2mdllZvZlM9tlZi83s6PN7EozuzP792nDztO8EadRh6NTh78TQlzJpQVSwmVTPdfovXirw4Kiowp/D6lTlReGHxLyQUP04TVTwntBeM4Sw4G5kN6oiwoXX/OGVB6v1Bifo7FCcmV/XlOuN+D1LKNsQKWLARO2m1cfkvQ5d3+jma2XdKik90i62t0vMrMLJV0o6YJBJ2lexwkAANTOLEN1ZnakpO+X9NOS5O4HJB0ws3MkvSo77GJJ14iOEwAAmGedRX4n7jhtNrOdwe3t7r492z5Z0r9L+nMze6GkGyW9U9Kx7r5Hktx9j5kdM+wizek4jTJUXUZ4LnX22qiPT2x7UqhtrJBADZL3ksKXiecacQkA8wHXtv7nyoX3UsJ5A2bYhe/HkauQr7HFg2dmVuG51POWPds1ZYbcpGG7CUNfcxO26zSg3OvXSHvQ92eah9z99Mh9i5JeLOkd7n6dmX1InbDcyBqYyAIAAOpkZcRpkp8hdkva7e7XZbcvU6cj9U0z2yJJ2b97h52IjhMAAGg0d39Q0v1m9pxs15mS7pC0Q9K2bN82SZcPO1dzQnX9JA1/T1jAsuSQ3MAQXOz5prSxjOKbAyzteXBq5yrD4nFbejcSZsVJyv+ZEjwmF95LCecVX6tIGC82Ey8awhv0cq7hxYOTlB2em2K6wSDtJ5+c6PGzlnuP10Ud3s8VcJmWZz+W8w5Jl2Qz6u6R9FZ1vvkuNbPzJN0n6U3DTtLsjhMAAKiFKeQ4DeTut0jqlwN15ijnoeMEAAAqNaVZdaWg4wQAACpmWvZ6pF03r+NUZV7ToJyfhFymaP7SoOcUO++obUzNhRgjZ2LpGw+M/Jh5t/TAnqHHLB5/XPzOMC8pVmogMQ8qmv8Uy30aZyHhURcPTpx+nr9+8xcPjuY1jfrZH+M7bPmJJ5IeUye+vDzyY8K8qNJLExSR71RLzes4AQCAWnFJ7ZpM9KfjBAAAKkeOU5lMw0N0kaHtsRbmHSfcNWpILjUEl9KWUcN5xbZELN23e+gxa1lqiHJx6/G9G+HrGYbOwu1CqC4WxptWCE9KXDw4pWRB52R9d09c3bnssEdCCYKBJQcmCM8Neq2aGJKbpmh4r+qFzZtSOX9M7vXJcapHKwEAAOZAM0acAABArbUJ1c2BccJzk8xGGzAsnwvPjRqSC69dPD7p8SnXKzy/4DFL99zb/zGYiqXd3+i7f/GZW3s3liOhMkkKQw/WfyaeKzgmHPmPhQYLkhYPTpl51zlZb7tuYbtZhXNGDM8RjpuBQe+TKsN4a2TmXaeOUz2CYM3uOAEAgBqoT44THScAAFCpOpUjqEcrAQAA5kDzRpxGzWsaOCV/QG5R95iEMgOF46KlAmLXG1RaIHreCXKqJC3ddU//41CaWLmHXO6TlH/tIiUMcq9uSsmCzp39z5VSeTxWskAaudr43OY7RS+XWIIgIa+p/eSTU2oVJhJ735Sd+1TzyvnDLM94kd9paV7HCQAA1IrLSA4HAABI1SY5vEzWNwRVenhuUEhskpDcoPMGC1aOeq6lr9wl1E8xhLd44nf2brSDsgPLkcrjQfmCgYsHh2UOEsoWJJUskPKZlROE7TqHVbxI6zCJoZz2vn0zbghmYl5CeMVr1jBsV6dyBPVoJQAAwBxoyIgTAACoK5eRHF6VkcNzxeMTKoEnhecGLcYbC9XFzhULx606b//H+0Jve/n2rwrNsnTv1/vuX/zOE3o3liPvrUjVcSk/E2+SBYNXhe1y7/lg/4hhu85hI864GyecMelivqHgOMJzDZb7vBC2S1WXOk6N6zgBAIB6cVdtKofXo5UAAABzgBEnAABQMVNb5DiVqm9uU0pe06BK3LlTJeQ1DcpFyh0XuX7s8bE8Jkke3hfmMt32FWFtW/r6/d3taMmCWNXxgrBsgSt4fPiQSMmC4mczl/M0xXynOiCvaQ0qvmfnpdr4nH2WXITqkpnZgpndbGafyW6fZGbXmdmdZvZJM1tfdRsBAMBsLas10U9ZKu84SXqnpF3B7Q9I+qC7nyrpYUnnVdIqAABQCpep7ZP9lKXSUJ2ZbZX0I5LeL+ldZmaSXi3pJ7JDLpb0PkkfHuGk/fenlABYdaowDBZ7fML+4nVaC5H9/a8XC8cVH7P8pV0C+glLFiyeclLvDlvqbS8PCC+H1cZzx0QWDA4VQoBJiwSnhO06DwrumnAx4GkJ2tje/1R17cB8Ct+zQXpG6e/ZBoXAy1Z1jtPvS/pVSUdkt58u6RF3X/k23y3p+H4PNLPzJZ0vSRvtsBk3EwAAzBJLrgxhZq+XtNfdbwx39zm0bzfc3be7++nufvp62ziTNgIAgNlzdRb5neSnLFWOOL1C0tlm9jpJGyUdqc4I1CYzW8xGnbZKemDsK5QRnotdY6HwIsbOGwvPLfY/xguhuvbNdwgYxdLdX+tuLz7r5N4dtpw/sF24vSJ8z4YhvOBwD6LRq8LnYbVyj4QnUmbblWXEKtCE55AqXEg7DNvljinjPb8qBD77S/ZphJZrUo6gshEnd3+3u2919xMlnSvpH9z9zZK+IOmN2WHbJF1eURMBAEAJ6jTiNI8BxQvUSRS/S52cp49U3B4AAABJ1SeHS5Lc/RpJ12Tb90g6Y+STrAztDypu2fdhhb7jtMJzxWHXlBlzi8NDde0bbxMwLUt33dPdXjz1lPydsYWBlyIhvEC0YKYkKfhshKGKQQsDryiG2dv9ZwaVPcOufeDAzK+BZouF7WILSFc6c3RG6hKqm4uOEwAAWLvcrdRw2yToOAEAgMqx5AoAAEDDNGPEyWx1pe6iIE4cXbB35Vzd+ybIayq2ZzGoEBuWFAgrx0ZynPz6WwXM2tKdd+duLz7nWb0bsXzBWB5UrNK4ijlPI+Y7xcoXVMCXDlbdBDRUSpmCWO6TVM/8J5fUJscJAAAghdUmVEfHCQAAVKpTx4kRp2pEKoRHw3OrFuOdIDyX286fNxd6C48L9rfXBfuv/aKAKi195a7u9sJ3Pbu7beFnISwXvhQsGDxA+NWYC9uFf222w9ICA8oU5KqKz37RUsJzKFtK2K5obha8HhFr1QEAADRM80acAABArbiMUF3pupXDE8JzsZlzxfsmCM/lKn8rH57zdeF2GKrrbTMUiHmyvOur3e2F50fCduFnbKywXUSsinhJwlAJUKWmh+3aM/6fz8zulfS4pGVJS+5+upkdLemTkk6UdK+kH3P3hwedh/+fAQBApdylZbeJfhL9gLu/yN1Pz25fKOlqdz9V0tXZ7YHoOAEAgMq13Sb6GdM5ki7Oti+W9IZhD6DjBAAA1gKX9Hkzu9HMzs/2HevueyQp+/eYYSdpTo5TP7G8pti2lM95SqkKHslrCvOYpGKOU5DXtL63v/WPNxWfATB3fH3wtRF8fuzgZImdFuRfuAW5GMHnzQqVw/PlCaaXv0FeE+Zd8T2akvMUqzY+D7lPneTwicdyNpvZzuD2dnffHtx+hbs/YGbHSLrSzL48zkWa3XECAAC1sDz5kisPBblLq7j7A9m/e83s05LOkPRNM9vi7nvMbIukvcMuQqgOAABUaqVy+KxynMzsMDM7YmVb0msl3SZph6Rt2WHbJF0+rK3NGXHKhiBzJQhCsRIEq0J1YRgvUkV8MbIwbxiOKwybehCSy5UdIDyHmmnffEd3u/W9L+jdEYbtxjlxEIbLlSkIQxKxxYal/Gc3qCJurV44bx5CEsBMxP7vS6iivyqEV20FkFk5VtKnsxIqi5L+2t0/Z2Y3SLrUzM6TdJ+kNw07UXM6TgAAoKamkuMU5e73SHphn/3/IenMUc5FxwkAAFSuPXmOUyma0XGyIESXUiE8NnNOklqRhXrD2XMLQSXWhIrgktRe37u9HGyTZIY6a994W3c7F7YLJH8VeiSMFoYaCn+RhrPsRp5hF4Q2WLwXdRe+h21xXe+OMIQ3o8Wvp2GlAGYdNKPjBAAAam2WobppqkcrAQAA5gAjTgAAoFKdApiE6qqRUiF8UOXwhf65UGHZgVwV8Vy18P4VwaV8XtO6z+8U0DTLh/byKlrB56o4rB39aoyUCsjnMRWOCUsghJ/XMF8qUpqgfeBArCVArdU134nkcAAAgAQrBTDrgBwnAACARA0ZcbJeGYKUCuFhyYLiwoi58Fz/hX3DsF17Xf+K4GH5AYnwHJrP/vnm7nb7ladFjws/Gbm/L8MwnEcW+S1W5A9vxCqMhxds87ci1pawKngu1D2HYbu6zKprSMcJAADUVsJ6c/OCjhMAAKiUi+Tw8q0MzcfCc7mh+4X++6VCtfBgiDM3e65/hfBwJt3yhvyQY3N+0cBwS4f23vGD3vutcGHflAV4C9XFc4sBJywS3N6/b/g1gAYJZ4+21q/vbkfDdhWqy4jTfPy2AAAAaoCBEAAAUKk6lSNoTsepGHIr7ssVtozMvNOABXxz20F4LjfDrvfYDZ+5PqnZQBOtu+KG7vbB156euy/80rHIwr65/blFfovHL/UeE5w5DM/V46sYKFd0tl2F6DgBAAAkYMkVAACAEdRlVh3J4QAAAImaMeJk6uUqxfKawumWYV7TQqHv2BpejiBWIbxYggCAdPDwhfidQc7SgvevHJ7+N2iQ7xTkbyw98kjyGYAmi5UmCD8vlXFynAAAAJIwqw4AAGAEdJyqEq0W3j+EF5YfkCSPLOabqxCeC9X1Hn/I3143VpOBJjv0U/nPxb4ffWl328LixUHVgVjQe+DXavh5b8/HoqXA3JrDRX7ronkdJwAAUCuUIwAAABiB03Eqk/WG6XOhushMOotXDlcYnouE9MIK4cvr6vFCA/PiwOG9z1UYnlMQtlsMZtW1ItXFJclys2h71cKX7rl3kiYCjdfe/1R3u7VhY4Ut6alLHaeGdJwAAEBdeY3KEVB4CAAAIFFlHSczO8HMvmBmu8zsdjN7Z7b/aDO70szuzP59WlVtBAAA5XC3iX7KUmWobknSL7v7TWZ2hKQbzexKST8t6Wp3v8jMLpR0oaQLhp5tJdehFSlBEORChPlKXshx8pRq4WGO04Z6DC0C8+Kov/y37vZjP/ny3h1htXDv/9XUssLnLfj82kEG0IH6qs+susq+adx9j7vflG0/LmmXpOMlnSPp4uywiyW9oZoWAgCAsjDiNAIzO1HSaZKuk3Ssu++ROp0rMzsm8pjzJZ0vSRsXjiinoQAAYOpYcmUEZna4pL+R9Evu/pgVh+Ij3H27pO2SdNT6Y71/OYLIgr+RhXylASUIFsNyBL1jjvj4tUntBbDagcN7nytr9xYDtkgFgmKl/9aBIFS3SKgOGMs8LPJbI5V2nMxsnTqdpkvc/VPZ7m+a2ZZstGmLpL3VtRAAAMyc59Ic51qVs+pM0kck7XL33wvu2iFpW7a9TdLlZbcNAACUqy2b6KcsVY44vULSWyTdama3ZPveI+kiSZea2XmS7pP0pqFnMvWGGiPVwj0SwvPC8H7STLr1Q1sEIMGBI3vb1g4+ox5+dnu7F1vx0Hrr+lun3j5gLWjv21d1E+RiyZWh3P1fFF/s/Mwy2wIAAJCi8uRwAACw1tWnjlNzOk7DCmCGhS0HFcAMHtOOzKpbXl+PFxeYd8f99r92t/f88vf17vD+i3J7IVS3sL53HBF0oN7qkhzenI4TAACoLXKcAAAAErjXp+NExTgAAIBEzRtxilUOD7fDvKbi9OZF67vdDvKaNl38bwIwXQeO6m3nKodbmHeY/1tv4UDvQHKcgHqrS3I4I04AAKBy7pP9pDCzBTO72cw+k90+ycyuM7M7zeyTZjb0bzA6TgAAoHLuNtFPondK2hXc/oCkD7r7qZIelnTesBM0JFRnvYrh4dRl6192wCOlCaR4CYJwG8D0HdzU7m6b9/+brvh5PerDhM2BJnCN1PkZi5ltlfQjkt4v6V3Z0m+vlvQT2SEXS3qfpA8POg8jTgAAoAk2m9nO4Of8wv2/L+lXJa38lfZ0SY+4+1J2e7ek44ddpCEjTgAAoM6mUP/yIXc/vd8dZvZ6SXvd/UYze9XK7nGa0YyOU26R33D2nPpu58J2xVl1rf5hvHYzflPA3LKj93e3D7Y3BPcEi23zOQSaafZ1nF4h6Wwze52kjZKOVGcEapOZLWajTlslPTDsRITqAABA9XzCn0Gndn+3u2919xMlnSvpH9z9zZK+IOmN2WHbJF0+rJl0nAAAwFp1gTqJ4nepk/P0kWEPYOAbAABUrqwlV9z9GknXZNv3SDpjlMc3p+O0ktvUSqgcHmyuynEKSxCs6+1fXkc5AmCWTvnJm7vbd/5FL7/zoPc+iMVyBACaI7WIZdWa03ECAAC15KrPIr90nAAAQLVcEh2nikTCc6mVw3MlCILtp//pv06zlQAGOGLTk93tx9uHdrcPLjTvKwtAvfAtBAAAKkeOEwAAQCo6TiXLwnIem0nXWn2sJBXXEg1v+8IU2wcg2dZNj3S3713ufSj3LWzodziA2pv9Ir/T0pyOEwAAqK+ajDhRORwAACARI04AAKBas1/kd2oa0nGyoHJ4b28s3ymsFr6qcngrLEcw3VYCSPOswx/qbu872Kscvu6/3lZFcwCUoSahuoZ0nAAAQL3VY8SJHCcAAIBEzR5xiobtgmOKlcPDx7Tq0fsFmua7Dnugu/3IwUO62w9W0RgA5SBUBwAAkIiOEwAAQAIW+S2ZSVoJq8UqhwebYdiu+DqF1cKpHA5U47kbeqG6PYdu6m4/2JCvLACr1WWtOpLDAQAAEvHnGwAAqF5NRpzWZscpt+Bv/q5c6K4e4VagcU5efKy7fc+GXjHMa/WMKpoDoAxNyHEys5dL+klJr5S0RdI+SbdJ+jtJf+Xuj868hQAAoPGsJiNO0RwnM/t7ST8j6QpJZ6nTcXqepF+TtFHS5WZ2dhmNBAAADeZT+CnJoBGnt7j7Q4V9T0i6Kfv5XTPbPLOWAQAAzJlox6nYaTKzI8Pj3f1bfTpW1YuWIwhLEMTLEeQXA5524wCk2LLQqxZ+3LqHg3vIcQKayZqR4yRJZvazkn5DnfymlcEwl3TyDNsFAADWkprkOKXMqvsVSc+fy9ElAADQDA3qON0t6clZN2ReuPXfBlCe1x//vd3t//W1GytsCQDkpXSc3i3pX83sOkn7V3a6+y/OrFUAAGBtadCI0/+R9A+SbpXUnm1zAADAmtOwRX6X3P1dM2/JpGzMX/i4jwNQisPsYNVNAFCC2hfADHzBzM43sy1mdvTKz6wbZmZnmdlXzOwuM7tw1tcDAAAVakABzBU/kf377mDfTMsRmNmCpD+S9IOSdku6wcx2uPsds7omAADAMEM7Tu5+UhkNKThD0l3ufo8kmdknJJ0jiY4TAACoTMqIk8zs+ySdqHzl8I/NqE2SdLyk+4PbuyW9dOAjvCbBUQAjadUl8QHAROryUU+pHP6Xkk6RdIuk5Wy3S5plx6lfxnbuV2pm50s6X5I2Lh45w6YAAICZa9CsutMlPc+91CGd3ZJOCG5vlfRAeIC7b5e0XZKOOmRLTfqpAABglZITvCeRMqvuNpW/suYNkk41s5PMbL2kcyXtmMmV3PM/AObKU77Y/QGAqqV8E22WdIeZXa985fCzZ9Uod18ys7dLukLSgqSPuvvts7oeAACoWE3GLlI6Tu+bdSP6cffPSvpsFdcGAADlqn1yuJmZd/zjsGNm07QxJTTHBhwTvnB1eRGBprnigS92t6/ed0iFLQFQmpr8nzsox+kLZvYOM3tmuNPM1pvZq83sYknbZts8AACA+TEoVHeWpLdJ+riZnSTpEUkb1ck5+rykD7r7LbNvIgAAaLyajDhFO07u/pSkP5b0x2a2Tp0k8X3u/khZjZtYGJKLhOdWheOC46xdj5oSQNPsXX6iu/3g0jMHHAmgCcxnmx5jZhsl/ZOkDer0fS5z91/PBoY+IeloSTdJeou7Hxh0rpRyBHL3g+6+p1adJgAAUB9uk/0Mtl/Sq939hZJeJOksM3uZpA+oE0E7VdLDks4bdqKkjhMAAMBM+YQ/g07dsTKUvS77cUmvlnRZtv9iSW8Y1kw6TgAAoAk2m9nO4Of88E4zWzCzWyTtlXSlpLslPeLuS9khu9VZK3eglLXq3i7pEnd/eOSnMK/awXahl5qLsbYFoAL3LG3sbt93YHOFLQFQlinkOD3k7qfH7nT3ZUkvMrNNkj4t6bv6HTbsIikjTs+QdIOZXWpmZ5kZGdMAAGC6Zhiqy12mk699jaSXSdpkZiuDSKvWxe1naMfJ3X9N0qmSPiLppyXdaWa/ZWanpDcTAAAgwnsz68b9GcTMviMbaZKZHSLpNZJ2SfqCpDdmh22TdPmwpiatmunubmYPSnpQ0pKkp0m6zMyudPdfTTnHTLmkdvZbi5UgyFUED0oOLOd/27YcbBOqAypx54HeuuJf3/f04J4ny28MgCbYIuliM1tQZ9DoUnf/jJndIekTZvabkm5WZ5BooJQcp19Upxf2kKQ/k/Tf3f2gmbUk3Smp+o4TAACotxnWcXL3L0k6rc/+eySdMcq5UkacNkv6L+7+9cLF2mb2+lEuBgAA0FfdK4evcPf/OeC+XdNtzhS0+2/nFvZt9w/hSfnwXGu5Jq8i0DC79h3X3f76t58W3EOoDmiqWVYOnybqOAEAACSi4wQAAJAoaVYdAADATNUkVNeQjpPnSw9kwrwmD0sQhKUJ2oVyBO2wVAG1PoEqfOXxY7vbex47srv9xCXH5I475c03l9YmADOUUItpXjSk4wQAAGqNjhMAAEAiOk4ly0Jx0fDccv+K4sXq4OFxrWUBqMB9j23qbj/+6CHdbf/2uiqaAwBdzek4AQCAWjKR4wQAAJCOjlNFoov89g/brVrkN6w2Htz3rZ/5vu720X/2r1NoKICYhx8+rLvtj67vbi9+m9JzQCPVaFYd30IAAACJmjfiBAAA6qcmI050nAAAQPXoOJVsJYepnZDj1O6/LUm2FJQjWFKwXZNXFKipez/5wu52+997eU3rHutlFCw+WWqTAJSoLjlOzek4AQCA+qpJx4nkcAAAgETNGHFy9UJ0raDLGikt4MvxUF0rrBwehu0OsuAvMEsHH97Q3V736EJ3e/1jvWMW95XZIgClcdVmxKkZHScAAFBr5DgBAACkouNUJpc8i8t5L20ruuBvuL2UX+XXloLHR8J2AKZv8eHe11EYnlv3eHDMvvzn8LE3v7y7feQl/zaztgGYvbqMOJEcDgAAkKghI04AAKDWajLi1JyO00r4LSx6uRyE4RaCWXHt3n5r5wfdwll2rYO97YUDve1Hf6oXHjjqY4QHgGlY/2hvOwzPrX+i99lbfDL/zbpwsCbftAAGY1YdAABAGst+6oAcJwAAgESMOAEAgOoRqitRWDk8rAQ+YhXxzn1B/tNyb0AuLEcQ5jsBGN/u935fd3vDI7396x/vfcbWfbv3mVzcly8f0iLHCWiMupQjaEbHCQAA1BsdJwAAgERQ7yaNAAAbfUlEQVR0nCriwVB+WEU8KEHgbeu7X8pXDm8d7N3XXuw9pkWoDpiKsEJ4LDwXbi88tZx7fOtA77ad8d3dbb/+1mk2E2i01qGH5nd8u5p21EXzOk4AAKBenBwnAACAdHSc4szsdyT9Z0kHJN0t6a3u/kh237slnSdpWdIvuvsVSSftVzk83A5n2y1HtpVf9NfC8FwQwgurFT9x7su624d/4tqkpgLoyIXnnoiE5/Yt9bb350N1dqB3X/jZzR8FYCCfjx5LXUacqiqAeaWkF7j790j6qqR3S5KZPU/SuZKeL+ksSX9sZgsVtREAAJTFJ/wpSSUdJ3f/vLuv/Kl4raSt2fY5kj7h7vvd/WuS7pJ0RhVtBAAAKJqHJVfeJunvs+3jJd0f3Lc727eKmZ1vZjvNbOeB9r4ZNxEAAMyS+WQ/ZZlZjpOZXSXpGX3ueq+7X54d815JS5IuWXlYn+P7/jrcfbuk7ZJ01PpjvRujzeUyBZkOFpQgCMoU5I6RpIXguEhpAg+OWVisSVAWmBMPv+3l3e0Nj/Y+P4tPRvKanuptt4JtSdLBoBzBEplNwFjac/D/WMnhtknMrOPk7q8ZdL+ZbZP0eklnuncz03ZLOiE4bKukB2bTQgAAMDdq0nGqJFRnZmdJukDS2e7+ZHDXDknnmtkGMztJ0qmSrq+ijQAAAEVV1XH6Q0kbJF1pnRDate7+c+5+u5ldKukOdUJ4v+DuCePvPrwcQbi9HK74m48O5soRtGKlCYLtg/2iiwBiwhIEi08uB9vDw3Nh+QFJ0sHg9lJve3FrLzVyafc3Jmov0EStQw7p3ZiDUJ1ptnlKZnaCpI+pk0LUlrTd3T9kZkdL+qSkEyXdK+nH3P3hQeeqpOPk7s8acN/7Jb2/xOYAAICqzbb/tiTpl939JjM7QtKNZnalpJ+WdLW7X2RmF0q6UJ2IWNQ8zKoDAABrnLlP9DOIu+9x95uy7ccl7VJn1v45ki7ODrtY0huGtbMZS664pJXFesPQWzj8GG5bEIJr50NtuQWAgxl3YUjOW0GoLphh99TZvZJTG3eQmgVI0hP/7WW52+uDCuFhqC6sCh4Nzx04mD/5UhiqC6L6xdmyAPLC/xO9HT+uLCXOqjOzEyWdJuk6Sce6+x6p07kys2OGPb4ZHScAALDWbTazncHt7Vnpoi4zO1zS30j6JXd/zGz0PGU6TgAAoHJTSA5/yN1Pj57fbJ06naZL3P1T2e5vmtmWbLRpi6S9wy7SnI7T0EV+g6H7ViScJ+UW/bVWcF8w264VPr7VP4QHoCMMzUnS4rf7z56zlPDc0oBZdWF4Lpg5u7BpU2/3I4+kNhtonNaGjb0b8xCeK5rtrDqT9BFJu9z994K7dkjaJumi7N/Lh52rOR0nAABQWzNeNuUVkt4i6VYzuyXb9x51OkyXmtl5ku6T9KZhJ6LjBAAAqjfDjpO7/4v6L+smSWeOci7KEQAAACRq3ohTO1IVPGHxXylfniBXszxcJDi3EHCw4G9QsuDAD78kd971f3/DsJYDjXHgR3qlOcKcJklqPRUszBvLa4pUBFdxId/gs+zhceHqAO05zOUAqhDJa/I5qBwun3mobmqa13ECAAD1Q8cJAABguFmvVTdNzek49Su3PuqCv1I+JBeG8SwIEQSzo1u5Y4Iw3xhFtYCmCMNzYWhOklr7ex+gMDxnKeG5YjmC8PMbCc95uAJAsLBpe9++aPuBpmitX993/1yE52qqOR0nAABQX0PWm5sXdJwAAEDlCNWVynuz5sICC7EZdsuR/ZIULvq73D8MZxY+PpghFBxTLCK+9JpeFfjFq3YKaJr2q17c3Q7Dc2FoThoQnjsYVgiPLNhbCK2HYbhYeC6/gkBNvpmBGZjr8FyJi/xOqiEdJwAAUGdWk8ohFMAEAABIxIgTAACoHqG6ErnkWUVU00KwP1ItfCE4pliOIBTmNQWpGJ47JMxrKlQ1jiDfCU20sK9/FfBcRXBJdiAsLxDZbkdKCxQrH6eUHAnzOoLH2+K63u6lfB4WUGdhCYKkvKZIRfGykRwOAACQwkU5AgAAgFSMOFUlHJaMpb7HhvSLtyNhvNhCwGEFArLusRbYS7+ntx0u2LsUhuoKIexYhfAwPBcu3jugHEEujDfBYr4Whu+L1wRqoK7huTpqXscJAADUDyNOAAAAw7HIbxWyoUlvBbNm2kHALFZRfEJheM4j+4uXD/krT+s95p9vnlKrgNmxM767t/1UZCZdOEPuYGFh3jA8t5xQITwMOxQ/uyPOpJvrysnABKLv7bqE5NxrkxxOKg4AAECi5ow4AQCA2iJUBwAAkIqOU8lWYqNBODea72QDciaCSuADq4r3OT6W71S8rxXEcd1706Db/+lFveP/5Zbh1wYq0ArymnSwl5dkuTIDYUXwQo5TpOxA+HnLVQgfUD5kkhIEg3I/wvIElCbAPAor368yYl7TvOT+MeIEAACQwpWf1DHHSA4HAABItDZHnAZNeYwtBpzrCYdhiP6nKZYjyJUqiFw/F8KLtxAo3cJ3P7e7nVu0N7ZIb6wiuJT/jOXC1mEIr/8iv5Ww4O/LukztRiM1MTyXM4dN6mdtdpwAAMBcIccJAAAgVU0KYDaj4+TeG85vBcPqKTPsBmV5RWfThA+KhO0Kb4DcjLtW//1a7j3GvvcF3e32jbcNaCQwGwsveE53Ox6eC6uAR8Jrxc9RuDBveF+sQngYzhtUOTxSLXxS1go/pcy2Q7ly4bkJ39dzGZ4L1GXEieRwAACARM0YcQIAAPXlIjm8MuFQ/qhhO2nAaryxMEJx/txw1ooU9Fvsv7/14uf3Ln3T7SNfD0ix8F3Pzt22/QmL8YZhu1h4rli0Mhaei+Q3TFTkcoYokolZCd9bTQ/PrTDFZ5zPm+Z1nAAAQP3M599Iq5DjBAAAkIgRJwAAUDlCdWVb+YVbJOcol6MU7G7lxwajpQrCHIbwGmEsetCiwGGcOhe/DkoQhNOugzwoC3I8Wi98Xu607S/eEb8mMMTic0/t3ThYKIOfy1/qvxhvSl7TqvyfWF5TpARBTnF/Qv5GNMdjihXByXfCpHJ5TWOoSy5TFMnhAAAAqZwCmAAAAKnqUgCzeR2n2NB/WJogEraTxqgwHg7Ltwbk2sdCiJGwRb5kQaTtyi++unzrl+PXBzKLz3lW78bBSOVvKb5Qb/geXO4fwsst2FsMISSE5zwlbFc0JwvwhpWefelghS3BvBs1PFf7cFxDVDqrzsx+xczczDZnt83M/sDM7jKzL5nZi6tsHwAAKIn7ZD8lqWzEycxOkPSDku4Ldv+wpFOzn5dK+nD2LwAAaCqXbD4GjYeqMlT3QUm/KunyYN85kj7m7i7pWjPbZGZb3H3PsJOtDO1bNCSXELaTJlsYeFCF41ioLheSi8y8C85ry4WLL/QeH1Z+Xt711XhbsOYsnnpK78bBWEXwwvs3NmOu3T8MFw3PrfqMJS7g20/NQhW5BVpF6A5p4bk1G5KrSXJ4JaE6Mztb0jfc/YuFu46XdH9we3e2r985zjeznWa284A/NaOWAgAA9MxsxMnMrpL0jD53vVfSeyS9tt/D+uzr2wV19+2StkvSUQub69FNBQAA/c34f3Iz+6ik10va6+4vyPYdLemTkk6UdK+kH3P3hwedZ2YjTu7+Gnd/QfFH0j2STpL0RTO7V9JWSTeZ2TPUGWE6ITjNVkkPzKqNAABgPpj7RD8J/kLSWYV9F0q62t1PlXR1dnug0nOc3P1WSces3M46T6e7+0NmtkPS283sE+okhT+akt+UO3+YDzRqvtOq44Ldk+Y7hdfJ5XgEg2xhHlTYjlb/SuPFc4XPNzflPDjv0pfvjDQYTbN4ykm9G2FeU0oVcClaCXzkvKbieVPymiLXW33cnGSTWtrfoK3167vb7QMHZtUazJFBOU1zk8s0L5+jGec4ufs/mdmJhd3nSHpVtn2xpGskXTDoPPNWx+mzkl4n6S5JT0p6a7XNAQAAM+daVVdxDJvNbGdwe3uW1jPIsSsDNO6+x8yOGXJ89R0ndz8x2HZJv1BdawAAQE095O6nz/oilXecpiMofhWEpUYO20lJFcZzCwPn1v4NQ3iFPPfwXGFILrxeLlQXbIfDucXzho8Jr7/Q/7zhtPSlO+8W6m/xWSf3boTv56VIVfCEhXmLt3ML16aG5CLnnTg8l2BuQiADELZrrtyiz/P6XpyX8FzGlJynNG3fXCl7ZGZbJO0d9oBKK4cDAABIqqpy+A5J27LtbcrXluyrISNOAACg1mY84mRmH1cnEXyzme2W9OuSLpJ0qZmdp85KJm8adp7mdZwiIbGksJ0UDzfkZsVFLh2beSflx/Y8EuoYNYS36jGRKuSRGT+LJ5/Y/zySlu7+Wt/HoHq5102SDgbVqHML8EbCc9HtQkht1BlzkS+9gRXBR/2inLPwQlfC561zWLiQd+8zF4btwse391Pcd65FXmvCc2OYTnL44Eu4/3jkrjNHOQ+hOgAAgETNG3ECAAC1U1Fy+Mia3XEaNWw34DEDi2Z2jwlO0yqMOUZn3/U/VW62UyyEt6qNkf3hktOthOMlLX5nr4D70tfvF8qXC8mFQ//hbDkpH5LLheoiM+ZSZssVTRKeG/RlmBLSSAwvjBwemeewRaa1YWNwo/cZbe/bV0Fr1rDEAqdzowbv7b7oOAEAAKSYaGZcqWrWjQYAAKgOI04AAKBartqMODWj4+Tq5UoUK2t3jxme7yQNKFWQku8UqTS+qinRsgWRCuGx8gXFdoVi5QwG5DXFzru49fi+7Vq6b3f88Ui2eOJ39m6EuQmxyt8DKnxPbWFeKakSeFIV8KIq85rmSKw0QbS0QXB8LvdJIv9pCoqL8dbuvVXXvKZQTZ5CMzpOAACg1phVBwAAkIqOU0UGLYi7IhaCUzz0sKpsQfd6CZXGVz0muF5YpyBsb2r5gthziZUzSG1jKBLqWzxuy/BrSFr6xgNp12mYXAhOyg+lh69bGJKLhWVj4bjifeOE5LrHFN7Lo1YCj5438cswIdQwcfikjHBGYhXx/EMSwnah4nljYTxCeKsUQ3IrahGaa0I4rgGa13ECAAD1EuYqzzk6TgAAoGL1qePU7I7TqGE7KRpySlokODbzriih8nju2mEV8kKUJhfGS5mVlztx4pt00Oy7fgrPb3HLM4Y+ZGnPg6Ndo2KLz9zauxH7KylcfLcoFoaLVdyOhOM6d0XeOCkL88bapGpDcvnDaxCeS712QuguGrYb57yDZuIlmPdFhmNht0EIyc0xOk4AAACJatJxonI4AABAIkacAABAtUgOr8BKTDgh1h/Nd5LiQ4UDqo13Dxk0vT8l/2nE3CcprZxBvhmRnKhVbUnIkcqfOLh25OIDLB57zMiPSS6nsGLUXK1BliIVukOD8txS8s5ipQVWXWeyXKbu7nHaO8Uq4PmH1DivaYqS8p06Bw4/WWJphFA0L2rC3+/Ir2/Z16tCQ96z4/Pa/A6a03ECAAD1RY4TAABAszRvxCmlau+gYdsxqo33Lh0pWVB8fGiSEF5RwihnNLRXlBTqy595ZKOGA4vC39c0w3DTMsZCt9Fq34POO62QXEkL8+YfMsW/MOswzD9GVfHeQ+O/q4FhvH7XLkHV4da5Dc/V4X1aBXKcAAAARlCTUB0dJwAAUD06TnNgnGHxcaqNd68xfOZdpymRyuPRNiUO7Y46y2x5ePixyGMxvFw7EsNmI06+i4cJpbFChTOWPBMu/6DI8aPPeBsrJJe7JuG5mUpdwDfpVMN/p0nhvCleb4STTfFU8/c9UOv3aKnqs+QKyeEAAACJmj3iBAAA5p8rPbpSMTpOAACgejUJ1a2djtOk+U6hUXOfpOT8p+7ho+ZBSZP11kfNjxpkOe33MKqk/KpBJs3xmHh69Rj5S8MeqwlzmVKf04h5GjPLNVkL+SITlCwYfNo5+k9pSq9jE5/TmkbHCQAAIIXXpo4TyeEAAACJ1uaI06Ah1ZSh8VFDeFJa5fDw8MSw21ghvX6qTsqbZqgwZlAIcVrG+T0mvG5jLcAbmlFpgfzD53OKeq1NsWRB6Wb0GhKeaygfUsJljqzNjhMAAJgv89QpHoCOEwAAqB7J4eVx9R++HatS7iRD45MuHhwzYCZaakive6pYSKzqN2zVocJJTPi7S3oNp1ntO3/x0Y4XIbm5kPp7m1VIr24LBk8T79k1rxEdJwAAUGPutfkDmo4TAACoXtWRj0SN7jiFw7sTL3A5aVG6cWbida89xptpwtl60xQND4Zq8oHpZ+Lf6ajPfZywRdUhufyJZ3NerFbj3zXhubWniv+fxtHojhMAAKgDr80f0DUoBgIAADAfGHECAADVclHHad7E4uWllywoSnmjjNXGKb4BJ1iYV6pP3HoqJv29T7yQMIvxTttUcyWxylzlMoUa9B6ujZr8ztdMxwkAAMynWD3GeVRZjpOZvcPMvmJmt5vZbwf7321md2X3/VBV7QMAACVx74w4TfIzhJmdlfUt7jKzC8dtaiUjTmb2A5LOkfQ97r7fzI7J9j9P0rmSni/pOElXmdmz3X15Vm0p9nAnGoqfVTXfSXvhE5diqMdfAXNvktdxnhbgzZ94Nuedomk+91HPRWivv7qMLKA5zGxB0h9J+kFJuyXdYGY73P2OUc9V1YjTz0u6yN33S5K77832nyPpE+6+392/JukuSWdU1EYAAFASb/tEP0OcIekud7/H3Q9I+oQ6fY6RVdVxerakV5rZdWb2j2b2kmz/8ZLuD47bne0DAABNNttQ3dT6FzML1ZnZVZKe0eeu92bXfZqkl0l6iaRLzexkSf3Gtft2I83sfEnnZzefuPLgx78ycaP72yzpoRmdu2pNfm4Sz6/u6v/8Bn+X1//5xTX5uUnNf37PKfuCj+vhK67yyzZPeJqNZrYzuL3d3bdn28n9i2Fm1nFy99fE7jOzn5f0KXd3SdebWVudN+JuSScEh26V9EDk/Nslbe933zSZ2U53P33W16lCk5+bxPOrO55ffTX5uUlr4/mVfU13P2vGl0juXwxTVajubyW9WpLM7NmS1qvTe98h6Vwz22BmJ0k6VdL1FbURAAA0ww2STjWzk8xsvToT0XaMc6Kq6jh9VNJHzew2SQckbctGn243s0sl3SFpSdIvzHJGHQAAaD53XzKzt0u6QtKCpI+6++3jnKuSjlOW0f6TkfveL+n95bZooJmHAyvU5Ocm8fzqjudXX01+bhLPr5bc/bOSPjvpecyp0QMAAJCkssrhAAAAdUPHqQ8ze5+ZfcPMbsl+XhfcV/slYczsd8zsy2b2JTP7tJltyvafaGb7guf9J1W3dVzTKq0/D8zsBDP7gpntypYoeme2P/o+rRszu9fMbs2ex85s39FmdqWZ3Zn9+7Sq2zkOM3tO8BrdYmaPmdkv1fn1M7OPmtneLE91ZV/f18s6/iD7LH7JzF5cXcvTRJ5fY743I8+v0f/vTROhuj7M7H2SnnD3/13Y/zxJH1enAulxkq6SNNMlYWbBzF4r6R+yZLkPSJK7X2BmJ0r6jLu/oMr2Tco6pfW/qqC0vqQfH6e0/jwwsy2Strj7TWZ2hKQbJb1B0o+pz/u0jszsXkmnu/tDwb7flvQtd78o6/w+zd0vqKqN05C9N78h6aWS3qqavn5m9v2SnpD0sZXvi9jrlf0H/A5Jr1PneX/I3V9aVdtTRJ5fY743I8/vfWrw/3vTxIjTaBqxJIy7f97dl7Kb16pTz6JJplZafx64+x53vynbflzSLq2NivrnSLo4275Ync5i3Z0p6W53/3rVDZmEu/+TpG8Vdsder3PU+Q/a3f1aSZuyPwbmVr/n16TvzcjrF9OI//emiY5T3NuzIdmPBiGCJi4J8zZJfx/cPsnMbrbOUjivrKpRE2ri6ySpExaQdJqk67Jd/d6ndeSSPm9mN1pnVQBJOtbd90idzqOkYypr3fScq85f7yua8vpJ8deriZ/HJn5vSmvn/72JrNmOk5ldZWa39fk5R9KHJZ0i6UWS9kj63ZWH9TnVXMY6hzy/lWPeq069rEuyXXskPdPdT5P0Lkl/bWZHlt/6idXmdRqFmR0u6W8k/ZK7P6b4+7SOXuHuL5b0w5J+IQslNIp1iu6dLen/Zrua9PoN0qjPY4O/N2v//15ZqiqAWblBS8KEzOxPJX0muzm1ku2zNuz5mdk2Sa+XdGZWfFTuvl/S/mz7RjO7W50FmUsvvz+h2rxOqcxsnTqdpkvc/VOS5O7fDO4P36e14+4PZP/uNbNPqxMK+KaZbXH3PVloZ2+ljZzcD0u6aeV1a9Lrl4m9Xo35PDb5e3PA+7Exr9+0rNkRp0EK8fcflbQy86ARS8KY2VmSLpB0trs/Gez/jix5VdZZdPlUSfdU08qJTK20/jwwM5P0EUm73P33gv2x92mtmNlhWdK7zOwwSa9V57nskLQtO2ybpMuraeHU/LiCMF1TXr9A7PXaIemnrONlkh5dCenVSdO/N5v+/940rdkRpyF+28xepM5w5L2SflaS3L0pS8L8oaQNkq7s/J+sa9395yR9v6TfMLMlScuSfs7dUxMI58Y0S+vPiVdIeoukW83slmzfeyT9eL/3aQ0dK+nT2XtxUdJfu/vnzOwGSZea2XmS7pP0pgrbOBEzO1SdWZ7ha9T3e6YOzOzjkl4labOZ7Zb065IuUv/X67PqzKi7S9KT6swmnGuR5/duNeR7M/L8XtXw//emhnIEAAAAiQjVAQAAJKLjBAAAkIiOEwAAQCI6TgAAAInoOAEAACSi4wQAAJCIjhOAqTKzQ7I1uxamcK7vMLPPTaNdADANdJwATNvbJH1qGkXy3P3fJe0xs1dM3iwAmBwdJwBJzOwl2crpG7NlUm43sxf0OfTNypbbMLNXmdk1ZnaZmX3ZzC7JlpCRmd1rZr9lZv9mZjvN7MVmdoWZ3W1mPxec72+zcwJA5VhyBUASd7/BzHZI+k1Jh0j6K3fPra+WrQ14srvfG+w+TdLz1VkY9P+ps4TMv2T33e/uLzezD0r6i+y+jZJul/Qn2TE7s2sCQOXoOAEYxW+os4jyU5J+sc/9myU9Uth3vbvvlqRsrb0T1es4rSy+fKukw939cUmPm9lTZrbJ3R+RtFfScVN9FgAwJkJ1AEZxtKTDJR2hzshQ0b4++/cH28vK/8G2cl+7cFw7OG5jdl4AqBwdJwCj2C7pf0i6RNIHine6+8OSFsysX6dqXM+WdNvQowCgBHScACQxs5+StOTufy3pIkkvMbNX9zn085L+0xQv/QOS/m6K5wOAsZm7V90GAA1iZqdJepe7v2VK5/snSedko1kAUClGnABMlbvfLOkL0yqAKen36DQBmBeMOAEAACRixAkAACARHScAAIBEdJwAAAAS0XECAABIRMcJAAAg0f8HumBFOxjfFMUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 756x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 312 ms, sys: 33.7 ms, total: 346 ms\n",
      "Wall time: 347 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "xy_dimer1 = NearField()\n",
    "xy_dimer1.coef_filename = 'amn_test1_planewave.dat'\n",
    "#xy_dimer1.coef_filename = 'amn_test4_dimer.dat'\n",
    "\n",
    "xy_dimer1.name = 'dimer_ex1'\n",
    "xy_dimer1.output = xy_dimer1.name + '.out'\n",
    "\n",
    "xy_dimer1.plane = 'xy'\n",
    "xy_dimer1.nx = 150\n",
    "xy_dimer1.ny = 100\n",
    "xy_dimer1.nz = 1\n",
    "xy_dimer1.xlim = np.array([-60.0,150.0])\n",
    "xy_dimer1.ylim = np.array([-60.0,60.0])\n",
    "xy_dimer1.zlim = np.array([0.0, 0.0])\n",
    "\n",
    "\n",
    "xy_dimer1.save_int  = 0\n",
    "xy_dimer1.save_etot = 0\n",
    "\n",
    "xy_dimer1.run()\n",
    "\n",
    "\n",
    "aspect = (xy_dimer1.xlim[1] - xy_dimer1.xlim[0] )/ (xy_dimer1.ylim[1] - xy_dimer1.ylim[0])\n",
    "size = 6.0\n",
    "fig, ax = plt.subplots(1,1, figsize=(aspect*size, size));\n",
    "\n",
    "x = xy_dimer1.grid_x.reshape(xy_dimer1.ny, xy_dimer1.nx)\n",
    "y = xy_dimer1.grid_y.reshape(xy_dimer1.ny, xy_dimer1.nx)\n",
    "z = xy_dimer1.Ecartint.reshape(xy_dimer1.ny, xy_dimer1.nx)\n",
    "\n",
    "pcm = ax.pcolormesh(x, y, z)\n",
    "#ax.axis('square')\n",
    "plt.axis(np.r_[xy_dimer1.xlim, xy_dimer1.ylim])\n",
    "plt.xlabel('x (nm)')\n",
    "plt.ylabel('y (nm)')\n",
    "plt.colorbar(pcm)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dimer scattering\n",
    "\n",
    "This is a simpler case of two particles of the same type scattering a plane wave.  In this case, two gold particles at 532 nm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-26T18:07:10.714304Z",
     "start_time": "2018-03-26T18:07:09.671584Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdcAAAF3CAYAAAD3mAEvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XnwJHd55/n3k1V90pJaLamFrl1JtjAjDg9Y3DOzgDDGQCBiwzjYYBiNYVcxXoyxPR4OExuO2TAb2BBgHDuz3g6EB3bwYCwzSGMz5pBhd70GgQQIIcQhDouWGnRAC8l9VuWzf2RmVVZVVlVWVp5Vn1dExa+OrKxv1S8zn3y+V5q7IyIiIuUJmi6AiIjIplFwFRERKZmCq4iISMkUXEVEREqm4CoiIlIyBVcREZGSNRpczWy/md1gZl83s7vM7FlmdsDMPmlm34r/nt1kGUVERFbVdOb6HuCv3f3xwM8CdwFvBm529yuAm+PHIiIinWFNTSJhZmcCtwOXe6oQZvYN4LnufsTMLgA+4+4/00ghRURECmgyc70ceAD4EzP7kpm918weA5zv7kcA4r8HGyyjiIjIyvoNf/ZTgde7+y1m9h5WqAI2s+uA6wB69H5uL2dWU0oREVnJI/z4QXc/r4p1/8LzHuMP/WhY+P23feXkx939RSUWKVOTwfUwcNjdb4kf30AUXH9oZhekqoXvz3qzux8CDgGcaQf8GXZ1HWUWEZElPuU3/H1V637wR0Nu+fjFhd+/44Jvn1ticeZqrFrY3X8AfN/MkvbUq4GvATcB18bPXQvc2EDxRERECmsycwV4PfBBM9sJfAf4FaKA/2Ezey1wD/CKBssnIiKt4gw9bLoQSzUaXN39y8BVGS+pjldERGY4ENL+S6U2nbmKiIisJKT9mWvTk0iIiIhsHGWuIiLSGY4zbGjyo1UouIqISKeozVVERKREDgwVXEVERMrVhcxVHZpERERKpsxVREQ6w0EdmkRERMrW/lGuCq4iItIhjqtDk4iISKkchu2PrerQJCIiUjZlriIi0hnRxP3tp+AqIiIdYgyxpguxlIKriIh0hgOh2lxFRES2jzJXERHpFFULi4iIlCiauF/BVUREpFShK7iKiIiUpiuZqzo0iYiIlEyZq4iIdIZjDDuQFyq4iohIp6jNVUREpERdaXNVcBURkQ4xht7+auH2l1BERKRjlLmKiEhnRFfFaX9eqOAqIiKdojZXERGRErmrzVVERGQrKXMVEZFOCVUtLCIiUp5onGv7K10VXEVEpEO60eaq4CoiIp3RlaE47S+hiIhIxyhzFRGRThlq4n4REZHy6JJzIiIiFQjVoUlERKQ8XRmK0/4SioiIdIwyVxER6QzH1KFJGmQlV0p4WO76REQK6sI4VwXXrio7eK77eQq+0nZV7jPa/mvjTidmaGp/CUVERDpGmWsX1J2lFjGvjDqjl7o0uZ+oZqdGpqviSAFdCKSrmP4+OsjIurq4j2SVWftCIU43qoUbD65m1gNuBe5195ea2WXAh4ADwBeBV7v7qSbLWLkuHiyKSn9XHVwkj03dPxRwC9M413zeANyVevz7wLvd/Qrgx8BrGymViIi0jmOEXvxWl0aDq5ldDLwEeG/82IDnAzfEi7wfeHkzpauQBZO3baXfQObZ1m1jG7/zhmq6WvgPgTcCZ8SPzwGOuvsgfnwYuCjrjWZ2HXAdwG72VlzMEmiHWUzVxaJ9ZExVxgt1oVq4seBqZi8F7nf328zsucnTGYt61vvd/RBwCOBMO5C5TGvooLEaBdrtoX0jP+0XQHyxdHVoWug5wMvM7MXAbuBMokx2v5n14+z1YuC+BstYnA4a5dABZfNo31ifBVu8PxjDDgzFaWwrd/e3uPvF7n4p8Ergb9z9VcCngV+KF7sWuLGhIoqISMskmWvRW13aeAr5JuC3zOxuojbY6xsuT37b2gmjLvptu0v7Rvn0m7Za0x2aAHD3zwCfie9/B3h6k+UpRBt4fZLfemurxTpE+0U9tqz5RNXCm05njs3S795e2i+as+G/u7tVXi1sZr9pZnea2VfN7D+Z2W4zu8zMbjGzb5nZn5nZzkXr2Oz/goiIbJyhB4Vvy5jZRcCvA1e5+xOBHlG/oJUmOFJwLUJn5e2h2oP20P+iPfR/WFcf2GNmfWAvcIQVJzhqRZtrp2zYBmtBvW0XHlY4JFltsc2pab+oYnutdJts2ga2xTpUelUcd7/XzN4J3AMcBz4B3EbOCY4SCq6r6GBgrTt4LrOoPKUd5LZ6DGDNKtgn6t5m83zeRgTgjdkvbN2r4pxrZremHh+KJyWK1m52NnANcBlwFPhz4Bcz1rNwo1BwFRGRzojGua51Avagu1+14PUXAN919wcAzOwjwLNZcYIjBdc8OpSxti1TXUW67GtnChtYHdYqJe0TXdle55WzcxnthmSvFc8tfA/wTDPbS1QtfDXRZVGTCY4+RI4JjroTNZrS8sBqgU3cNsX091rru7X8f9gpa3aU2bTttdTttC7q7LSQu99C1HHpi8AdRHHyECtOcKTMNUuLN7xO7LwVSb57oWxBnZ3WV3C/2LZtNuv7tjLD7WgWm1zPtdLPcP9d4Hennl5pgiMFVxER6ZSwA5WuCq5pLc1YW3nmX8ZvVfCsea222Y6erTeqwP+6lG123W2sRf/ntWpdqtTBvgnuMKw4cy2DgmtLNRpQ6zrJWPQ5OXf0QoG2gweU2q3ZrlrH55Sy7pr//6V22itbh5pOqq4WLoOCa6IFWWutAbUF33ehrPIt2ekLZQcdOqDUpspMtW3b3bzy1LA9tDrQytoUXBve2WsLqDV+z1W/08oZZ84gu/K6tznAtiygrrINVRKYCpzcrfdxLQq0Ld8Xog5NLTtJy6DgKiIindKFS85td3Dd5Ky1ou9WRZlXnhJxhTZTC0xtsctUkVnmWGdZ21Jt0xfWlM22IottcfZawgxNtdje4NpQYK0koHYokK5bhpmDTY6AqGriDDUH1Ka3pcpmWJr+ziVvM40GWvVHWMv2BddNCaolfo+1y1blbzq1Yy882OQMtFudyZYdVJesb533rqVAb/Pobe2ddnOlbbfcD27Z9q82VxERkdJVecm5smxXcN2ErLWE71C4PEU/O+/nLWtfhYkz6FxZbMYZ99aNja2qt+6c9Wa+N08ZVt0uF/3vCraPlprJVlBl3NhkFC2qItYkEluuLQG1lgnvy/qu89aTPpDMCXJzA+aSg8JaY2MXrLcVqpgLeJWAOm/5qrcXyHeiBisN6xq/pWBgK3G7aawttiVVxKoWbouaM9Ymp34r9NmrfNac9YfHj6/+uYs+Zu/e6I777GcmB5NVAu2SA1vhjKBtgbbK7WZq3bkD6tzlxs+Hx44t//wVBHv2TD4x7/9aILssJbstOdA2PjZWZmxHcBURkY1Qx1VxyrD5wbXFMxPNWUn1n7tGplp2hjrPvEwm2Ls3+yrE4fwstvL22On1ZpSjEmtu26W2qy7LVuMstewMdZ6s7XQim60oky2cxXYpe21B+6s6NG2Rui/mXVlAbSiY5jV9cB4dMNPlngq0S9tj4+WmlXLQTFv3YFRXZ7YyAmrNwTSPQgF3xXbawlXGJVQT197ZqaH2V00i0bSaMta6x4iWPp9rxvraFlAXSZd1ItBOt7XGB4F1Oz7NvG8Vbe6tvqEBdZnM7Wf0Yo7stoXzXKsNth02N7iKiMhGUm/hprQ9a20qW12wni5lq/PMzUJCz8w2ivQqnvu+Fip9usL08snrHcxW55nZfiaaGpZksWW3za5ZTVxbFXET7a+uDk0bq6721XWq80aWHmBtIw6M03IF2rwdn6aWnVykXYF23Qn1F1YDT1UBb+J2kwiPH5/cbuYNB0sUaJst+1KLa3/OOmpsf3XUoakZFWetVWerawfUhRMAVDeusM3C48fH42bTP13ejk+wcqDNXEfJyqrRyB1QYaMy1WUWtseuktVCOR3mCgawxmZ1qlAXMtf2V1yLiIh0zOZkrg1fm3WhMrLWErLVxDZkHdPS33lhFpvRHhs9nX8Iz3ixcjPZModfbXov4LJN1H5ANHNYIm8Wm2Nbyd0W28Yq4pqqhjUUZwNVNbZ0rXayHFPLwXYfGKclv0Wwd+/4IJn8jgsmoxg/PSfQTr1nWqXXNC1yLdVF8/5uaRPCIpknaJA/0K7QtJAryLY1wEIUASuk4Lohqpqvt3BQzRlQLQgYPvro0nJsq/DYMXr79gHgYTg7j3GY8fuX0DZbmiLXUl3UngqjbciCaDltP9kmTtDS+11JgTZXkG3BTElN0PSHG6KKwFqo+neFgJrQgXG55DeaCLKJiWrj5ZNSRE8vyGrHCxUrbJETtjyT6Ke2IQsCMGP4yCNFSrh1wmPHJrPYZYF23ixQCzLZsrPYTejg1IXewi1uqBQREekmZa5lWrcqOG/WmtFJKXq7staiho8+Su+MMyZ+w4VZ7KrXlU0rqfNdoUu+TWepGc8ra13NTPNCIiuLnZfB5riYxNIq4hVrRDo7TaKrzbXTyu68VFpV8IJqYFBQXUc6qOQOtJCrbXbypWIHtJW3oYzOSQqo1ZhuXoCMQJunPXadKuICbbBdrCJWb+EOK7tX50oHxRxtq9MBNXlNB8gS9Xrjg2EYYr1edN998qAJi9tmE3Oy28JWaEvNDKhzgqysZ/joo6MAO3Nytqg9NudQr6qy2K7pQnBVm6uIiEjJlLmuq8j4wnnvzdG2OjobVpVepYZHj9I7cHb0IF2ll2SxqcxjYZUxzLbRRm/KV5BVxzhn1XDMy1bj9w9/9ON8ZZFcRlXEZ5wx2k4sCPJVE2e1xc5phy2zirhL7a8aitNBlV2AfNn7Vm1bVWCtRRJ0+uecgycHqemDYhjOb5ud7sQyek9GsM1jwcnXvO2DeGjN9DrMAgYPPbR6GSS34SOPRAEWwH1+NfGiKmKYGyjLbi/tUvurK7h2x0qBdZ1ewStkqzNtq6nXFVTrM3joIfrnnRs9MI8CbXIASv+PpgJtYhRw5wXbvHKecDGTsY4DanJ/8MCDxcogK0n203SQhalMNiuLXWE87MKsc0PbX7swzlXBVUREOsM1FGd7LZ12LlG0mk9ql2R6/YPnYW5gUaYwkcVO9TBOJP/LmV7GK5rXSxyYzFZT21U6W8WMwf0PrFUGKSgIom0iVRU8sV1M9yRecTzs0uw14z3zdKn9tc0aC65mdgnwAeCxQAgccvf3mNkB4M+AS4HvAb/s7pX1uKj0SiPz3rNmUB0+/PDSYko1Bvc/QP/8g6PgOhto4wUzAm16OE8h806ysgIqRM+ZRWUBBkd+UOxzZW3Dhx+md9ZZ45OuVFXwws5O86qIi9igKuIutLk2mQYNgH/t7v8IeCbwOjO7EngzcLO7XwHcHD9uh1U7oVhQXmC16ECpwNoC/f74FgTQi27W641uBDZ6nl4vWi65xf/LlW7pdfR68XOpz7UguvV6489Nl1E1Ho0bPvxw9nZAtO/P7P/pNvNlc0YTndyXNUa/0is4rS3qLVz0VpfG9jh3P+LuX4zvPwLcBVwEXAO8P17s/cDLmymhiIi0kbsVvtWlFW2uZnYp8BTgFuB8dz8CUQA2s4Nz3nMdcB3AbvZmLbL8c0s8Oyt29ZwFM+jARDva8OjRNUonZRrcex/9iy+KHqSr8IbDyeriUVtsRttsYlkVcdbsSYFNtqWml0uy2uS5XsDgnsMrfDupUlLz1DvrrPH+nVzukIwq4rR5Mzmt2v7a8aphTX+Yk5ntA/4C+A13/4nlnIrN3Q8BhwDOtAPVtr6X3c5aZHyitMrg8L0A9P/bS2CYMaQiFWhxx+hNts2OpP7HYaqdbYqlt6ekLTU9OcRUQE3uD/7++0W/olRpol01GLfFxh2dFg7TgdzzEa9VxA6Ne22jRoOrme0gCqwfdPePxE//0MwuiLPWC4D7K/nsKtsUVg2si8YogrLWNuv3wZIDY5gdaN3jnqKprDaRPnClEtqRrG3HbDagjjLXYJwZq521tYZHj9Lbvz964D6TxWb2JM4zFja2MDB2/SLrXrxPYJ0a2/ssSlGvB+5y93elXroJuDa+fy1wY91lm7DO9IZZHRFmVpLdcSn3+6VRg29/F3b0o1u6o1N/qhNTrze+TTy/5JZnHf0e7NwR3fr90TKDu7/T9M8jiyT794Iaq7k1W9PmHKcKXY1rlXU0JMQK3+rSZOb6HODVwB1m9uX4ud8B3g582MxeC9wDvKKh8omISMs43RiK01hwdfe/hbmnEVdX9bllnIXlbmedeG3O2ShkVgUTmCZU7wjfGVXD2jCAIK5qC8PZ6uJ0G1u6yniRrKrfUZtrvN0kGW583/uqDu6CZP/uHTh7PD46qSKeaoOFjCri6fbXrlf3bpjGOzS1VpExrWl5OzDNBNpoWQXW7hje+U0Aek96/LgHeJgKhh5EgdbHB8yJ9rJ5AXZ6wv1RoA1mA2pv3LlpeMfXy/haUpPhj348vgJTyOT2MN2hKUsN7a/t6tykq+K0TqVZa1pWb89FbSepjMRWDerSGsM7vk7ws1fGD4aj4GdJxuHJgSwVXJdlGdOzLUHc1mp4sk31g9H98PavlfV1pEbJfu9BPMtXuoNTasje3F7EEyvr/nCbZbrQoWmrgquIiHSf2ly7apXsMWc7a+ZzSc9AXWNzY/iueBjMwLBBlD14GFX/W5J1hA5Bkrmmtp/kdNwymhRS24kHAfTGTQreVztr1yX7/ah6OF1LMe8iENPtr0smmChjcok2TOrvruC6cXINvRktnLMDU7oquIVd3mU1/vk7or/P/Fks7txkwyEWeBRkAYaOrdihyXuT1cLeMzwez+r9AD53e7lfRBphFsRVw+kxrSvO5CStsBXBNXdb67qdmBYumzoTnUMXsN4c4c6AINnuhgaDEMKk7czxdHBd1KEplbmmOy15P8DjmZjCHUGjV+CQ8gweeoj+OedEARayOziR8dy8y9RNZaRtyDzLoA5NIiIiJVOHpg2xcg/hrAn5J1Y41TtY1cEbZ7i7B6eSauHoUmAWxu1lQx+3v7qPxzhOCxhXCyftrIAHhvcDwn6SuZoy102S3laS3sOwuP116v2FrgHboWE5anPtkrI6MY0WmTOedaoD08z0Z7IRdnziVk696GkABAMnCAyL5x22wEdVxPmDq406LXlghDsCwn702s6//kJl30MaYEY0D9GSZfKOh50Kmm0Ijutw6r10XFEbH1wrnxdzTiemzMdZzBjc/0C5ZZJWGO6OD2qnHQ8gGMYnVQPHwnGbazKR03RvYZ8JrtH9sGf4joDhjvYfYGR1g/sfoH/euQBYGEy2v6Y7N8UyOzdltb/m1aGew2228cFVREQ2SxdCuoLrEjOZ76pVwjB5pQsNvdkaez56CwDH/vtn0AsCfBAPo+g5wTA+PISkhuXEb0xqjFM9hD1gdD/cYYR9Y+9Hbqnle0gDRscLxu2vDMevT80/PPGenL19Opt5apxrR6wz3aCqhCWHwZ4AD5ygFx3IgqGPA21qqMXo2urJpmM2njGxb4Sp4Drc2f6DixSXHBf6B89jdNYV2OQE/ym5xr1u0rSIFZ8TmNl+4L3AE+NPew3wDeDPgEuB7wG/7O5zJ4Hf2E6GFtha7a3rvh8YX28z6cSU7h2sjkxbY7DbGOwxhrvj286A4a7kZoS7AsJd4+fGj43h7oDh7oDBrmD0/kF8ky2RdeyAyeNL2vQFH4rQHOfvAf7a3R8P/CxwF/Bm4GZ3vwK4OX4819b/giIi0i3uVvi2jJmdCfwz4Pros/yUux8FrgHeHy/2fuDli9aj4LqKJWdzFgTzx7YmRtfiVOaxLQ5c/3ejbDO5jbLYXUGUye4MGO60+BaksluLM9jJ9x+4/u+a/lpSB5tzvMg4fuQ6/kwdw1pRw1dAMrFZkRtwrpndmrpdN7X6y4EHgD8xsy+Z2XvN7DHA+e5+JPp8PwIcXFTG7W5zLVr1saitNXluXiem1HsGP7y/2OdL5wz3jDskBT0IBtHzwZBR+5HFZ9XJybUH43n9w74R7ojv76yp0NK4wQ/vj9tdmZxcIt25CbLHvcL8aRE7zFm7Q9OD7n7Vgtf7wFOB17v7LWb2HpZUAc9biUwppYfw3AWVsW6jwe7o+uYA3gMfxJ2ThqTGvMYLj4Kr4fFFdsI+DBVct9PEMSOjc9PUzE3lfW6+GZtq56R6/VXiMHDY3ZPu+DcQBdcfmtkF7n7EzC4AFmZHG1ktXFs1RVZVTVZ1zXQnpqQjgmyNC//g7xjsIbrtHt+GuxhX/U7d0ssNdkfZ73BPtC7ZIunjRbpD0xy5qoYzqofX0UTVcFXc/QfA983sZ+Knrga+BtwEXBs/dy1w46L1KHMVEZFOqWHi/tcDHzSzncB3gF8hSkY/bGavBe4BXrFoBQquq8o1if/kRdCBiYx2cOQHFRRM2m64O/rrPfBT0f2wDxY3nyVbSHLc8B7jauGdjNpcZbskx4v++QfHY6ItY1rE6UklogXHj4tO6N9GFX8Nd/8ykNUue3XedWxvcC1pHFeuttZ0kFV18NYa7I2OCL2T47ZUG6SCa9L0Goz/eryHhjtguHtDDoxSTBCkAiVM92lKm3vFnI2gifs7J7PdICsIz+shnDanh7A6NG2vy//NZwH49rueRXAq7jl8mtTE/fHfVG/hcQ9h56d+67P1FVbaKdULODm+LO05PHddsxdSh4yr5bRxZqcOnGcqjRIRESnZxmWulfZaW7bu9DVb05StSkq4JxxPyN+3UbXw6BqvQXRa7j0Id8T3d3bgVF2qZdN9OKaG5WRcjm7ifcmY1663u2ri/s2Xazam0f1o2cF9RyoskXTC3gFhL2p09b5h8ZjX5Oo4o0kk+o7HwdV2LWhgk60wuO8I/QseGz0YDsfHl5ZsGrVeZacD5wcKrhWZaW9V9iqxXY85xakgakz1fkCYXEQ9TC4vFwfUvhP0oyxk557TDZRUWivVrmoWzLa7JovluVpOJ7X/eKo2VxERkZJtX+Y6ZwhO7ikP57WnTs0lnLmsCHDmnhM8Go+7OXWqTzhIhk3EPYjjzDXoh+zcGU1CvG/3yQZKKq0z75iSng4xvdyyXsN5pzhs21SIqhbeEGVc11WBVmL7X/It7GNXAHC8F3J6ELW/hlPBdUd/yJ6dUXXwWS/+VgMlldZJn8yPAudUpJk313DOITq1tp0W1fLigYJrdTZork0p39l7jgGwo7eLE4NoNxzGVzzpxTPv7O4POHPXCQAGDZRROiKw1To1df0qOdVP3F8KBVcREemUGuYWXtvGBNc6r8qwbAjOqKdwalamwb33VVwq6ZKDux8FYGcw5Ngg6jl8Ooyqh3cEURqyt3+a/TuPA6DZqAUYHUf6F14wM1vToh7DsKlTIbbXxgTXRmRMeSiSx/m7HwEgMGd3L7pA66k4uO4cBddTnLcrCsIKrrLUaNxrc2ndKMmpOo5vUuZqZmcDFwLHge+5t6XbWMMWzco07zXZegd3/GR0/9E4uJ4Mo91xVxC1sJ7ZP8G5Ox6Nl9pXa/mkQ9KzNSWP07M1daEOdVVdb3M1s7OA1wH/A7ATeADYDZxvZp8D/r27f7ryUpYl75Vw8g7DESnogh0/Ht3fZXsAOOk74sdxD+H+8dRyCq6Skg6gwxV6M83rMdzGyfkXsA6cLyzLXG8APgD8U3c/mn7BzH4OeLWZXe7u11dVQBERkRGn+9XC7v7zC167Dbit9BLVbGlHqCLtqGp7lSUO9qNq4ZCAHfHM/Sfii7futqha+ED/0dFyIgutOhxnibljXTuW4TZplTbXJwOXpt/j7h+poEybYeYKFiJj7/ypJwDw63d/gyDu/XEs3AXA3iCajWl/79hoOZEJZRxTOnuFHOt+m2vCzN4HPBm4k3E/MAcUXFNsuq1WQVWWOKf36Oj+7iBqa90dt7mmXxOZkRrqt5GdlhbpwNfNm7k+092vrLQkIiIieWxQcP2smV3p7l+rtDQdsPQariIr2B+cZOjRNnUi7i2cZK77A03WL/klNWfzJpOQeuUNru8nCrA/AE4SXUzP3f3JlZVMZAvsD5zTHgXRnR4dFJMOTfuDDpyeS2ds1LVdO7Br5A2u7wNeDdxB9XNvdNOiHsLqPSxz7LUee+MZmYK43Wy3haPXRDJNX2Ju+rUGZ2mq3IZN3H+Pu99UaUlERERy2IRJJBJfN7M/Bf4LUbUwUO1QHDN7EfAeoAe8193fXtVnlU69hCWnXdbnMZa0kUWpyGNs/JqIZNig4LqHKKi+MPVcZUNxzKwH/Dvg54HDwBfM7CZ1qJJNs8N67IhPxnbEg/N3xB1TdqhaWPKanl9YGpcruLr7r1RdkClPB+529+8AmNmHgGsABVfZKAEBAVFw7VnynI1eE5FuyjuJxHnA/8TsDE2vqaZYXAR8P/X4MPCMij5LREQ6ZJPaXG8E/l/gU5Q6g+VcWY2WEz+nmV0HXAewm701FElERFphg3oL73X3N1VakkmHgUtSjy8G7ksv4O6HgEMAZ9qBDpzHiMwKCQnj88Zk9EQYn5aHGvUmMqsjV8XJ26jzl2b24kpLMukLwBVmdpmZ7QReCXRnKNC2zfMphZ32IafdoxsW3ZLHrpl2RLoqb+b6BuB3zOwkcJrxDE1nVlEodx+Y2a8BHycaivM+d7+zis8SadJJH/AP8bnYCU/OdeNJJHzQTKGke7bthL4DXzdvb+Ezqi5Ixmd+DPhY3Z8rIiLt1vkOTWZ2qbt/b8HrBlzk7ofLLljnhB7l2PNeE8lwzIccC6MNJ7lYehjPLXzMVC0scyw6pmzD8aYDX3FZ5voOiy61cCNwG/AAsBv4aeB5wNXA7xJ1QNoKHoa6Mo6U5mho/CS+SHpyVZxT8VVxdtjxxsolIutZGFzd/RVmdiXwKuA1wAXAMeAuoirbt7n7icpL2REez7Bjc1NYkUlHw10cDaOhZNOXnOuZegtLeTbmijiwEZkr8ZSDb62hLCIiIguZb0Cbq6zBfft68MnKHhru4+gwylyPxdXDe3WRdMlj6viS1JxthQ2aRGJ7hb7+9VgVZGXKb387Gln2g8F+fjTYB4w7NCVBNiQYLffOn3pCA6WU1kofU7bx+NKBr7xdwdVDsGDqqehi7PPnAAAgAElEQVS/ZGVe0Dz08SzsIhnuH0RDxH94+iweHuwB4GTc5rorbnM97T0CzdIkeRTpIbzgPb4NPY4rlqvbq5ndPD1Dk5kdqqZIIiIi8yXtrkVudcmbuV4GvMnMnubu/zZ+7qqKytS8jAw3et5XuxB60jtPF0+XKUdOnw3Ag6f38ZPBbgBOhtHuuCuIxrkmmaxIaeZVIedtr21Lu24HEuu8wfUo0ZjWPzKz/wL88+qK1CHu4wAaBExczGfVQCxb5f7TUbXwAyf3cWywE4BT8WQSO4No8ojjw1Opd7TkoCbtkD72TAfM9GubaMN6C5u7D4D/2cz+JfC3wNmVlaqAStpOV5W0U2iYqyzxwxPRjKJHT+3h2CDKUE/HwXVHHFxPDPuEo16RD9deRumYrHbSbezs1BJ5g+sfJ3fc/T+Y2R3A66opkoiIyAIdOGfIO3H//zn1+DaiGZu2UjLTybxpEN3DaJam5KzRnf5FFwIwuPe+zPfIdnnw5GMA+MnJ3ZwYRLvhMIy2p14QbV8n+jtGVcXKXAUYHUcIw4msdNEY1zpnZqqtl/GmBNetl2wwq1Q5aziOLPDQsSi4Hj+1g9ODKICGYbS9BEG0vZ3q9zk9jF47q4EySkds4bCZTWpz3Xgeer722nkdlcIQer3ZZZPX1LlJYkf/6goePR5NFHHqVJ9wEGWsHgdXi4Pr6X6PwTB+7a+uYP9LvtVAaaVV0seUorYwGDdBl3cREREp2fZlrknbRNY41kLLxWeBZqmzSmZ7DKvXnsR+cnw3p45HPYR9EOCDOGONM9cwzly9F3AizmpDN/Y3UFZpsfQxJfTJrDbreJP53GwGnNlu2pbxrYkOHE63L7iWaNG1XUedmqIHNZZK2u7kP+zET0bbhp02giS4Jv3f4hYE7/uoyvhkqGYFYaKT5PiplgW+qm3YONfOyN12WofRuFcbtZH0L7yAwX1HGiyUNO5Yn+BU3HnptGHD+PkkgI4yVyMcxPeHLdmmpTH9Cy+A4XD8RI620426hmuaguuGSa6Qk64Knng9NVtTumo43QlKWezWC44HqeAKo2uiJ5tGvK14ABZnteFQ242kFJ2VqcLOTLVO9t+B3UEdmkREREq2vZlr3svPzZvEf2Z9s0N0kraQiQklZGt9553PAqB/zAiiq8phg3Hmmvwdtbn2IJ7Dn3BgfOcd0fsv/zefravI0ibu2W2t0xljVlabub5w6mE3jlGG2ly3wrLZmoC4Ojm5v6FtILJU/x/GVcFBPCe/DRm1uSanZslxw3vRbbyc2l23WvrYsSQQltrW2sYOUwquGyhpd4XZbDXdFT5pdwU8mOw53L/gsQAMjvygpkJLG/RORH+DU+OM1IYQJH1URm2u8cMA4qvQTQRh2S7J8WJmiE0SYFvQ3lqrGnoLm1kPuBW4191famaXAR8CDgBfBF7t7qcWrUNtriIiIpPeANyVevz7wLvd/Qrgx8Brl61gI4Prum0HM+/3cGnVyNJqmNDHZ45hqOrhLXPfG59N7zj0jkP/xOStd9Ln3KaWPR7d7nvjs5v+OlKn5HjhPjqOZI5tnWqTHT1Xg9rba32N2xJmdjHwEuC98WMDng/cEC/yfuDly9az3dXCeTsrFZFUDcPssBzZOv0T43bW4DQE8fhVC1MdmqYmkcCc5KI41h+PhzWdl22nnMeOpSf6eWdlarNqi/uHwBuBM+LH5wBH42uaAxwGLlq2ku0OrkWlr5KTNeZ1Uc/hMBhdLad//kEGP7y/8uJK83rHoXc62laCQSq4DtPjXCe3JQ8giINrOIAwDq5DdWzaGv3zD87vIbxousOsttlVtbEjU2zNNtdzzezW1OND7n4IwMxeCtzv7reZ2XOTj8tYx9ISKLiKiEi3rBdcH3T3q+a89hzgZWb2YmA3cCZRJrvfzPpx9noxsPTC3BvZ5gpRNcc6VR2Z78/R7ppZJZN1VpnVPiIb6UevfTb9Ez669U6k21VDeqeSm8e3+PHJMFo29d7k9qPXqt11ayTHihzDb4pUCcuYu7/F3S9290uBVwJ/4+6vAj4N/FK82LXAjcvWpcw179VvsiwbljMxHWI8R2wQRlXDoIupb4n+iShoBkm18NCxVJtrcpI12+ZqeDLPcN9G41zTE0/IhksFVfdw9eE3y1ZfdufPOuTsmFSyNwEfMrPfA74EXL/sDQquq8rRCSrzajlZl6Nzp3/wPAb3P1B6MaU9+sdDgsE4uNrQx52Yho6NrnQSvyE5XzMb1S2FQ8Pik7FgaKPL08lm6h88L7qTp3Zry9pboZ4Zmtz9M8Bn4vvfAZ6+yvsVXEVEpFs60KKm4LrESpewm+45PD1j02ixVM/hjW31luMvfwYA/ZOOnQ4J4ivb2MCxVPXeot7CyX0LDO8nV8gxLAxG69/z0Vuq/zLSjPSY1ukewim5pjvs6FzCWTS3cAvkDo5Fx7ymh+VMrG6qanj6cnSM21+TKiBVD2+W3onoHx1VCYejgGqDcCK4Mue4aKng6oGNrulqvSDanDQkZyP1D57XXHXwCrocnOuw8cG1DDNXy1nWCSqrc9P0tV+nl5eNcvqFV9E7mQTXcDKgDh1LTrYWBFdSwZUgGLWz+jDKfIP4Wq+nX3gVOz5x65yVSOfk6cSUHvtadg/hlre3AqoWbovMS8llL0i84Oofku45PFpdxhVzpmZumug9LBujdzIkOB3//4dRcE0OjDZMVe8tGpJlNq4WthDvJVXEQTRUrJdsRx040kh+86qCp0wE1RWy1s5nnM30Fl7ZVgRXERHZDEb2lElto+C6gpn22+l22ulpEbMuRzczFjbKXgH655zD4KGHKiq91Ck4OYyyVcCGwyhbTVXp2TCVuS6StLlaavhNEK+rF43pClTzsTH65507P/PMmuZwmVWqeHMu2/nMtyYKrjWY6Nw0PbkETEwwId1mT39S9PfUYBRcCcPJdtZwqlo4kXUCNqoWNkjaXAPHQhsFaA8DPP5c//wdFXwrqU1SJbxksojMKuEcQW9jAmMHvsZWBdcyeg5nZq8wm8FOTeq/sP01dTDtn3MOgDLYjrKT0ez6FoaQZK5hfLCc6MSUHDyXnFAl21UwDq5mIYTBaPYmC0MYRllsB445kiHZ70eBNWvIzSqXk9vETkwpGoojIiJSNgXXLbJonGyq/XWiitgss/ewdFPvSY+HQXTJRxs6DNPVwuE4O5iehH1Rb+FknE4w7jmMBVG7bTL7RBhAXEXce9LjGd7x9RK/ldRhYe/gPMNuaqoSbk21ckuKscjWBdcyhuXkWkfG0JzZFWVN9h/d7R04m+GPfrz4/dIKvSf+DAB2ejgbUCHaliaqgue0uWYZzfaVbocNo+mbbNyOb8F4fb0nPA6A4Z3fXOdrSU16B86eDY7LxrMuGnpTQScmWd3WBddKTQfkORdVn2h/ne7glJrgv3fgbAAF2Zazk1G2OtGbc5gKtNOXGFzUjjaz8lRQTdd4uI8vjZMOtGGADdV7uAuS/XumFmPR1IYF2lk3KmMFcLW5ioiIlK8DwbWRU1wze4eZfd3MvmJm/9nM9qdee4uZ3W1m3zCzX6iqDLnPxBZUm8xdx6KqlqxqntHFkKfOUEPP1ZYizen/9OVxljqM2ltPnY5ug2H0/0xuyTLD6eeX3PKsYxB/dnI7Hd36P3VZ0z+PLJLev6ePAYtm7sp6/xwLj3MdrhI2L36rS1P1R58EnujuTwa+CbwFwMyuJLr6+xOAFwH/3sx6c9dSFw8XVrksDdQzbSLjHWemyic5YKaW6+3fj7TU6dOpwLYkoI4CpuPDYXTzcHxb9FxSzZwE3HQwHQzmB1pppYl9Omu/T16qop11Ra2qEk74GreaNFIt7O6fSD38HPBL8f1rgA+5+0ngu2Z2N9EFaj9bSTlWuZzc6iufP3vTaBnPboOFyV7Emti/dfr/zcXRncFwsnPScDi+n5ieGCAt6387TD2XanN1knUbZgGjI8VUm34ycxPu9C++iMHhe1f7clKPBW3wK80bnBFYNzVj7ZI29Hx4DfBf4/sXAd9PvXY4fm6GmV1nZrea2a2nOVlxEUVEpC26UC1cWeZqZp8CHpvx0lvd/cZ4mbcCA+CDydsyls/8Odz9EHAI4Ew7UM9PtsrMTcnykD17U+Y65k+T2DvrLACGDz9cuPhSjv5FF0YZK8xmq6mMYWLsYvJ6Is/Fracl28bQcUttR0PiTJbMLLZ/0YUADO69b/XPlFKNqoOz+lgkd1fpTZ6hrGrcVlYHQ+3Vu0VVFlzd/QWLXjeza4GXAle7j7agw8AlqcUuBio9IuQe97rOulaY4H9uFXGsd9ZZCrBNS9pTYSKgerq6bXoygDUPmDPrgGgIDkRNCzNVxmRfQEIa0zvrrKXDsTyj3XWkjHbWTakS7kBwbaq38IuANwEvc/djqZduAl5pZrvM7DLgCuDzTZRxrjI3zqxJBUYfM6/NJRxlsVK//sHzoo5CceeidAckQk/19E11bhoOJ8a6ehgWuyUdm9I9S6fXn+4slXSCipfrHzyv6Z9va4322YwOi8Dofzx6Pi2rV/Ccdtay2lpbm7V2SFPjXP93YBfwSYvOqj/n7v/K3e80sw8DXyOqLn6duw8bKqOIiLSMoUkk5nL3n17w2tuAt9VYnOhzV+k5vGBqxLnrmveedBts+ox1Xi/iZD5iqVX/vHOjO8NUlprIW9WX8VoRmVdYSqqpzWbbZj25mk4w+h6DBx5cqwyyogXVvQvbWHNWBW9dptmBr6sZmtaxagenee9ZMEwnestURyegd8YZAAwfeaR4+SWX/jnnRFWsMNs5qUinlKIHwqmTME8H1JgFwWygTa77mgq0/XPO0WUNa5Dsp8BqJ1w5h9xETy/ZnlZsyupCoLYODE9UcE0ps3PTwvXlyWJhIltN74zpTLZ3xhkKsBVJ5n4dtXMm0nMIk6N3ZxntYBbMriejxmNmOxkOpy4OEW1zbq65qyuWFVgzr2pTMFvNZRPbWbe9t3CXlVlFvHB9BbLY6G2TmWxv3z4Aho8+mq/MslRv//7xxPswU623ckBdtyPc9Puzgi1MdFGcPpBbEExl2tHCvf37GR49ul75ZGQUVPN2UkxbYbvRRBHtpuAqIiKdog5N22TRxdLJUUUcvTi+P6+KGLKriYHevn3KXssyTHVSXzYOcYX2seilYkeGiW0nK5OdLsvU5phZZQwaB1ui3r59q3dmW3H7iV4uN7p0pko40YHiKrjOUaj9dZ0q4uT9WZNNwNJq4oSqiIvL1flk1JGpxA4nOc1bjwWW/dnh1HY4p8rYgmC03WCmNvwCRr8fc9pVIX/HtnVPzDaxnXWKMtcNUGhy/6JZbPLe6MXJ56fbYqcy2fHbo/cHe/cSHkvPzyGL9M44g2SCh5ECballH/wmLOiZPl4kI7tdktVmBVqdnOW3MLDmbVsdr2Du55TdK7jTOhBc2zBxv4iIyEZR5ppD4SriBdlrst7cVcSJeb2JE6ksVtnrcknW4VmXiktMTBix4uW95rynkHnrSW0n02WZqTJelsXGQ3fUvJBPb9++SrPV6OVqakO6WiWMq1pYymiDTcyrJoaFgTbYuzdaXEF2RrB372wHpTKulznnPZVZsJ3MVBkva5uNJ/tPTjZ0gpYt2a8WTrQPhcc4rxT4NnCSiKU68BU2J7jmCGTrf0TBi6vnbIOFBdnxmoE2ORiAAu3ot1jU43eVttSaD24Lt8F5vYhZsJ1NXM0ndUUdAHedoMXS+9DIopmCqg6s25StpnRlbmG1uYqIiJRsczLXRI62zvVWv0b2CrnaYaPFlmQn89aTJ4vds4fw+PGF5dg0E9nqvIy1omy16jGJuTPZZW2z6eUnqorjP3v2jNr0ty2LnchaK8hWx4tUl7XWoq5yaW7hzbTWHMQ5g//SquJF1cSJeYGW+ECZLLbhgTbYu7ecgJqx/OzL9e/0uYNt0bbZOYF2q7abIm2qUG5Azbm+Uj6n2AdUu/6ULlQLb2ZwrTh7HX9MtVls7s/J2qgXtcvCRLDdxAPlTKa6To/fFgbURdZtw5/Z3tLba2rdm3yCFuzZU2lAjRZbcbspGLzatn2uzVGHJhERkbJZS2vF0zY3uNbQezj6mDWriKM35/6c3J+1bN1TVcabkIVktqsW7QXcsWx1nqJZbOb7ppfbsCw2/R1GSm4DrStbLfx5xT+ons/pkM0NrjUrXEUcvTn6u0I1cWLpZ67aNku3DpR1BtSF713FugeigieMKwXaZUN4sqZWjF/r1PYzHVAr6lCkoFr25zbzsavY/OBaU/tr9FFrZLHRCsb3CwTapZ+dp20WJjOSXbvjO3Ev0YYPljOdk7oQUMs+AC0Yy5p/FcU6zOXv+DQVaFvQy3itzBSqD6gFPqO0z+0YdWhqi5qqiMcft0YWO15JNQfNrM9JW5LdTgfb0cslB925Qx9WmPhh/HS+5Ra+J68mzuTXDLa5A23ejk8w94o8MwGuoqC7diCFQv/LJreb2oNqk1mrhuKIiIiUS5nrFiste02smcVGqyg/k01kZgrrWPOi0kWqgTuVrS6yxnaT63KIqfXmmlJxqvPTSDyPMVSw/UBl7aeTb13jKF/idlNr1tq27b2ltiu41tj+Gn3cilW0i1c2vl9Ch5bxqlboELXosxft3Kt+94LjCGttV+3KAabgdrO0/0Dejk+LyjDvN8+7vdR8MtSKzmwTq+pA+laVDnz17QquUHv76/hj1+zsNLmyycdrfJeVs9sclz2bUWGHoHVmVGqqw0lj1uwwt2p7bOb78my7ZQWNNvW4LXm7aSywtmD778rE/dsXXEVEpLum5wdvqe0NrjVXEY8/toS22NmVju+v+Z0KtdNOl6Fijcz/24Iz9lIVzGJLmbd43vsaUlkWWMF32+aMtWu2N7hCowE2UWmgjT5gzdVl78yll7tAGaYWWn8dK65zY6ywH5Qxb3Fa67aj4iuvcNUNZ2kt3A9ULdwFDbXBjj++wkAbfcDscyV81zw7/Crfp8rgp8CaQ9mZ7PR656xz0f9m1f1hE3vMNhpY27wPKLiKiIiUS5lrlzRURTxZhIqz2PEHzT5XwXev5Kxb2Wr1yq4mTtY5XnDldTeuTX0K6ilE0yWYzymvR3mFFFxbqrZAO/7A7OcbPuFYZSdvy4D+jVCguST3cLOaTu5W1uA20IqAmtC+UAoF17SG21/nqT3QTn748mXW/b2amldVB5Hlqhgbu+xzppWxP7bwf92qgNo1HfjpFFxFRKRT1ObaVSWOGy1bqTM9laWhrGBj5gLuigL9EkrZXjfo/9XqbLVLv3MHJpFoV+Roo5ZucB76zG0brP19W/r/7AwPC/2G27adJlq/jyb/z47tF+bFb0vXbXaJmX3azO4yszvN7A3x8wfM7JNm9q3479mL1qPgmkdHNr5NDLilfZ+O/A87Y80ryWzadpro1PfS/jDPAPjX7v6PgGcCrzOzK4E3Aze7+xXAzfHjuVQtvIoWDNdZVRtmWMqryaE7UkBF0262cducp/UBNEvX9wmn0g5N7n4EOBLff8TM7gIuAq4Bnhsv9n7gM8Cb5q1HwVVERDojuirOWtH1XDO7NfX4kLsfyvwss0uBpwC3AOfHgRd3P2JmBxd9iILrqjqYvWYpe/rCMj+3pA+q53MkUuJ+UeaUiOvqZGa6yKbsF+t9jQfd/aplC5nZPuAvgN9w95+YrbbtKbgW0eLexGXq7IFlUw4gXVPDOPHObpNtoP0iNzPbQRRYP+juH4mf/qGZXRBnrRcA9y9ax+ZGhrqoo0w7dLTX40bS/6JdNvB/Ye6Fb0vXHaWo1wN3ufu7Ui/dBFwb378WuHHRepS5iohId1TcoQl4DvBq4A4z+3L83O8Abwc+bGavBe4BXrFoJQquZdmQttjO2bAz8o2j/aIZG71feKWTSLj73xL1m8pydd71KLiWaUvaYhu30QeODaQAW68t2D80/eE2U6At3xYcNDaW9odqad9onUa3cjP7bTNzMzs3fmxm9kdmdreZfcXMntpk+UREpIXci99q0ljmamaXAD9P1DCc+EXgivj2DOD/iP92m6rFitMZ+eZRFluObd03HKwDX73JLfvdwBuZ7Pd1DfABj3wO2B+PJ+o+DU/IT7/V9tD/upht/706kLk2ElzN7GXAve5++9RLFwHfTz0+HD+XtY7rzOxWM7v1NCcrKmlF0geUbd9JEvo9RNvAYvp9OqWyamEz+xTw2IyX3ko0ZuiFWW/LeC7zVCOeC/IQwJl2oAN9x0REpBQdOOJXFlzd/QVZz5vZk4DLgNvjuRovBr5oZk8nylQvSS1+MXBfVWVsjW1tg9IZuMyzrftEQvvGQmtO3F+L2js0ufsdwOhqAmb2PeAqd3/QzG4Cfs3MPkTUkenh5CoEWyNrp9qUg4sOGFLENgVa7SP5KLiu7GPAi4G7gWPArzRbnJboWsDVAUKqMr1ttXk/WEb7STHOulfFqUXjwdXdL03dd+B1zZVGRERkfY0HVylo2VlvVWf0OtuWNpm3PbYto9V+Uxoj39Vtmqbguqm0M8s2q/vkU/tbvRRcRURaSMGw2xRcRUREStSRDk0ta5gQERHpPmWuIiLSKerQJCIiUjYFVxERkTLVe3WbotTmKiIiUjJlriIi0h1OJzJXBVcREemWDgzFUXAVEZFOUW9hERGRsnUguKpDk4iISMmUuYqISHc4ELY/c1VwFRGRDunGOFcFVxER6RYFVxERkZJ1ILiqQ5OIiEjJlLmKiEh3qEOTiIhI2Ry8/VM0KbiKiEi3qM1VRERk+yhzFRGR7lCbq4iISAU6UC2s4CoiIt2i4CoiIlKmbkx/qA5NIiIiJVPmKiIi3eFAqHGuIiIi5epAtbCCq4iIdIuCq4iISJm8E+Nc1aFJRESkZMpcRUSkOxxcE/eLiIiUrAPVwgquIiLSLR3o0KQ2VxERkZIpcxURke5w1yQSIiIipetAtbCCq4iIdIorcxURESmTroojIiKylZS5iohIdzga5yoiIlI6zdAkIiJSHge8A5lrY22uZvZ6M/uGmd1pZn+Qev4tZnZ3/NovNFU+ERFpIfcocy16y8HMXhTHoLvN7M1FitlI5mpmzwOuAZ7s7ifN7GD8/JXAK4EnABcCnzKzx7n7sIlyiojIdjGzHvDvgJ8HDgNfMLOb3P1rq6ynqcz1V4G3u/tJAHe/P37+GuBD7n7S3b8L3A08vaEyiohIC3nohW85PB24292/4+6ngA8RxaaVNBVcHwf8UzO7xcz+bzN7Wvz8RcD3U8sdjp8TERGJVFstXEocqqxa2Mw+BTw246W3xp97NvBM4GnAh83scsAyls881TCz64Dr4ocnP+U3fHXtQtfrXODBpguxgq6VF1TmOnStvKAy1+FnqlrxI/z445/yG85dYxW7zezW1OND7n4o9Th3HFqksuDq7i+Y95qZ/SrwEXd34PNmFhJtPIeBS1KLXgzcN2f9h4BD8fpudferyip7HbpW5q6VF1TmOnStvKAy12EqeJXK3V9U1bpjuePQIk1VC38UeD6AmT0O2El0VnYT8Eoz22VmlwFXAJ9vqIwiIrJ9vgBcYWaXmdlOok62N626kqbGub4PeJ+ZfRU4BVwbZ7F3mtmHga8BA+B16iksIiJ1cfeBmf0a8HGgB7zP3e9cdT2NBNe4B9Y/n/Pa24C3rbjKQ8sXaZ2ulblr5QWVuQ5dKy+ozHXoWnknuPvHgI+tsw7zDlxdQEREpEt0VRwREZGSdTq4mtk/NrPPmdmXzexWM3t6/LyZ2R/FU1d9xcye2nRZE12d9tHMftvM3MzOjR+3+Td+h5l9PS7Xfzaz/anXWvk7lzHdWtXM7BIz+7SZ3RVvv2+Inz9gZp80s2/Ff89uuqxpZtYzsy+Z2V/Gjy+Lx9h/y8z+LO600hpmtt/Mboi34bvM7Fkd+I1/M94mvmpm/8nMdrf9d66cu3f2BnwC+MX4/ouBz6Tu/1ei8UrPBG5puqxxuZ4HfArYFT8+GP+9Ergd2AVcBnwb6DVd3lS5LyFq3P974Nw2/8Zx2V4I9OP7vw/8fpt/Z6JOE98GLifqOX87cGXT5coo5wXAU+P7ZwDfjH/TPwDeHD//5uT3bssN+C3gT4G/jB9/GHhlfP+PgV9tuoxT5X0/8D/G93cC+9v8GxNNsPBdYE/q9/2Xbf+dq751OnMlGth7Znz/LMZjka4BPuCRzwH7zeyCJgo4pavTPr4beCOTA6nb+hvj7p9w90H88HNE49Sgvb9zKdOtVc3dj7j7F+P7jwB3ER1YryEKCMR/X95MCWeZ2cXAS4D3xo+NaBjgDfEibSvvmcA/A66HqPOnux+lxb9xrA/sMbM+sBc4Qot/5zp0Pbj+BvAOM/s+8E7gLfHzbZ1GsXPTPprZy4B73f32qZdaW+YpryHKsKG9ZW5rueYys0uBpwC3AOe7+xGIAjBwsLmSzfhDohPDZN67c4CjqZOvtv3WlwMPAH8SV2W/18weQ4t/Y3e/l+j4ew9RUH0YuI12/86Va/31XJdMo3g18Jvu/hdm9stEZ3svoKTpq4qoetrHKiwp8+8QVbPOvC3juVaU2d1vjJd5K9F46Q8mb8tYvg3d5dtarkxmtg/4C+A33P0nUTLYPmb2UuB+d7/NzJ6bPJ2xaJt+6z7wVOD17n6Lmb2HqBq4teL232uImlqOAn8O/GLGom36nSvX+uDqi6dR/ADwhvjhnxNX/VDS9FVFLCnv2tM+VmFemc3sSUQ7zO3xAfRi4Itxx7FWljlhZtcCLwWujn9vaLjMC7S1XDPMbAdRYP2gu38kfvqHZnaBux+Jmwbun7+GWj0HeJmZvRjYTdSE9IdETRj9OKtq2299GDjs7rfEj28gCq5t/Y0hSmi+6+4PAJjZR4Bn0+7fuXJdrxa+D/jv4vvPB74V378J+Bdxj9ZnAg8nVSoN69S0j+5+h7sfdPdL3f1Soh3/qe7+A9r7G2NmL1f6/u4AAAL/SURBVALeBLzM3Y+lXmrl70xJ061VLW6vvB64y93flXrpJuDa+P61wI11ly2Lu7/F3S+Ot91XAn/j7q8CPg38UrxYa8oLEO9b3zezZOL7q4lmrGvlbxy7B3imme2Nt5GkzK39nWvRdI+qdW7APyGq27+dqO3n5+Lnjehit98G7gCuarqscbl2Av8R+CrwReD5qdfeGpf3G8Q9oNt2A77HuLdwK3/juGx3E7Vhfjm+/XHbf2ei3tffjMv21qbLM6eM/4Soau8rqd/2xUTtmDcTndzeDBxouqwZZX8u497ClxOdVN1NVOO1q+nyTZX1HwO3xr/zR4maklr9GwP/Fvh6fGz7v4h65Lf6d676phmaREREStb1amEREZHWUXAVEREpmYKriIhIyRRcRURESqbgKiIiUjIFVxERkZIpuIrUzMz2xHNL90pY13lm9tdllEtEyqPgKlK/1xBNgzlcd0UeTTl3xMyes36xRKQsCq4iJTGzp1l0gfbdZvaY+OLRT8xY9FXEU8GZ2XPN7DOpi2N/MJ5CDjP7npn9b2b2WTO71cyeamYfN7Nvm9m/Sq3vo/E6RaQlWj9xv0hXuPsXzOwm4PeAPcB/dPevppeJ5w6+3N2/l3r6KcATiObK/v+IJpz/2/i177v7s8zs3cB/iF/bDdxJdAFqiKbK+70qvpOIFKPgKlKu/5VoIv4TwK9nvH4u0WW50j7v7ocBzOzLwKWMg2sygf8dwD6PLlL+iJmdMLP9Hl1I+37gwlK/hYisRdXCIuU6AOwDziDKMKcdz3j+ZOr+kMmT3uS1cGq5MLXc7ni9ItISCq4i5ToE/C9EF2j//ekX3f3HQM/MsgJvUY8juhqJiLSEgqtISczsXwADd/9T4O3A08zs+RmLfoLo8m1leR7wVyWuT0TWpEvOidTMzJ4C/Ja7v7qk9f0/wDVxViwiLaDMVaRm7v4l4NNlTSIBvEuBVaRdlLmKiIiUTJmriIhIyRRcRURESqbgKiIiUjIFVxERkZIpuIqIiJTs/wc7wcUDk+20yQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 898 ms, sys: 89.8 ms, total: 987 ms\n",
      "Wall time: 994 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "xy_dimer2 = NearField()\n",
    "xy_dimer2.coef_filename = 'amn_test5_pw.dat'\n",
    "\n",
    "xy_dimer2.name = 'dimer_ex2'\n",
    "xy_dimer2.output = xy_dimer2.name + '.out'\n",
    "\n",
    "xy_dimer2.plane = 'xy'\n",
    "xy_dimer2.nx = 250\n",
    "xy_dimer2.ny = 250\n",
    "xy_dimer2.nz = 1\n",
    "xy_dimer2.xlim = np.array([-80.0,80.0])\n",
    "xy_dimer2.ylim = np.array([-60.0,60.0])\n",
    "xy_dimer2.zlim = np.array([0.0, 0.0])\n",
    "\n",
    "\n",
    "xy_dimer2.save_int  = 0\n",
    "xy_dimer2.save_etot = 0\n",
    "\n",
    "xy_dimer2.run()\n",
    "\n",
    "\n",
    "aspect = (xy_dimer2.xlim[1] - xy_dimer2.xlim[0] )/ (xy_dimer2.ylim[1] - xy_dimer2.ylim[0])\n",
    "size = 6.0\n",
    "fig, ax = plt.subplots(1,1, figsize=(aspect*size, size));\n",
    "\n",
    "x = xy_dimer2.grid_x.reshape(xy_dimer2.ny, xy_dimer2.nx)\n",
    "y = xy_dimer2.grid_y.reshape(xy_dimer2.ny, xy_dimer2.nx)\n",
    "z = xy_dimer2.Ecartint.reshape(xy_dimer2.ny, xy_dimer2.nx)\n",
    "\n",
    "pcm = ax.pcolormesh(x, y, z)\n",
    "#ax.axis('square')\n",
    "#plt.axis(np.r_[xy_dimer2.xlim, xy_dimer2.ylim])\n",
    "plt.xlabel('x (nm)')\n",
    "plt.ylabel('z (nm)')\n",
    "plt.colorbar(pcm)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scattered field from dipole emission\n",
    "\n",
    "Even though the dipole emitter is not included, the field is still so intensely focused in the hot spot that the near field is entirely dominated by it.\n",
    "\n",
    "Also, near-field calculations from dipole emission require much larger VSH expansions and therefore take much longer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-26T18:07:41.141637Z",
     "start_time": "2018-03-26T18:07:20.079871Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdEAAAF+CAYAAADdiSCgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X2QZmV55/Hvr7uHeYGZAEGUMKwDCboxiVFE1JgXFTVILMgfcYusMSS6SyUxruZlfaO2UrsbtzTuqkkl2WRKibrRGIMaiatRMRJ1S1BAERANqCgjGCRCHBgYZ7qv/eOchp6he6b7TD99zmG+n6pT/Zxz7ufpa05N99XXfd/nPqkqJEnSyk31HYAkSWNlEpUkqSOTqCRJHZlEJUnqyCQqSVJHJlFJkjoyiUqSBiPJRUluT3LdMtr+dJKrk+xN8guLnN+S5JtJ/ngy0ZpEJUnD8lbgrGW2/QbwK8A7lzj/34F/PPSQlmYSlSQNRlV9AvjOwmNJfjDJ3ye5Ksknk/zbtu3NVfUFYG7/z0nyBODhwEcmGa9JVJI0dNuBl1TVE4DfBf70QI2TTAH/C/jPkw5sZtLfQJKkrpIcBfwE8DdJ5g+vP8jbfgP4YFXdsuA9E2ESlSQN2RRwV1U9bgXveQrwU0l+AzgKOCLJ3VX1ykkEJ0nSIFXVd4GvJXkeQBo/fpD3PL+q/k1VbaPp/n37JBIo9JxEkxyd5OIkX0pyQ5KnJDk2yUeT3Nh+PabPGCVJayfJXwGfBh6dZEeSFwHPB16U5BrgeuDctu0Tk+wAngf8eZLr1zzePh+FluRtwCer6s1JjgA2Aa8GvlNVr03ySuCYqnpFb0FKkrSE3pJoki3ANcAptSCIJF8GnlZVtyU5Abisqh7dS5CSJB1An925pwDfBv4iyeeSvDnJkcDDq+o2gPbr8T3GKEnSkvqcnTsDnEZz788VSf4QWPbAb5ILgAsAppl+wia2TCZKSdKK7OTOO6rqYZP47J99+pH1L9+Z7fz+q76w+8NVtdwVkQ6qzyS6A9hRVVe0+xfTJNF/TnLCgu7c2xd7c1Vtp7kBly05tp6UM9ciZknSQVxaF399Up99x3dmueLDWzu/f90JXzluFcPprzu3qr4F3JJkfrzzTOCLwCXA+e2x84H39xCeJEkH1fdiCy8B3tHOzP0q8Ks0if3d7bTmb9BMXZYkCShm60FL5fam1yRaVZ8HTl/klH2zkqQHKWCO/m7N3F/flagkSSsy9+CHtvTGZf8kSerISlSSNBpFMdvjSnv7M4lKkkbFMVFJkjooYNYkKklSN0OqRJ1YJElSR1aikqTRKHBikSRJXQ3nLlGTqCRpRIpyYpEkSZ0UzA4nhzqxSJKkrqxEJUmj0SxAPxwmUUnSiIRZ0ncQ9zOJSpJGo4A5x0QlSRo/K1FJ0qjYnStJUgfNAvQmUUmSOpkrk6gkSStmJSpJ0oAluRnYCcwCe6vq9KXamkQlSaNRhNm1ubHk6VV1x8EamUQlSaPimKgkSR2swpjocUmuXLC/vaq2L/JtPpKkgD9f5Pz9TKKSpBEJs3VI3bl3HGiMs/XUqro1yfHAR5N8qao+sVhDVyySJGmBqrq1/Xo78D7gjKXamkQlSaPRPMVlqvN2MEmOTLJ5/jXwbOC6pdrbnStJGpUJ3yf6cOB9SaDJke+sqr9fqrFJVJI0GlWHPCZ6kM+vrwI/vtz2dudKktSRlagkaVTmXPZPkqSVa+4THU4nqklUkjQikx0TXSmTqCRpNOZvcRmK4UQiSdLIWIlKkkZl1gXoJUlauTV8FNqymEQlSaMy58QiSZJWbmi3uAwnEkmSRsZKVJI0GkWcWCRJUldDuk/UJCpJGo0qBrVi0XAikSRpZKxEJUkjEp/iIklSF4XduftIMp3kc0k+0O6fnOSKJDcm+eskR/QdoyRpOGaZ6ryttt6TKPBS4IYF+68D3lhVpwJ3Ai/qJSpJ0uAUYa66b6ut1ySaZCvwc8Cb2/0AzwAubpu8Dfj5fqKTJOnA+h4TfRPwcmBzu//9wF1Vtbfd3wGcuNgbk1wAXACwgU0TDlOSNBQu+wckeS5we1VdtfDwIk1rsfdX1faqOr2qTl/H+onEKEkalqJZgL7rttr6rESfCpyT5GxgA7CFpjI9OslMW41uBW7tMUZJ0qCE2QHd4tJbJVpVr6qqrVW1DTgP+Ieqej7wceAX2mbnA+/vKURJ0sAMrRIdTsfyA14B/HaSm2jGSN/SczySJC2q74lFAFTVZcBl7euvAmf0GY8kabiG1J07iCQqSdJyVGUi3bJdmUQlSaPisn+SJD0EWIlKkkajwKe4SJLUTQbVnWsSlSSNRnOfqJWoJEmduHauJEkPAVaikqTRmH+e6FCYRCVJozI3oE5Uk6gkaTSqYNZKVJKkbobUnTucmliSpJGxEpUkjUYzsWg49Z9JVJI0Kj4KTZKkDoa2YtFwamJJkkbGSlSSNCKOiUqS1JmPQpMkqQMXW5Ak6RAMqTt3OJFIkjQyVqKSpNHwKS6SJB0CJxZJktTBWiy2kGQauBL4ZlU990BtHROVJGlfLwVuWE5Dk6gkaVTmaqrzdjBJtgI/B7x5ObHYnStJGo+a+MSiNwEvBzYvp7GVqCRpNIpmYlHXDTguyZULtgvmPzvJc4Hbq+qq5cZjJSpJGpVDrETvqKrTlzj3VOCcJGcDG4AtSf6yqn5pqQ+zEpUkCaiqV1XV1qraBpwH/MOBEihYiUqSRmRozxM1iUqSRmUtkmhVXQZcdrB2JlFJ0mi47J8kSYdgSMv+ObFIkqSOrEQlSeNRTiySJKkTZ+dKknQIhpREHROVJKkjK1FJ0mh4i4skSYegTKKSJHUzpPtETaKSpNGogd3i4sQiSZI66i2JJjkpyceT3JDk+iQvbY8fm+SjSW5svx7TV4ySpOGpSudttfVZie4Ffqeqfhh4MvDiJI8BXgl8rKpOBT7W7kuSBO3s3K7baustiVbVbVV1dft6J3ADcCJwLvC2ttnbgJ/vJ0JJ0hANqRIdxMSiJNuAxwNXAA+vqtugSbRJjl/iPRcAFwBsYNPaBCpJ6tXQlv3rfWJRkqOA9wAvq6rvLvd9VbW9qk6vqtPXsX5yAUqStIRek2iSdTQJ9B1V9d728D8nOaE9fwJwe1/xSZIGpprbXLpuq63P2bkB3gLcUFVvWHDqEuD89vX5wPvXOjZJ0nDNkc7bautzTPSpwAuAa5N8vj32auC1wLuTvAj4BvC8nuKTJA1M4bJ/AFTVp2DJPwvOXMtYJEnqYhCzcyVJWh6f4iJJUmeTmCDUlUlUkjQqjolKktRBc6vKcJJo74stSJI0VlaikqRRcWKRJEkdObFIkqSOhjQmahKVJI1GMZlHmnXlxCJJkjqyEpUkjcqAhkRNopKkERnYfaImUUnSuAyoFHVMVJKkjqxEJUmjYneuJEkdudiCJEkdFFaikiR1U8CAkqgTiyRJ6shKVJI0Ko6JSoej7NfxU3Mrf99y3yM9lJlEJUnqYlgL0JtEpTWSqQd+8Guulq5M9z++4L01O7HwpPEYUCXqxCJpjdTsLDU7C5naJ6HeL1NLJ9BMNe+VNChWotIay1SouSkyNddUpAdo175o3rNG8UmD5gL0kiQdggH9RWkSldba9DRh9oDV6HwX7v2vp6fXOkppwIZTiTomKq2xTE/vM/aZqewzRvqg8dJMNe+RNDhWopKkcRlQd66VqLTGZu++G6an9umyfZB2MlHTlTvVvEdSow5hW2VWolIPsm4dVdWMjbZ3rjyoS7ftws26dX2EKA3ThBegT7IB+ASwniZHXlxVv7dUe5OoJGlUJrx27m7gGVV1d5J1wKeSfKiqLl+ssUlU6kHWHwFzzczctHOG5mfpznfzZmbmgbaS1kRVFTA/frKu3ZZM2yZRqQ+bNpH2z+nauxeApP05nQqZmSEb1t/fVtICh1aJHpfkygX726tq+8IGSaaBq4AfAv6kqq5Y6sOWnUSTHAP8AHAvcHOVj5OQJPXg0MZE76iq0w/48VWzwOOSHA28L8mPVtV1i7U9YBJN8n3Ai4FfBI4Avg1sAB6e5HLgT6vq4x3+EdJhbe/NX2dm2yNJAt/b0xycXxt3epocse7+CnTvzV/vKUppmLJGt7hU1V1JLgPOAlaeRIGLgbcDP1VVdy08keQJwAuSnFJVb1mFeKXDytz3b2Zqeorce19zoO3WZWYGNm5g7ugjm/2bewlPGqYJ3aoyL8nDgD1tAt0IPBN43VLtD5hEq+pZBzh3FU2fsSRJDxUnAG9rx0WngHdX1QeWarySMdHHAtsWvqeq3ts9Tunwdu8Jm9gwHaZ2NveBZk9Tida6GeY2r+e+4zcCzc1qkuZloveJVtUXgMcvt/2ykmiSi4DHAtcD8xOKCjCJSh3tPGmGmtrI+vXNj+HU7iaJzq2fYfcx67j7xOa4SVTaz4CW/VtuJfrkqnrMRCORDjM7txWZm2FuXbP038yu5sdx76Zp7j1uip3bmt8Ux/UWoTRQA0qiy10799NJTKKSJC2w3Er0bTSJ9Fs0SyKFZmGHx04sMukhbuYH7+ae+zZTU83fsuvuacZ59hwZdj2iOS9pEQOqRJebRC8CXgBcywNjopIOwdMeeSMf+s5jyWzbjbuhTaKb4b4T9vKcR94IwFd6i1AaoAkvQL9Sy02i36iqSyYaiSRJy7BWiy0sx3KT6JeSvBP4O5ruXGCyt7gkOQv4Q2AaeHNVvXZS30vqwznHfo5PPuwHufeeLQDUdFuJbimOetg9nHPs5wB4Iz/cW4zSII0wiW6kSZ7PXnBsYre4tDe5/gnwLGAH8Nkkl1TVFyfx/aQ+PGn9XRy/eSdfPbJZmSh7m8e5zB45y/Gbd/Kk9Xcd6O2SBmBZSbSqfnXSgeznDOCmqvoqQJJ3AecCJlFJ0mAsd7GFhwH/kQevWPTCyYTFicAtC/Z3AE+a0PeSenHM1CY2zewh65q5ejXTVKJZN8emmT0cM+Uj0KTFjHFM9P3AJ4FLgdnJhXO/xaZe7XPZklwAXACwAX/ZaHzunNvFrr3rqD3NLS5p15+vPVPs2ruOO+d29RidNGAjnJ27qapeMdFI9rUDOGnB/lbg1oUN2oeobgfYkmMH9HeJJGliJvwUl5VabhL9QJKzq+qDE43mAZ8FTk1yMvBN4Dzg36/R95bWxBW7j+b2nZuZvqfpxp3Z1fx1XTPT3L5zM1fsPrrP8CQtw3KT6EuBVyfZDezhgRWLtkwiqKram+Q3gQ/T3OJyUVVdP4nvJfXlku88nru/fSQb72yS57qdzfHMhruPPJJLvjP/IIn7+glQGqqxVaJVtXnSgSzyPT8IrFXlK0kaidFMLEqyrapuPsD5ACdW1Y7VDkx6qLvs66ey4bYZNn672V93T/ObYea+UNMzXHbsqQCcxLV9hSgN01iSKPD6JFM0s3OvAr4NbAB+CHg6cCbwezQTgSStwN6vHMWWW2HjHc0tLjO7monv6+6ZJnNTfHfDUX2GJ2kZDphEq+p57SPQng+8EDgB2AXcQNPV+pqqcsBG6mDzzeHI2/ay/s49wAMP5V539wxTe9ZRU8udsiAdZkZUidIutXfhGsQiSdIBpUY0JippcjbfspcNt9/L1M7mmQ7Z01SiU+tmmNq9nsxt7DM8abhGuNiCpFW28bZdTN11D9zbjojsbZJoZmaY/t4eNs42f277AF9pPwOqRKf6DkCSpLFaVhJN8rEkZ+93bPtkQpIOD1P/shN23k3turfZ7r2v2XbdCzvvZupfdjZtJO1jfly0y7balluJngy8IsnvLTh2+uqHIx0eZrY9Eu65p0mau3c32/f2NNvu3dS998E998A99zRtJT2gDmFbZctNonfR3BP68CR/l+T7Vj8USZIO4hCq0ElUosudWJSq2gv8RpJfAT4FHLP64UiHiV27qPua6pNqnyc61/yEp3L/a4BmYTBJQ7TcJPpn8y+q6q1JrgVePJmQpIe+2v09au9eqDlqdt9H9NYsZJrmPMDu7/UQoTRgA5qdu9wF6P98v/2raFYwkiRpbY0tiUpaXbVnD8zO7tNte3937lTTnRtmH2gr6X5DWrHI+0SlNTZ91FEwO9ckzVpiKYVqztdcwexc8x5Jg2MSlSSpI7tzpTVWs7P7VKALu3Tn9zOVhQceNPlIOqwNqDvXJCqttfmx0LbLdjFNIp2/9WXq/vFR6bDnU1wkSToEJlHp8HWwKnSfdkCm5qg5py9I9xtQEvUnU1ojmZ4m09NLJ9CaW3S27nzSzfT0GkQpaSWsRKU18qDEeYDbW/aRqYNWrdLhIjgmKklSdyZR6TC0VOU5qfdJD0UDm53rmKgkSR1ZiUqSxmVAlahJVJI0LiZRSZK6GdKYqElUkjQuA0qiTiySJKkjK1FJ0ngUg6pETaKSpFFxTFSSpK5MopIkdTOkStSJRZIkdWQSlSSNSx3CdhBJTkry8SQ3JLk+yUsP1N7uXEnSeEx+du5e4Heq6uokm4Grkny0qr64WGMrUUnSaOQQt4Opqtuq6ur29U7gBuDEpdqbRCVJWkSSbcDjgSuWamN3riRpXA6tO/e4JFcu2N9eVdv3b5TkKOA9wMuq6rtLfZhJVJI0Kod4i8sdVXX6AT8/WUeTQN9RVe89UFuTqCRpXCY4sShJgLcAN1TVGw7W3jFRSdK4TPAWF+CpwAuAZyT5fLudvVRjK1FJklpV9SmWN5EXMIlKksakhrXsn0lUkjQuA0qivYyJJnl9ki8l+UKS9yU5esG5VyW5KcmXk/xsH/FJkoYr1X1bbX1NLPoo8KNV9Vjgn4BXASR5DHAe8CPAWcCfJpnuKUZJ0hBNdmLRivSSRKvqI1W1t929HNjavj4XeFdV7a6qrwE3AWf0EaMkSQczhFtcXgh8qH19InDLgnM7WGLNwiQXJLkyyZV72D3hECVJQzGk7tyJTSxKcinwiEVOXVhV72/bXEizYv475t+2SPtF/9ntMk3bAbbk2AENM0uSJmbyT3FZkYkl0ap65oHOJzkfeC5wZlXNX5IdwEkLmm0Fbp1MhJKkURpQEu1rdu5ZwCuAc6pq14JTlwDnJVmf5GTgVOAzfcQoSdLB9HWf6B8D64GPNssUcnlV/VpVXZ/k3cAXabp5X1xVsz3FKEkamOBiC1TVDx3g3GuA16xhOJKkMTnck6gkSV2lhpNFTaKSpPEY2OzcIdwnKknSKFmJSpJG5bCfWCRJUmcmUUmSurESlSSpqwElUScWSZLUkZWoJGk8JvQ0lq5MopKkcTGJSpK0ckNbO9cxUUmSOrISlSSNi2vnSpLUzZC6c02ikqTxGNgC9CZRSdKoZK7vCB7gxCJJkjqyEpUkjYvduZIkdePEIkmSuii8xUWSpK6GVIk6sUiSpI6sRCVJ4zKgStQkKkkajaEtQG8SlSSNR9WgJhY5JipJUkdWopKkUbE7V5KkrkyikiR1YyUqSVIXBcwNJ4s6sUiSpI6sRCVJ4zKcQtQkKkkaF8dEJUnqakCLLZhEJUmjMqRK1IlFkiR1ZCUqSRqPwolFkiR10TzFZThZ1CQqSRqXub4DeIBjopIkdWQSlSSNSqo6bwf97OSiJLcnuW45sZhEJUnjUYe4HdxbgbOWG45jopKkEamJLrZQVZ9Ism257U2ikqRROcTFFo5LcuWC/e1Vtb3rh5lEJUmHkzuq6vTV+rBex0ST/G6SSnJcu58kf5TkpiRfSHJan/FJkgaoqvu2ynqrRJOcBDwL+MaCw88BTm23JwH/u/0qSRIUxPtEAXgj8HL2nS91LvD2alwOHJ3khF6ikyQN0wQr0SR/BXwaeHSSHUledKD2vVSiSc4BvllV1yRZeOpE4JYF+zvaY7ct8hkXABcAbGDT5IKVJB02quoXV9J+Ykk0yaXAIxY5dSHwauDZi71tkWOL/unQzqbaDrAlxw5nIUVJ0mQN6Df+xJJoVT1zseNJfgw4GZivQrcCVyc5g6byPGlB863ArZOKUZI0PkNagH7Nx0Sr6tqqOr6qtlXVNprEeVpVfQu4BPjldpbuk4F/raoHdeVKkg5jzs5d0geBs4GbgF3Ar/YbjiRpUIpBPcWl9yTaVqPzrwt4cX/RSJK0fL0nUUmSliss72ksa8UkKkkaF5OoJEkdmUQlSepgYBOLfCi3JEkdWYlKkkbFiUWSJHVlEpUkqYvJrDzUlWOikiR1ZCUqSRqPYlCVqElUkjQuA7rFxSQqSRoVZ+dKktTVgJKoE4skSerISlSSNB4FzA2nEjWJSpJGZFj3iZpEJUnjYhKVJKmjASVRJxZJktSRlagkaTycWCRJUlcFNZwli0yikqRxcUxUkqTxsxKVJI2HY6KSJB2CAXXnmkQlSeNiEpUkqYthLfvnxCJJkjqyEpUkjUcBc94nKklSNwPqzjWJSpLGxSQqSVIXNaj7RJ1YJElSR1aikqTxKCgXoJckqaMBdeeaRCVJ4zKgiUWOiUqS1JGVqCRpPKpcbEGSpM4G1J1rEpUkjUpZiUqS1IVPcZEk6SHBSlSSNB6F94lKktSZKxZJkrRyBdSAKtHexkSTvCTJl5Ncn+QPFhx/VZKb2nM/21d8kqQBqmoq0a7bMiQ5q81BNyV55YHa9lKJJnk6cC7w2KraneT49vhjgPOAHwF+ALg0yaOqaraPOCVJh5ck08CfAM8CdgCfTXJJVX1xsfZ9VaK/Dry2qnYDVNXt7fFzgXdV1e6q+hpwE3BGTzFKkgao5qrztgxnADdV1Ver6nvAu2hy06L6SqKPAn4qyRVJ/jHJE9vjJwK3LGi3oz0mSVJjst25K8pDE+vOTXIp8IhFTl3Yft9jgCcDTwTeneQUIIu0X/RPhyQXABe0u7svrYuvO+Sg19ZxwB19B7ECY4sXjHktjC1eMOa18OhJffBO7vzwpXXxcYfwERuSXLlgf3tVbV+wv+w8BBNMolX1zKXOJfl14L1VVcBnkszR/CfZAZy0oOlW4NYlPn87sL39vCur6vTVin0tjC3mscULxrwWxhYvGPNa2C9JraqqOmtSn91adh6C/rpz/xZ4BkCSRwFH0PyVdQlwXpL1SU4GTgU+01OMkqTDz2eBU5OcnOQImsmulyzVuK/7RC8CLkpyHfA94Py2Kr0+ybuBLwJ7gRc7M1eStFaqam+S3wQ+DEwDF1XV9Uu17yWJtjOefmmJc68BXrPCj9x+8CaDM7aYxxYvGPNaGFu8YMxrYWzx7qOqPgh8cDltUwNaDV+SpDHxKS6SJHU06iSa5HFJLk/y+SRXJjmjPZ4kf9Qu2fSFJKf1Heu8sS53mOR3k1SS49r9IV/j1yf5UhvX+5IcveDcIK/zSpYZ60uSk5J8PMkN7f/fl7bHj03y0SQ3tl+P6TvWhZJMJ/lckg+0+ye396jfmOSv28kjg5Hk6CQXt/+Hb0jylBFc499q/09cl+SvkmwY+nVeNVU12g34CPCc9vXZwGULXn+I5n6fJwNX9B1rG9fTgUuB9e3+8e3XxwDXAOuBk4GvANN9x7sg7pNoBtm/Dhw35GvcxvZsYKZ9/TrgdUO+zjSTF74CnEIzU/0a4DF9x7VInCcAp7WvNwP/1F7TPwBe2R5/5fz1HsoG/DbwTuAD7f67gfPa138G/HrfMe4X79uA/9C+PgI4esjXmGYhgq8BGxdc318Z+nVerW3UlSjNDbBb2tffxwP38pwLvL0alwNHJzmhjwD3M9blDt8IvJx9bzge6jWmqj5SVXvb3ctp7vOC4V7nFS0z1pequq2qrm5f7wRuoPkFei7NL37arz/fT4QPlmQr8HPAm9v90Nxed3HbZGjxbgF+GngLNJMwq+ouBnyNWzPAxiQzwCbgNgZ8nVfT2JPoy4DXJ7kF+J/Aq9rjQ10+cHTLHSY5B/hmVV2z36nBxryfF9JUzDDcmIca15KSbAMeD1wBPLyqboMm0QLH9xfZg7yJ5g/A+fXevh+4a8EfWUO71qcA3wb+ou2CfnOSIxnwNa6qb9L8/v0GTfL8V+Aqhn2dV83gnyd6kOUDzwR+q6rek+Tf0fz19kxWuGzTapr0coeTcJCYX03TPfqgty1ybBAxV9X72zYX0txv/I75ty3SfgjT04ca16KSHAW8B3hZVX23Ke6GJ8lzgdur6qokT5s/vEjTIV3rGeA04CVVdUWSP6Tpvh2sdnz2XJohkruAvwGes0jTIV3nVTP4JFoHXj7w7cBL292/oe2yYYXLNq2mg8R7yMsdTsJSMSf5MZofjGvaX5RbgavbCVyDjHlekvOB5wJnttcbeo75AIYa14MkWUeTQN9RVe9tD/9zkhOq6ra2S//2pT9hTT0VOCfJ2cAGmqGfN9EMPcy0VdLQrvUOYEdVXdHuX0yTRId6jaEpXL5WVd8GSPJe4CcY9nVeNWPvzr0V+Jn29TOAG9vXlwC/3M4gfTLwr/NdIT0b1XKHVXVtVR1fVduqahvND/hpVfUthnuNSXIW8ArgnKrateDUIK8zK1xmrC/teOJbgBuq6g0LTl0CnN++Ph94/1rHtpiqelVVbW3/754H/ENVPR/4OPALbbPBxAvQ/mzdkmR+AfczaVZwG+Q1bn0DeHKSTe3/kfmYB3udV1XfM5sOZQN+kqbv/RqasZkntMdD81DVrwDXAqf3HWsb1xHAXwLXAVcDz1hw7sI23i/Tzjge2gbczAOzcwd5jdvYbqIZY/x8u/3Z0K8zzWznf2pju7DveJaI8SdpuuS+sODank0zzvgxmj9iPwYc23esi8T+NB6YnXsKzR9PN9H0YK3vO779Yn0ccGV7nf+WZgho0NcY+K/Al9rfbf+HZgb8oK/zam2uWCRJUkdj786VJKk3JlFJkjoyiUqS1JFJVJKkjkyikiR1ZBKVJKkjk6i0xpJsbNdOnl6Fz3pYkr9fjbgkrZxJVFp7L6RZ/nH2UD+omqXWbkvy1EMPS9JKmUSlVZLkiWkeBL4hyZHtQ4p/dJGmz6ddAi3J05JctuAhzO9ol04jyc1J/keST6d56PxpST6c5CtJfm3B5/1t+5mS1tjgF6CXxqKqPpvkEuD3gY3AX1bVdQvbtGvjnlJVNy84/HjgR2jWgv5/NAunf6o9d0tVPSXJG4G3tuc2ANfTPOgYmiXifn8S/yZJB2YSlVbXf6PcMTscAAAA90lEQVRZUP4+4D8tcv44msdFLfSZqtoBkOTzwDYeSKLzC9FfCxxVzcOwdya5L8nR1Tyw+XbgB1b1XyFpWezOlVbXscBRwGaainF/9y5yfPeC17Ps+8ft/Lm5/drNLWi3of1cSWvMJCqtru3Af6F5EPjr9j9ZVXcC00kWS7BdPYrm6RmS1phJVFolSX4Z2FtV7wReCzwxyTMWafoRmseKrZanA/93FT9P0jL5KDRpjSV5PPDbVfWCVfq8TwDntlWupDVkJSqtsar6HPDx1VpsAXiDCVTqh5WoJEkdWYlKktSRSVSSpI5MopIkdWQSlSSpI5OoJEkd/X8mavltqy8PCAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 15.6 s, sys: 3.3 s, total: 18.9 s\n",
      "Wall time: 21 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "xy_dimer3 = NearField()\n",
    "xy_dimer3.coef_filename = 'amn_test4_dimer.dat'\n",
    "\n",
    "xy_dimer3.name = 'dimer_ex3'\n",
    "xy_dimer3.output = xy_dimer3.name + '.out'\n",
    "\n",
    "xy_dimer3.plane = 'xy'\n",
    "xy_dimer3.nx = 300\n",
    "xy_dimer3.ny = 225\n",
    "xy_dimer3.nz = 1\n",
    "xy_dimer3.xlim = np.array([-80.0,80.0])\n",
    "xy_dimer3.ylim = np.array([-60.0,60.0])\n",
    "xy_dimer3.zlim = np.array([0.0, 0.0])\n",
    "\n",
    "\n",
    "xy_dimer3.save_int  = 0\n",
    "xy_dimer3.save_etot = 0\n",
    "\n",
    "xy_dimer3.run()\n",
    "\n",
    "\n",
    "aspect = (xy_dimer3.xlim[1] - xy_dimer3.xlim[0] )/ (xy_dimer3.ylim[1] - xy_dimer3.ylim[0])\n",
    "size = 6.0\n",
    "fig, ax = plt.subplots(1,1, figsize=(aspect*size, size));\n",
    "\n",
    "x = xy_dimer3.grid_x.reshape(xy_dimer3.ny, xy_dimer3.nx)\n",
    "y = xy_dimer3.grid_y.reshape(xy_dimer3.ny, xy_dimer3.nx)\n",
    "z = xy_dimer3.Ecartint.reshape(xy_dimer3.ny, xy_dimer3.nx)\n",
    "\n",
    "pcm = ax.pcolormesh(x, y, z)\n",
    "#ax.axis('square')\n",
    "#plt.axis(np.r_[xy_dimer3.xlim, xy_dimer3.ylim])\n",
    "plt.xlabel('x (nm)')\n",
    "plt.ylabel('z (nm)')\n",
    "plt.colorbar(pcm)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-26T18:07:41.694070Z",
     "start_time": "2018-03-26T18:07:41.149390Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdcAAAF3CAYAAAD3mAEvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvWuwLUd15/lfVfuce7kIvbgSEpIMwha4sQwNVktgwDyEDJYZ0x/sCRxtB8buUbTH9rjb7bDBxAQxHxxtt2dwMzHd03MD40c0NqZp90AztLGQzUOgh4V4CYR5C11ACAECpKt7z9mVaz5UZe2s2llVmVVZz71+ETvO2fVYmVWVO1etlStXEjNDEARBEIRwRGNXQBAEQRCWhihXQRAEQQiMKFdBEARBCIwoV0EQBEEIjChXQRAEQQiMKFdBEARBCMyoypWIziWitxHRp4nobiJ6NhGdT0Q3EtFns7/njVlHQRAEQfBlbMv1DQD+mpl/EMDTAdwN4NUAbmLmKwDclH0XBEEQhNlAYyWRIKKzAXwMwJPYqAQR/QOAFzDz14joYgDvZeanjFJJQRAEQWjBmJbrkwB8A8AfE9FHiOiNRPRoAI9j5q8BQPb3whHrKAiCIAjerEYu+5kAfo2ZbyOiN8DDBUxENwC4AQDoyP6P7F18QT+1zAvsV7zQF5Les1/khzFLev5ZHHzpKw8wcy+d8kte+Gj+5reS1ud/+ONn3s3MLw1YJStjKteTAE4y823Z97chVa5fJ6KLDbfw/baTmfkEgBMAcOTyS/ni/+1Xtw8K+Lsnkk66N3rsnwd5blPVLwNcOnOPFy8/ud4I+twsz+meV77mnnAFFHngWwlue/elrc/fu/jzxwNWp5LR3MLMfB+Ae4lIj6deC+BTAN4B4JXZtlcCeLuzUCp9WkLEW5/RKF9T188U6xyq2D6e21TuaRsGqHevv5UltO0p1NlC0Oc2kWuaGmNargDwawDeTET7AL4A4FVIFf5bieiXAHwZwM84SeqoTHtlSg1uSnVpSS/PawH3xZuqa+54e6ueT69WLjDPZ9i2zn38BIzn1ulZ9f4cGAmrvgvpzKjKlZk/CuAqy65rh66LIAiCMH0YgJrBmMHYlusoiOUzD5bwnKjn8nqZSbc0i3aJNN2yHp7VlJ6TgliukyFYRz2d9jV7enXH9/Sc+laWvvjWp5MyHljp5uIn1KnPhh6elX5O8jzcWLRy7dx5j9SGJDLZk8DPaWoKNCR119Za8fakdHPxO/B7GExhlYtpcWvLz2NoZctgJCMlP/Jh0cpVEARBWB4y5joCU7ZWd+ENvFeCTtsJJ6sTdW1iYIsguFXbs0W7JHz6hqCWYkBLdigLlgEkM2hEy1Cu1EFxBZ5nKXSgh99m70q013HjALIDdXhV9zGo0gVE8TpQ1890VnDm6Z7PYsj+bw6W69ir4giCIAjC4liG5epL15c7sVD9GMBbFNxCXcozdr2OlhZPUIsW8GsrC3lEIQk63SmAy7gPGJhFQNNuWK4d03MNkgoxdFq1KX0CQmT/tBPG1Z9dI/C9CPqcKgtZ8CcwQdId9lxHH1SHTxNE9CYiup+I7rLs+00iYiJqzE+8G8pVEARBWAQMRtLh48CfANhaNYeILgNwHdK0vI0s2y3c8u0quHU68lve3JiTi3fqQwRBIjirrrGF7F7m2S6RplsbYk2KrvNV9eFDPzcGkh7LZOb3E9ETLbv+EMBvwXExmeUp1xZ9SbDVU4Raeo3cDazkpq40XXG9jnZjcmGnEbm2D1HC6EX5tk7cP5aSbc9xIrrD+H4iW8K0EiL6KQBfYeaPkWNDXZ5yFQRBEBZLmri/Ew8ws23BGCtEdAzAawH8uE8hy1Guvh6NrusXDsRkkh1MhYAW5VKs0xAEnzsZ0JW8JXpBv4nerPDyPeo4Z9WpDQz2XAjJsK7C7wdwOQBttV4K4E4iujpbl9zKcpSrA1NRqEvqHHphBgp0Ts+wawceVPFOKCPVFPBpR0EXXWipbKeQtJ8BqAHfi5n5EwAu1N+J6EsArmLmB+rOk2hhQRAEQcggor8AcAuApxDRSSL6pTZydsJybT2nq3O53WUshplE7C7xmfVpHQ1m1W4JX+CDaqDXFJSOMloHPQWmT7cwM/9sw/4nushZpHIdI8/wIJ3yjo4RLkF5hryGPju1kFNl+s2BO9PfQg/PzvbMvBWuKWPiijZN3D/9l6tFKVfvDqz1PNh259mFzbSTCEjocdE+FOiUgp9C1KVNZziU4t3InH4H6s1A1nn5WXk9nxZW7dCKVs2gbSxKuQqCIAjLRizXQfHMldkq0YT/OZuTp2P19M1QFl4o63RKFumQhE4u0fQ82owLtnk2i7J2A0ZWB7FkW7iLd5mFKFcHvOfBdimrr+kfu9NoZ6c8p9anB7rspvs3pvK1lzPt30gw5W+7Tg/ZrZTtRDIxMQjJDCa67I5yFQRBEBaBjLlOAZ8UmW2eV6A35am/cYdiFsFG0//dNtPmGgK7bn2stF3JKzx4rucWnoXGe9wisjgkMuY6Jn0q1I4d+dSU6BzmdQa7ZzO41lFxuT8+oQ09TMUZur2OpcyDLXpeluNwvr7H03UVExKevlt4+jUUBEEQhJmxLMu1L4u1heU0h7mbUyDofZqDy3lCtHM9ugh2ENOnezQgff3u2lrEnS1a8/yGc6bqKk5XxZm+Xbgc5drQtvp2/3bthJegPHtXREFSUi5XWfriey/cO/AmQe5lLnU6Tujo6Xar2AygaHtCxlwFQRAEISDM8xhz3Qnl6mwVDmCt9m2hztYym8m81um/L9fTadWyUCkLAwdObYkf8TcQymrumqTfOx1hCyt27pHbfbNo5eqkyDx+iL4/2tklQgjNxMdA564o29D3DJ1QSSecKzqxn4ZP+wyV37lJyXm7jPXxE1ayaga/3kUrV0EQBGFZpPNcxS08OCFdwENaqoNbpxN88QtxD8a4rKl6FoK5KH3KbJIVyrLNBfod3oqeHm8fCThcrEhnl7Gjq3h4C1bGXAclpAvYp7Nso1B3LSnCVF25U1WKoQh1fV4dfZOspvN7dqu2omsxgTJfuVyvb85gb0XroGT7RqbiTImAStW3AY2xcHsIxlA8IS55SgpzqLr0rWRCKryAs3RmswB9qItuo3B9lK2WH8KaFXZFuQqCIAiLIZmBYl+2cg04rupjsfa9tmywsgdiiq7cKd6ntgxxLe5Ly3UbUw05vuvDGC70TeFVwhxO9Zxy4zI+6mTBAs5RxaGRJefGJJBSdVWoQynTsRTCED+dIMFMC1KYUyOUazhk2sMRhkIbCRqwZTvUQSm6lOWjZJtkgXhwBaskoEkQBEEQwiFTccag4a0xlLXap6U61cjaLkzJlbuEHM5j0SY7kF2Oi/syUDaoujI6nW3Uw6fMrvemvLuj5eka8NToKpZApy2Wo1xrGm1Tgw7u/u3Dndy+mCD06XIVpTkPfO6xq7uxXkZTJGy7dhM6MrjPCGgvZRvAxdvkLnaS0/PwDIMkoEkQBEEQQiPzXIeibqw9gNXa+IbcVySxn+ig5U6lvBAWaSSBTr2gAiQUGMLCbSs3VJl1R/om1agPLHIT7hINTOSWgGKM5f2YIRmaxiLE2KrTDzDw+GzbZtr7SjA9ye+iOIdQmEN3wm0Ys44+z6BKEbu0Ad/E9K50eTahxoQ9A4LbKduWLl6XtIrO03Z2kEUqV0EQBGGpkKyKMwZ1b5RDuYCdLOdmMd4yh5DhVk53GV0s0ynNd51SXaoYwvptep59u5jtcvu1eNsm5vexZt0CjJoFNVmgLoFOQ1mvDHELO0FEMYA7AHyFmV9GRJcDeAuA8wHcCeDnmfmgWU77sdUQLuCQCrXdIuzTc93a6OrOndK0nl0iVORulRzXdhFCCRfr439OH1ONfDJXBVG0HZSsSyRx38xhnusUavjrAO42vv8+gD9k5isAfBvAL41SK0EQBGFyMAiK23+GYlTlSkSXAvhJAG/MvhOAFwF4W3bInwL4p41yasuofqMl4vq3LTI+NedXySA0iijIaKxPxfG+S+T5fHyIiBs/vtflcp2u57a5X0JKqHvc9fl0aWP2+oT/PbRpcz7H+/YrtUJanN+2f9g1xnYL/zsAvwXgMdn3xwJ4kJnX2feTAC6xnUhENwC4AQBWx8+xCm/tCu7oQm5qc37K0L+zCEHIaNwuyqyrIpQOoF+6TqWpdj+2P7dN23W1aJrak+/98MnDXHmv9H4HOVsyHMdj27iJ+2QObuHRlCsRvQzA/cz8YSJ6gd5sOdT66Jj5BIATAHD0+y/Jj+k0ttplvmy9ZDcZnhaoL31MXxlDcYrCnA5dszb1qYB9AmxCTCtK61RVl6rj3RVu43ioKaOyfjUyarR03XkuiwCEhCGJ+5t4DoCfIqLrARwFcDZSS/ZcIlpl1uulAL46Yh0FQRCESUFIZCpONcz8GgCvAYDMcv1NZv5nRPSfAfw00ojhVwJ4u6vMPtzATefWj/eGsVTbjH22YSqRuFOxyoX22Ky7ti7VtpHIfU21qWprvhZt2/y9Lu5iF1dxbURxSxdxWp+KQgPRt+VKRG8CoD2rV2bb/gDA/wDgAMDnAbyKmR+skzNF2/q3AfwGEX0O6RjsH7mc1CpoqWWwkmswga/cdL97EEXooCEf/ANZ2geNDBEoFfozNFO7zjbPq2tg0RBtu801V5fTfF1uv62afWjup6zn1nRyLr/1mfMnAF5a2nYjgCuZ+WkAPoPMMKxj7IAmAAAzvxfAe7P/vwDg6jHrIwiCIEyXPt3CzPx+InpiadvfGF9vRepdrWUSyrUv2riC+3AB1++rEWrg4/Zs86YeCpfrmUoiiT6ZQx1D4HOdpjvRtw2UXa4hAoea8A2k8oledsm1XBcA1uQu3nKRm/W0ltsQ6NTCTdwXzNTVLXyciO4wvp/IAmRd+UUAf9l00CKVa60buMV5QytV147Hp8Po2tkPNe6bluV3buS1YuY08M2NuoRrbHqudZ10XXsyFVWIhQA2sqoPdBkLrTu2fD1NY9TlOtcpUnO/73hspbJsiCQeWsF2TH/4ADNf1eZEInotgDWANzcdu0jlKgiCIAghIaJXIg10upa5+RVtUco1pBu4jyjgprfqpsAcV1yODRl00GSldrWafa22oQMq2kRHDmGJDnkfbPeg7hptlrtvO9HWko/bNS3Hqxjn+bk+buTyseY1tLVi21qwQNEgbeMibrKkQ8Lw9/x0hYheijTY9vnMfMrlnMUo1zauYF/F2sYF3Nb966YgfRRu8zGhprPY6hVSmbRVGi7X1yb36AKiI3Parl7j647t2h4UaKuddRnTLch2GN/1SYjhOz7aVtHWjrk2uZGxrTMrlWzDAG7/8QbU66o4RPQXAF6AdGz2JIDXIY0OPgLgxjRLL25l5n9RJ2cxylUQBEFYPuk8116jhX/WstlpSqjJspXrAK7gkBZrW3ezT9l15bfFx1Id0h09FZl9dQRjJM7oEuUb4tlrS81sX9pFGMo16ZIkYqikEK5WrK+r2CfQqW0kcZ9IbuGxCKRUQypUwP6j7apQuyrS0C4cs9OrSghh4qN4Qo1Jh6brQuF9MYV74jseWmarvZR+rczNirYLda7mNu7jOkVaPqbKXRxCyTaNxXpHEs8vmL13lqlcBUEQhEXCGHZd1rYsT7mOYLWGdAF3TTjR1lJ1CTCpitArn6vr2WSpM1Mrq67yGnq01OqjT+f32t73vaqLkPVfc3VzfNkZaLZJ04p1batVVFnATUFTTdasjyVb5U6ucxU3BTu5WrBTt16VuIUHxkOx9u0KDqVUuyab8FGmTVlvyucoUK1iLZe9VV/iSqWlf8BtFGnbrEF1hIw+DcUQ7mbXe+midOoUrvNzLrWZCJv7a64tGoFzBek7/l/lai7X0yU5BFA/VtvVjdvkKu4i33uqTtWgbWCYgUQs1wGxvUz1ZK2OpVC7KNNcRoUyrCvDNs5VJ0srVlNWVb1i401YVby1V9Wv7lpdx3ZDWp0hF+/uwhCWdNl6apqGUqUE6l7AnNoMNtZsWo/m9mkry6xb2fpzUbSb4+sVbpWl2ZeSrZPvG+xUvlqnqTo9MQe38PRta0EQBEGYGcuxXA1CWKw+VmUIi7WNtepjodRF8TZZC6ZloM+vmnZgWq1alinfdNdpVKmcfHvFGJ1NbhOhrUe7G26aY699W861eX+N/80x9oJrt6L9ldtNKi87L3MPRyi2n6YkF3Vt3db+qqzYqjHdJmvWxVrWZZXldRmPtVmwdbKDRBH3RBrQNH27cHHK1Uex9qX8hnYvF85tGPCoUqx1Cqus/MwfrW3ctewWjlAcf41R/FFHlCrcqo54U/diHavH0rr/2JvSq3Vd1aULvczv7TBQVnu/jc36+fo823Kb122HmQAquoUTx3Zpfjextb+ygqpyOZvllnHJxuTizm3jKm7rJq5UvJiGgu1zyblQLE65CoIgCMuFMY8x10Up16GsVh/Lsi8ruI2lYZZrs1ibAkhMCzZhyt/S9d84s1CJGCtSBdmR8X/B8tVuMSZEhjwFeyRpRFxp1RTuX4VFUKYyCClA7tvKuvVAXwsBNAbI6f0V7sjIPJ02wU9lK0k/1y0LtvRXy1ZMeZtZZ7YrEQOWdhmBC96UunYOYKutV1mv1vvRYNFWTb1pcuemx+pt7u7cOsu4bZBTVVBwXWKKXWQxyrWLYu06vtplzNZVXl3n2XWtVVt079a0G0o7MwBIVAQiAjj9IW0pP2LEkcrlrSJVlFdwFaYTwhUT1irtJCNG7vIr11N3wFGpzLprtSlYZVHcXdjq1EZOWRN6/Nd2j6vue6FDN/7Pn6utzWTPdRWpzTbjRSonu80KaXtRTFAq3cjESLLDzHZJlN4P3SZtz8Z8GdKu5vKYf9X4aBV1CxWUla1t3NRVIZrnN43FNr0o1LmefebC9ouMuQqCIAhCcIZecq4Ny1CuzbEU6feB3LYugUttLdZQie/L7uAqt1x+vOG21REknESFwAJtse5FCVaRwopSK2SlLQbL27O2WpUR5LJWUV6eGZiS19NSx7q5hCZVEchtaWMBt3WbtbEIQrqjbR6Pqn36eZr323Tlmi5b/VzNNhPVPOu8vYCxLlkwiYrywCbdbuJIIc5kl70eGmLaPBcjUMocnvClLmFEOerYZnHa2rTPvNUqN3GT29lHXpX12ieSRGJEXMdYXRVrVWRjCFm2812zJzV1nK6D/lqxmq5c2xhp/gNXm44vURuFvBcn2IsS7McJ9jMFCwArUoiMTrOoXFPX3pojrDg9Zk0R1tl2U+nmdS259rajPktRo1mZZVl196d2eknN1Iq6YzUh3WahpiO53out9gH7y5h2/5puXP0MzPL081yRyl/CAN1muNBuzPMUR1hHEdYqwgHFm3JXACUxmAlx1v7iSCEuydcyzfrk90BFALl34K5twKbIbFN7XMZOXcdiu7idbedORcHOwS08/RoKgiAIwsxYnOUa0mrtErzU1gKucwP7RBPb5unZoJLVmgeUoOiWM+WsoygPQIquvXdL5jr7hCACcPR9F6X/W6xqE5sVEoELFpR5bFuXqct5od3PTXQtp+l8m8VqbnNx0ZvP5/Tz79tsB3CQfboQATiSfQBA3XQZgM2wxKpkBeu6KCYo2gTUIVKAiqCIAUeLvul3aJaXHq/367rb3cSuc1dd562WLdhQsoa0XmVVnBFoq1hDRQX3ObZap5htSkb/IEzlqH806ThmUVZ5zGtViqx88LkPWOo8DGZHrDn+ofO2tmkXczpmlt7AvMNkbFyV2O5YQo/Buspr65oOcXwV5ToVlKjhBtZtpuyi165ckwd+9NtB6uaDfvFT2Ue/8J178/H8GO2uXqsob9DrXLFm03Ky+7Fxv9oVa9XvMD23OvK4r7FYl3FTnyjiJln6vyGUrAQ0CYIgCEJAGJJEYlBcrNYuQUdkeVPtKsclcKnK9VRryWZ/y5GOCmlEJFnSwq0itRWE9K3nfHP7oiZC2RJ6/K1nA0jy4Kg84CHCxjLJrFeg2i3cxfprcjXbOoQ+rNM+rkFbrXnAkeFqjUhtuVy/+qzvtq5D35S9MOd/8LHp9alti1Ax5YFNeu6sTsFYtlirrFSgmPrR3GfLXwz4BTrV5RB2DUqyyWkbkTx9tTcMi1GuJqGjeX0Va5ux1bKcKhku7ihbUnNkUw3WSZS6geMkP1x3kqd+7OtbdZoLZmd+6W1nYa3SCNKIo/T/CHnSARuFMcMWji0dlWzd19P4a528LtdQd4/Ml69VlOR/T17zkHd5U6H8Enn0fRdhTenLGTOwTvT/hFWsAComMtma3lNytZrxDzYFXKfUXMZiXVy7tmQTvoq66ryqzE59Modo4UUqV0EQBGGhsAQ0jUIoq7WNHJfgJZekEK5Wa3muYVmGmbJQv1kmKgIz0jdwpBYIXnQSp7erMVtMKyp1FwNKpdZDxEVLoG6+ZlU0chmd0GBrO/qxWMtUWalV5VZdly2y2txnRgWvonQu85Tdv23RAXSrG5+AR0A4XKdeEKL0dxNRaU54lZdAD0FgO3rftPZcLEZXC9bHtVsVLBVq/mtfMCSgaXBCTbkJERXc5Ar2dUmXx3XMyM26rDN5VqW4uO/8n/zMdgUWiO780+jiVUFZFP6a97VGKbm+MVcp3K3jajoJH9eu6/hx3bb82sxDaHNMOi6fxtyOEf07NHvX3YPzAXzznU/Jt+3HCVZxsrUwhUnexvRmKipYfUxZwZbPrxqLLY/DAvXRyOb5+pw6127IKTp9IZbryLSZchNasbbNtmQLXrIp0aqUgOUfSKwi7EUJHvMTX94ufAfQyuDsmy8AAESKC8E7tnmcwLbya1KwLuOrhSlAHRVo3f7G70bZ5UxWZl21Yv3uc7/hXNcl8diX/UP+//o931cYewaKbaJw71CyCMvfeXuFoLI8TdVasqYl22TFAsVApz4UrLBh0cpVEARBWBYMsVwHpSlfaxur1dfy9R1fdRmjraKc8Ny0uGyuzVMv/LrkugQK1peZTKA41ppaJYqjopWgf9ClMUnTFdcUMVwuqw0+FqwtCnrbot1cr5loQ2NLILKrrF6cen6Ovv9xhe1mpid9/9YqShcA4Hr3f12iCJexWFfrVcvwsV5tdahbz7Z8PX0hynUgyvM2XR5uUwcVUrF2Uaq10y0MxWrOM9RZcyLiSc9VHRutNB5/69nWH6sqjYNFlCkgy5hklWLtOw1iU9BSVfYk69QjY98SA5VCcurHvo7zP/hYAJliNdaX1Yu3ryKFtbFKTxm78rOPoepy0mPqFWz5/CYFax7vWs9ynYd0DUv6Q0EQBEHoAYkWHoEQ0cFNVmvT+a5Waxf3iRnIpJfrAjaT+3chmjMUX33Wd3H57cfy7/qHmzDlk9VVvgapXukzw4gINRk7cb99atEmX7RpxcYlK/eLV58aosqzx/QKHf/QeYgQp+vLZrFOa0TZlDiuXQDARtNUG1vi/6psTq7Ta1wTTLhGIO86i1KuIeahhowM1ue7RCRX1c8F0y0sirUdpkL5wTv2ccgRIkS5azhhSjM9EQCoXOnqTqccXezyLF1dW67totz+NhHQxgLkYKyiBDGZSfgV9kjh01d1XZdmd3ngR7+dTvdSgMq1T/N5dUMJwPZ4rG2qjauC3ZbdXsFWna/r2SssY66CIAiCEBSJFp4ANsuzuL8fq9UlR7BtXxWbxAebXME6OURqiaSuYGA3Jvf3zaevOsCVd8ZYqxiHeQ7TGIiSgvUKZBHFxhu9T+IIX09Fmzmx2hWsz9dW64oS7BnejruemVTKEtzQ1mt5uUO9bJ25fB3QHIwGbFuxdYn/mbfnZ3eZ/+piibpED/eBKNcB6aIYfaOLq8ZYXcdmXdzT5QbNnK7OkcuIN7J2dXJ/n9z1zARP/wgDifkTifN1PjcKUuUKto62iSNsOLuJtTvYGH81FeuROM229LFnbK+/KrTDfLk9+r6LAKRtJU07anf32vqfqmjhOmVZ5SJ2maJjYsu81CY94q6zGOUqCIIgLJ+5TMXZibwCrgFFm/3V57exWonYarXq7XX7gI1b6VBFOEhiHCQxkmxOnU4yLoTnY8/IllWLUksvjbJtZ+XpyO4Qn7ZE2TXsZRHlH3uGEqu1R04//750Pdjsd3uY/Wa1ixiw9CfE+fa6fbbt6T4jkBL2/qR8jj6vjvK5Q0XBV8HZPWzzaYKI3kRE9xPRXca284noRiL6bPb3vCY5i1OuLgkbisdXN7JyI2rKwGQrsypHMFV0muVt+li9ruTB4QoHhyscJjGia++tvTahO3c9Mx2PjLJFwcuY60rqY8b+2OpWrqOMsQ5DdO29OExiHByusE7SFamY3fuAqhfwJmWZb69JsFOnYP2HyhoPD4qeGtfm48CfAHhpadurAdzEzFcAuCn7XsvilKsgCIKwXJiNdJMtPs3y+f0AvlXa/HIAf5r9/6cA/mmTnEUpVxfXRl0ggXOUL9gaHVxl6Zrby2+gWpb+xKQQaxdkNtk/jhTiiDOXT3oN51z/2fqLFYLyqR9ZbyI0kSaXKFut6d9w7t+27uKy9Wq+sUfE+NSPrIe4ZULGOdd/NvvdAnHE2WeT0EP/5gu/ezT3Hea27e2b8pus16rzmvrLMqZLeoE8jpm/BgDZ3wubThgtoImILgPwZwAuQjqv4QQzv4GIzgfwlwCeCOBLAP5HZm41v6TOpVvXUFyn3KT77DLqfgimYq1002T/x9jkrY0jhWMv+UJlvYV++exVZ3DpbXtbb8AFd17HSOBOGHmObX+B9BqE4Xncy+/GqXc/CXGks2KprX4GAEDGtC4uTtsyo4WLU3Hs2ZFsy9S1iSCuSy7RFHncFx0TVRwnojuM7yeY+UTHKm0xZrTwGsC/ZuY7iegxAD5MRDcC+AWkvu3fI6JXI/Vt/3aTsC5vTK7n1inWOkVekFGjWMvjLiaKCeqFX3GrqNAbJ695qDCXUQc5acW6eX7147NtqQqoKiwokGeRigtjTCeveahz+UJ7jr3kC4j+7hIA1XNZzRciIs4XXDfX27VhU7rp9uok/1V0nbvav/XaOVr4AWa+yvOcrxPRxcz8NSK6GMD9TSeM5hZm5q8x853Z/98DcDeqki3OAAAgAElEQVSAS9DCty0IgiDsDn1GC1fwDgCvzP5/JYC3N50wiXmuRPREAM8AcBtKvm0isvq2iegGADcAwOqCcwr7mqbe1CWccLFAXSOSK0PnsR39Vxg/A29tf+h5jS9KwkDkeWQzVtpyJZUnwbcRUxqhW7UEWRNNspPc6knHWfWrs2Ttmg7a+3TsAxdugmwqoli19QoUXcRNSR3K26ssUdeEFOVjXZa16xNGvxmaiOgvALwAqfv4JIDXAfg9AG8lol8C8GUAP9MkZ3TlSkRnAfgvAP4lM3+XHH0KmY/8BAA86gcenzeDJsVa3Ff83te0m6q6ld3AupM2Fz9fkZI1WSeIVlhPuP3RTvNfCy7hiiaulWOVEjUT7dtYGXK1gr3n6odr6yWMw0PPux/nf/CxWHO2BqzarL60NhyKZSVijr0Cm2T+TdmbCjIc13+1rX5j4qLk5wgz/2zFrmt95IwaLUxEe0gV65uZ+a+yzV/PfNpw9W0LgiAIO0I2V7jtZyhGU66Umqh/BOBuZn69scvbt+1e5nbWExfq1mc1ZfvKAranbqwihRUp7Edr7EdrsVonzj1XP4xHr87gSLTe+uxlS7lpl7H+rCjZmnIBpBarubaq+VlRYshI24qWbyv70aszYrVOnG8955vYj9b5eszaY1WeHqNxiUSv6uPK02TaRLVPKWdwz0kkgjCmW/g5AH4ewCeI6KPZtt9BC9+2ps4l3HZOa9U55nl16cjq6lqYmwhTwSb5+qzC9LnrmQmu/miCBNEm2hMREtZrvWbfCz9s4/lWZFIqftdTs4ou6ChTyBGUcYzC7f84DnNxQq+sSAERoNTGvRox59HBmqZx+rJb1owerluDtXy8XZZ91RzbeUD/K9YwBlgzNgCjKVdmvhmVo09+vm1BEARBmBKjBzSNjWvCiDpL10VWGZ0YwjxXBzWtSOG+Z3/Hqf7CNDgWH+CQYySZFXrIMSIQFEWZ1aGKFqqR8KEqQMlsH7ERzJRbqIbFukcJ4kzOHiUAHhX2AoVeuO/Z38FFt5wDlT07RVm/oOe3NlhovpHD6TmbtV9NN+l8ApRkVZxBKUfturqEK+U5RAfbZDdtNxuz/YcginWOvPeHH4Uj0SGOGp89SrK1U5OtcdcYnLl4qz/6GO0K1op1L5OpZevyjmSf9/6wKNY5cd+zv4PyogvAdv9QNV7o0gfVjeO6jeUW5db1ry7DYl2ZQ0DTzluugiAIwryQMdeJ0zaQyXa8z9uaAmWTwovpzqreLoV58LdXPhrX3WWkGEwjmQCKcrfuBgXFtu3bbLmCjcClvSjJLdkbrzwrzIUIg2POdzdXcOE8QM5NmWynP2zOO1zcPnyAki+pBTqtOtlYnHL1ycZUdU6TS7jJDVLpEjYbBGVZV0o/nvuf/aC9ksIs0AouV7KZglUggNJRmCQfb8/G2RqihctjrHF2nijW5aB/92fffEG26lKqWHW/4DP2CpSTQTiMy5YSS7go5LqE/sIClasgCIKwbKZmTdtYTECTLy7uV3PitUuEcHmgvxygomEmJCrCoYqRqAiJSudHfve532h5NcLU2MuSPmgiHdmLTf5hcwWd8ic/xxIZbBKRyqKDhSXw3ed+A8rSP5StUlu/4hpYtPHAudVpikNVcwhoWpRyDekS9h1n9YmYy388SYwzyQpnklWeW1RYBu/6oXPSKF/Uj6PXRQs3npNFEr/rh86pPFaYH2sV4UyywmFSfPmuoukF37Z9c7xhQKC6PzOPL36vLrtPRlgVx5ud7NEr3+gcp9+UMVe4sWVesslJFOEgibFOIqyTCNG19zrXX5g/ZqrDpmOajhOWRXTtvVgnEQ6SGIkiJKpotQLbU2iarNitMgL3gUPCaK9YRbkKgiAIwoxZpHLdnuTse361XKD4llhelxUoJl0vf4+IEUcKccTYixKsYoVVLHmEl8iNV56VR36m63amo64mpnVa/pRRWe5icx1QiRJeJqs4nWYVR5x9VD4UULBYy31NQ9/kO92w6thyUomh4Q6foVhMtHCbB+zjOqk71nTV2I4vyGFCTMgT8+9dd49fpYVZkXCEQ45xyJtE+puVb9xeqrRCTpiQYCMnGjI6QxiUs176eRze+IT8e1lRlsdgmakwnaYufWFVUn/bcVq2D70rW5nnKgiCIAg9MIP3ysUrV6eE+x5PqsqatVmtW8FNpbfPFSmcdi5ZmCO3Pn2FH7xjH8AmMYSeVpNuq48qB1ILd+MSjrDOrNePPKPPmgtj85j9M1gbCUbK2dz0/4ppY41m/VA5IYRPtqVyQv98O7ktXyekLEK5tvH/N4091HV6LtNxbNMqimu2KpnXuiOcUatsjd6NOzhNzF/vFo6ztpggAlghIsIhE9ZKu4YPeqy1MDann38fzr75AgDp9Jx8rJ22V4UxFaxGK8O67Emm0nVJiVjF0BmaxC0sCIIgCIGZg8W8M8q1i0vYFnVXJ2MrIUC2fxWp9EMKqyjBd30uQJgtX7z6FK68M84zKekl6OKKdVzLRGy46SLgs1ed6auqwsQ4Gh8CANYUY80R1ir92PL4m9Zsk2u3ydJ0dQ2PAUMs19EJNQWnirJLuDy+qhXrKko7UVOxSoL+3eK81an8f1elqtFKOcnH346EqpYwcXQ/ceEt56brPGRNIM/olvVZehxWt6y2i577jqOOMu7KAGagXBc5z1UQBEEQxmRRlmvXYCZN1dJyrmXYLFgAmRWbYD9aO9VTWA4feNoRXP/J7+Tfk3wpser3WzPxv173VfII7ya6z1BJcd1n5Wih2uasuixFZzvWp4y+kDHXiRB6Co7zOShFC2du4pPXPNRwprBEHhM/AgA44FW6eAMigJLG9VzjLLp4n9YARLnuIieveQgX3XLOJkuT7lPYPjZah6uruDytxzwfGHnZN1GugiAIghCSYRPwt2WRyrXvQKa2pGt1zuCVS+iFY1E6L3WPExxyjIQjJFmKiDLmPNg4W7NV1m3dbVIvmAJqhhJCMJugpomzGOXqMhbqO17aB5FnpKiwHP74yd8HAPjnn/kSIigcYoUYKnUPA/nfGAoxqdwtt0dr7FOCNz75iaPUW5gGad8RNx7XXn73cVd9nLAg5epL03irbwNhpi37w5ZJRcZbhaPRIXRjOeQVoNMiGi0oAmOP0iCWfUrSc4Sd5uQ1D+H4h84DYKRCNFddarFeqUsQUtWc19GQxP2CIAiC0AMzMI4Xr1xDuHl9ZZgJtvWbZfpdphULwFE6SGeYm9YrNlGZ2mrdz8ZYj0aH6TnCzqN4k2M4/U6F766EyLQ0brYmsVwHZ4xgprxTZHNblkUlAiKV7tiPJCBFSDMuJVgjISoGNGUvXzEpxLRxC+/RWoKZhJy1ivLVcnRCfzbcxCGYRVDTxFmcchUEQRAWzgwU+c74KUNECjfJUKD8c6hiHKoYB0mWcDtz5wjCv3nS07CHJJtes06XoMsihNMoYZVbq3uUYA8J/s2TnjZ2tYUJoJiw5ggHSdq36H5G9ztNuGana3v+YHCHz0CI5WqhbQPSypOzMZA1E+JoM7f1wRc8EKyOwrzZpyR9CeM4H1s94M2+GCrLyIR8vyA8+NwHoG66DAdJ2nUnKkqnyBBvLR7iw9DrsXZCEvcLgiAIQniY239cIKJ/RUSfJKK7iOgviOiobx0XrVy75hSOOrwJJirCoYpwsI5bRfMJy+YoHeIorXE0OsxcwGl08H7mKk4jhNfZR+a4ChsUEw7WMQ7WMQ5VhES168Z9+raqPnOJGeeI6BIA/wuAq5j5SqSZO17hK2en3cI+kcKmq9ilQUXEYE7LWGIDFLpxzEgkEWfr/R5w6v7dLynVY5JAQjCIiPO+K+JmA6A8ZcY901L6d5JRwP3XaQXgUUR0COAYgK+2EeAEEZ0H4PEAHgHwJWaWPH6CIAjC8PToCWTmrxDR/w7gy0j13d8w89/4yqlVrkR0DoBfAfCzAPYBfAPAUQCPI6JbAfwHZv4730LnTPkt0XSX5EsxEeXWSEwKR+M0MEXSAAia33jCs/H6e24BFBBlJqyey7qH1D2sLdbfeMKzR6unMD2Oxmuo/TMAgMRITGNz3eZz8KewTFxAOjoDjxPRHcb3E8x8IpedGpIvB3A5gAcB/Gci+jlm/k8+hTRZrm8D8GcAnsfMD5o7iOhHAPw8ET2Jmf/Ip9C+qHLz+oSf+0QK68asz9FyYmzGa/ejBMdWaSf5LWfJwi5wlBQQHWIvcwcfcpqUPZ2Co9L9glDirL2DvK85UPFWxiYNM3nlBdbu4jYLr6fbdbnOl9KO7lNqHmDmq2r2vxjAF5n5GwBARH8F4EcBhFOuzHxdzb4PA/iwT2GCIAiCMHG+DOBZRHQMqVv4WgB31J+yjc+Y69MAPNE8h5n/yrfAuWCzZvU2m8VqRhZHxFiRwipSOBof4mgslquwzaMJiJHgdL60nF4dh3GUFI4OtdCwMCuOrc7kS1eeTvbydIhbOYez4yMuuocVU8FK1cfPB+p7zPU2InobgDsBrAF8BMCJ+rO2cVKuRPQmAE8D8EnkMY5gAItVrlWUxzVMxbrKGvwqShXrfrTG0fhQlpkTrOyRTtSftpu8MwSwT5TvFwSTk9c8hEtvOyv/fkArRIpTJVuaXamVpu63JrV0XBd6dj0z8+sAvK6LDFfL9VnM/NQuBc2dKkvWVKyraKNctWJ9VCzTKAQ7RxADRvalBHrMPlWsR3pcGFuYN9v9StaVZ29oppJVKI6PjruaTSCmOD2ohOvs41uIaKeVqyAIgiC44mq5/ilSBXsfgDNIF9NjZp5VNnHXzEyacuRwnUs4It5YrqSwIoU9UjgSrbtVWlgsexRDMQPZcIJesjAmwh4i7JFYroId3a+sVYw1KShSQLTtBgaK46lm9LBvPmGfyOPemYHl6qpc3wTg5wF8ApuhoZ2knBw7IkYEzv8CwCpK8o+swylUEYGwR1HeUWymcqXbo6l0ZMLk0P3KKkqw4ihdFSfZBChFzAClijUi3gTKzN0dDMwmcb+rcv0yM7+j15oIgiAIggNzyCjrqlw/TUR/DuC/IXULA+h3Kg4RvRTAG5DmZHgjM/9eX2WV8c0FvDUVB4w9UtijJNs2/bcsYXj2aAXFvD3cAEKECHu006m/hRp0X7NHCdYUY43t6YA+QUuzC3JakHJ9FFKl+uPGtt6m4hBRDODfA7gOwEkAf09E72DmT4Urw/3YytV1KpL5m//HUIBEfQoVRCDEeWqb9E9MJC5hoZa4NDpnKtXydpvS9FGmRBNN3j9xnJQrM7+q74qUuBrA55j5CwBARG9BmusxmHIVBEEQhL5wTSJxAYD/CdsZmn6xn2rhEgD3Gt9PArimp7KCERm5YCNSstScIAi9krqCs2hzUgAveonunDl0ra5u4bcD+ACA9wAYIvzV5q8o3E4iugHADQCwd8HZA1SpGZU37AQqS0cmCILQF2nKwyj7fzcUK4BFRQsfY+bf7rUmRU4CuMz4filKi9VmSwSdAIBjVzx+Bu8xgiAIQme6r4ozCK6vOu8kout7rUmRvwdwBRFdTkT7AF4BIOhUIJ8B+ioL1Nxe9X/ifIuFXUSBkXD6UeD8u5pD7yGMRoKo0LcUEvZX9EUmPl41CWZqh6vl+usAfoeIzgA4xCZDUy/+WGZeE9GvAng30lDbNzHzJ/soy4YyVotwPb6wIgUIhxxhxTGOsGRoEuwc8hoKKlekiZlVBwqH0naECnRfc8gxDnUSiYp1XX3kzYYZKHzXaOHH9F0RS5nvAvCuocsVBEEQps3sA5qI6InM/KWa/QTgEmY+GbpiU4XNdRL1XDHS6yemb39rFSMCY00Kh1GMYWLAhLmhwDhkhcOsRSWZ/00RAQxEJEMKgp1DTufOr1Wcf7T1CqD4P9My0h6azF25AvgDIoqQRgt/GMA3ABwF8AMAXoh0hfbXIQ1Amjz5gsGOT0Y3SJ3AX4EqE2IrENYq6wwjIOIIhxzhjFpBlKtg45ATHELhkItuYcUASCFiaTeCnbRfAQ45XSh9zRHWKtoaey27e83E+74KdzJJ+2dCrXJl5p/Jlpr7ZwB+EcDFAE4BuBupy/Z3mfl077WcAOY4LOdvhxvWiEpLGpi39gwEocwZJDhkxoG2WLPt+dirLPogVPBIsgcAOJ3s4UCt0oXSVapkgaJyLSvR2Y2v2liA5Yos5eBrB6iLIAiCINRCvIAxV2GbOtfwOguNV6rokrn0trMAACeveWjAmgpT55AZp5lxOrM2ksztFoOhdntlR6GGS287C6cNy1VbrOVo4fKw1aKYgfUtyrWCshuYsuCldG3EbOw207GqdDwArClCpDgbB9GBKaJchQ0PM3A6m0oBbIJU9ihBAkIiClawcGp9BKfWqXI9UHHlFJzN8NX2vFe2HD8rZmC5SjiiIAiCIATGNXH/TQD+j2zuqd52gplv6K1mLWC2LyWnLc8q9NtbRNx4bOG8zEWs3wITy5viAcXzfTsUeuU0Rzil9nCQWayH2dKEe5xgnxIgOhyzesJEeehwP7dcEyOfsG0WhI872NWarYoyHjKT0xzGXF0t18sB/DYRvc7YdlUP9Zk8ttB2/dEh8YmKcLBe4ZGDPTx8cASnkxVOJ+KBFza8/p5bcErt4TTv4RTv4xTv42F1BA+rIzjF++l2lX5ef88tY1dXmBCnkxUePjiChw+O4JGDPRysV0j0uKvRH5XdwYt6yecOn4FwVa4PIp3T+jgi+m9EdE6PdRIEQRAEO7yJGG7zGQpXc4qYeQ3gfyaiXwBwM4DzeqvVQFS5ke3HbtzFphu5Cp2xKWLeHN+tusKC0Fbrw7yP0yp18R1ylhiAYiSR0VokrkkwSOeuZv+DQIxsQKH6eBPX5BGSsL8brsr1P+p/mPlPiOgTAH6lnyqFoyoBf3k6Td35QL0SrSKOFCKm9O8cBgiEQTnNezjNK5xWe7lS1WOvAACVtiFBKBMRY3+VJhhJVOQcI1LGx01cNXY7mqt5Bl2qa+L+/6f0/cNIMzYJgiAIwrAsRbnuGj4Rwya5hUoMyqzmmBT24/Qt86ybj+PB5z4QsqrCTDngGAe8QoIot1hVNnBwwEBECgeZReuaC1tYPufefBwPrZPcYkyiyHuJTBtzS+w/B2fgzgwDujSetiHomgicf/aiBHtRgv04wYoUViTuYSHlNV/4OA4RZ2txrqAQQSFCwulHIcIhr7L9MQ4R4zVf+PjY1RYmQESMVfbCvh8neT+j+50mmvqwrn2gsGFnlKsgCIIgDMXi3MI+EcBtjreh3xhNV3JEjFWUWqwrCUwRDLTFesgrJEx5IoA8aIQjRFA4zH6eMTNiCRkWMlaR2liYEdIUq9m+yDBeu+QT9o0UHjyyeAZOwMUp1zIhxiN8ZUTEW590u3SQAnA6m35zkLuFtxOT6AjibIP4mAQAaR8SUbTpU7DpX7yifwO4d8eMFJ7DCNvilasgCIKwMGagXHf2fbjJZeI7cE8Wa9V8q9To5eeE3aVstepApgNepRHEWRo77To+4DhPNCHsLmbfUe5j9IeyjytsWUy9zCSXq5tB+sPFWK4u02d0I2pK4t9nVK/inX2f2Xle9ZkvA0gTSBxynCpVRJZ81TESRPk4a5onNsKrPvNl/PGTv2/wegvToO++I1SksEQUpyxGuZrowXX31IZ+x7dFcQTFMu66q5xS+wCAA15BMSHJHEe2TjMNaErnv8ascEAr7NN6uMoKkyNNvt//y/nUg5kI8xhzFTNKEARBmBc9u4WJ6FwiehsRfZqI7iaiZ/tWcZGWa5muOYbrZNSeAyos9aSXgbr0trNw8pqHvGQJ8+d7yaPy/xPotTMr3m85zqPLDxEjZsYZyLjrrnLpbWfhQG36E7NP8cU1yndy+YQ1w0QLvwHAXzPzTxPRPoBjvgIWpVxd0xa6LJ6uF04HimO0LmWYirjwQyDCWtWtXyEsled9/Ay+lxwFAMSuU7K42FYSjvC8j5/BB552JHT1hIlzoFZYq3jrhV3TZrzU3BZivHUpY61EdDaAHwPwCwDAzAcADnzliFtYEARBmBfd3MLHiegO43NDSfqTAHwDwB8T0UeI6I1E9GjfKi7Kci3Td7Ym00Jlpk2WFNqs4QpKM6jkRAAUcOEt5+L+Zz/oXpgwa769PoY9Shdw2KMkX9TBhYTTiGKdaxhIeqypMCUuvOVcAMBaxVhzhLWKsoxMeljB8IxZptX4unCnHsy0KbjT2Q8w81U1+1cAngng15j5NiJ6A4BXA/hffQpZtHI16TLuql3BrjK2GnT2da3STnJNEVYyJWdnuPz2YzijkoKfKAI7d3x6us4hxzijVrj89tQt/MWrT/VRXWFCnE7ScXbddyim3DUMYGsMVtM0Xtp2buvo460ZPY+5ngRwkplvy76/Daly9UJ6eEEQBGFe9BgtzMz3AbiXiJ6SbboWwKd8q7gIy9V0TXQNatKu4broYDPQqewaNrdr9P8RMcDp3wMV4+j7LsLp59/ndpHCbDkSpfNTN4FtEcAKCnFtBHqhDWXWawxGHMl8113g6Psuwqn1xv4xrVVzW5U16mKlmrKqXLy7FMxk8GsA3pxFCn8BwKt8BSxCudZhjqMWlJyBdoG4TMspK1OtoE3XsOlGrooc1o1RJlcsm2d9bI0zCtlY6QaFdOUbFzebXkQ9TXGX5GO3z/rYGrc+ffE/4Z3lewebqHDdz1Ql6df9SdmdyxaF3ISrS3jU8daey2bmjwKoG5dtRH6ZgiAIwqyYQ4amxShXV3ewzzmmpVt3rLZOTQvYdJMkBfdear2ukxiJIqz++vsBAGe99PNedRfmQUwqtzRVKcSh/L1ZFiPCRp7zfFlhdjz019+P9UGEOFugdRUnWT9U019heyjKhqs129bVO4iLeAbKdZEBTeWVHkKFl5uNsuzeVWUFmn3K3xUTEhUhUYRDFWOdRFgni3wMO891dz1UXLUkHTktHJMulm7/lImgtlZduu4uyfS1RNZJhEOVvoCnn2gTKYyavqahb6pSfD59JPN2nMvQELf/DIX06oIgCIIQmJ1UrlXukLZzu8oWbdmKtcmJI8Z+nGAVK6xiBXXTZc71F+ZPlXVqO6bpOGFZqJsuwypW2I8TxBHnrmHAiBguW7ClfqdxHmvgPnBwek7cH4JFKddyg3JxDW+Fr0O7VJqPr2vMdY0xIkYcKezFCY7EaxyJ11hFMn62JK7/5HeQoDjh34b5Qlb+NJ4DQgLC9Z/8Th+XIIzEKlI4Eq+xFyeII4U4Us5TtspDYqaL2J5feNPXmYrV1X3MNWX3RhfFKsq1f1zexAoNr+pNr6Fh2zpL0so1Kv54zr75gpZXI0yNQ44LK96kyjAdddWW6GbeYrT1yc/Jjkn0+VsvcdHWNB9hvpx98wWbl2+jfzCDKav6FdcX/I3CdavT1KxW6vgZip1VroIgCILQF4uZiqMpT5kxv+s3tXJy/vLSclX5hquWoitnc6qatmMeE4FBWdSn/g5IQv+5o6N3dZL9QxWnuYENixXQ1kf1u21xX+rZSJgAipDotqmQvx5fd9dDuPHKswJfjTAUOkH/6QR5hLkiRj7cavRZNkuyzmp1WVouhEt4UGYwFWdxytUHV6VYtVpOIUMTuy+mbirWXLkaGZ2EefKiux7OXbQFxZoFJRXT1qVaMalxVMVg49jNmHykz6GooGBfdNfDAIC/vdJ7dSxhZMzUhuZ0q3ylLXZbGL0qc1MZ35iSqSFJJARBEAQhNDNQrosZc60LLNqOInaQ55xfs+rN0L7ddDdbl6/jCBfdck5zBYVJ8YJPPIIzag+njc8hx1hn7mHtBtafJIv0rYsW1sekx+lzN+u6atm6vDPZ5wWfeGTs2yF4cNEt5xQC2TTl/qEq97lLH1TlFStP6aminDSirn8dxOKdQbTwzluuda7hYjpDhwUAStvrUibazlVMWHOEyFCw9z1bplnMgVPJfu4CBpCPsZoRwaYL2OYiNomolNRf/2u4iNO2m66UY7bjQ0j08Fy46JZz0kXQLRHiwKafqJrr7DSeWmFcVCX5r5MluLPzylUQBEGYEQOnMWzLopRrnRW6HUVcbYm6BjaZ0cPleWhNwU3lMpXhIlyreEEO++Vz5Z0xHkq2f0pm0JJtfmr+v80tt2XNqvxYHdCUWqyAIiOqKa8T465nJp5XIgzNmiOsVVwIaLLlCW6icgk6xyGxOqu1Lo/waNbuDJTrKF04Ef0BEX2aiD5ORP+ViM419r2GiD5HRP9ARC8JWW7bDCLFMPVq2b6ygO0J4WuVuogO1AoHaoXzP/hY7/oKw/GE2x+Nh9dHcEattj6HHOEwGyc1x1vXHG8lXAeKKRELiz2A0nOMcdd07DWVbyv74fURPOF2iRqeMud/8LE4UKtMwaafugxdbmOjzZmYXGXZZE8FSdxfzY0ArmTmpwH4DIDXAAARPRXAKwD8EICXAvgPRCQDSIIgCMIGCWiyw8x/Y3y9FcBPZ/+/HMBbmPkMgC8S0ecAXA3gFlfZrnNX0312N6/tvKbEEvocIE0uYXMZl+tWfjtdI0IEhko2x60pwlkfuBAPPe9+11sgDMDxD50HAHh4DaxILwWnEDu8Grsk7LeSD2PY3zc3aRWjPDju+If28cCPfruxTsKwnPWBC3E6MbxWJVdwpfXKOnXhtjvWOf1hTcKI4qyL+muQAKh6pjDm+osA/jL7/xKkylZzMtu2BRHdAOAGAFhdUJy60jSGWs7aVKdg07KKrayc6alqjLUyoxMIEW+U8NZ5piLP/kR/dwnUC79iuxXCwBz/0Hk4UOlPJwIDUZZVhwkqTwqyvRBDXUYmFxTb5WrZm0jlzdi9AuUvAqJkxyf6u7Q7O53YX7LLipWNZ6rxjRB2XrSk5pzRp96U2OmAJiJ6D4CLLLtey8xvz455LYA1gDfr0yzHW28jM58AcAIAHvUDj5/BrRYEQRA6M7B7ty29KVdmfnHdfiJ6JYCXAY/5DjwAACAASURBVLiWOX9POgnAXNj0UgBfdSvPnqLQ9dy0TvXHlV3D5Yhjl/mt+jjTetXbzb/lcxKOkLz7STj2ki/4XJoQmEtvOyvP/6rJ0xMS8h99leu2K3Vyi67F1Eo263npbWfh5DUP9VIvoZlT734S4nXqeYipehk502qtWg+6jGt0sG/CiDZ0Pd+tkAHK6MgobmEieimA3wbwfGY+Zex6B4A/J6LXA3g8gCsA3N62nDqF56IAbcfWKVgTfY6tDmUFq0lKU4d0eYoJ6yTG6fUK33v7PwIAPO7ld/veDqEjV9xxBI8kKOSENlFMw65pZSvfwMxbrbnijiP47FVnBq2XAHz97f8I0QHj6GoNAFjFSdqOsjzjmoLL16JYbUrUJX9wXcKIsrxtOdNyCc+FscZc/y8ARwDcSKlmupWZ/wUzf5KI3grgU0jdxb/CzM6T9Zqs1/LYZ9PcV5fMTeZ5VWO95nabgi1WcrvxJipCoigLp0/3f+ddV+Cc6z/rdmOEzjz1wyucUZvnrsc+N2uvRttZlUbCHNuNSJWsbMJTP7zCp35kPUbVdpLvvOsK8AGBmZEo3QdFQKS9HdVtxqZEqxRr1ZipTyYm1zmtNgaxWJGtyyqWqx1m/oGafb8L4HcHrI4gCIIwJ0S5Dk9d9K/9+PaZm1yn55Sn4mxcvkUi2o4e1pPCiYBVrIA4PWsvTqBuugzRtfdWXpvQnSvvTMc40wn+VTmAN9brlKiKWFYc4co7Y8ngNADqpsuwt06Q7GmLNd2uf9dNfUCVG7bNtJvyeT5Wq1W2x7Sd0NDQBbZgccrVhs/c13S///xX1yk9tu3mPrOO5r4I6Q9Plx1HaUDE0fddhNPPv8/1VggePP0jEc4ken3WKBv/bqdAQ7qLXdcN3q5DBEWMQ07XgX36R1I5H3uGfXqP0I2j77sIpxNGHCnsx8nWbxqw9CcVyrNun4s7uEkZz2pO60yihaf1qi0IgiAIC2Axlms5MQRQ7x4uJ/VPj69+HapLLmFar2a5VWWW3c2mzCooi1DVZWp3tGLC2TdfgP0oDVCRRAFhuPLOGGeSOLXyAKw5zpeQK0ZxbpYKc6VNXleg+MwbjyXeBFnp8lhnccqm82SLDVx5ZyIu4kCYCUZOrzcBcHGkUjew5dnV/farEvLb9re1Wpvq02S11iWc6AsJaBoZl+hh18xNgPsUHVPJ1p5TMXWijN6uFaz+X5e7VhGirMM8/qHzRMF25Afv2McjSZYkP+ukEt5kPDLnkQKlDDseitNVIZsrJ9UeV4oK1go2H90nZCsuJQDi/BoOOcIP3hHj01cdONdd2Ob4h87DWsVYq+yFC5uX8AiM7bWL7FS1i66RwTbZPpmY7HVqPKQfRLkKgiAIQljEch0YW6BSXfSv7Zwu7uWm+a+67Co3tO2N1SV4RSdpzyoh1msLLr/9WP7/I2pjrZpu3zqr1baUYCic3cBGHSJwyXpNa6mtV0WbubpxtqTd5ben3cEXrz5VFi80kFutRo5nF3xdtHWWp09kcPnc8vF1AVUu5/eOKNfhaVKwgF1BAm7JJfT5gNsUHZcoYrNsn+vUCb7XiHLPn6JUCZx783FExPjWc77pJXcXefytZ+ORxNKZWDLkqFLnaaYctNF3ZGVUMcSQV73QrFIXcd4+LSkSgfR+fPVZ3+2lvktCr7WsmHCgKF+PdW2+kHG7NaSblFtTWsM+FavNFTyltV6nwuKUqyAIgrBgWNzCg8EVb22u0b9WmQ3Wa1lGk/Wa1mdzHmCPFra5iusSfK9VhFWkUusVQMTFubTH3v84AMCpH/t65bXuKufefBwAcGpdfBbl1IYamwu4vFxYFX3Nc60aSsi3m7sZUBxn+83jy6kcKb83APDgcx8IVvelcOz9j8NpI8A6XyoOm4hgbclW4RI57GJ1+kYGd0kWYWMUq1WU63D4JoYAmpNLbLuLmxWsebwtyX+6f3OupmkN2TI66b9C+iMuJGnn7Uwv6qbLoJiwevGXa+UunbNvvgBAes9OG4vSF+9ftJn2YonQLXc2rgrW5xgblS7giv1bGb/KL3f51J54a7Fus67H3v84rCKF7z73G63qvSTW7/k+RMQ4nWwSb5j32bx/2iXcFOXtMr6ZHrd9bmjF6utOttG3siXMw3KVJBKCIAiCYEBEMRF9hIje2VbGYixXGzbXLuAe/Wvb1mb+q6YqyKlKRlrXjZvXzEVqWrqqJMf8X8tZc4R1EuPgnU8BADz2Zf+AXUJP7j/I0hnqoJOIty1X62pFGba3eddtVnk1Vk050KhOpi0vta1OBZeyJcJ5y3IFZe7k3U1W8s3sNwMA+4cJVnGCRKlNQpeKqP+yFWizCttYrbY202U+q22bq9U6WiBT/xNsfx3A3QDObitgUcq1STFutvlF/5a3mc+1LhtTut9tHNZWD7NcU8GmspC7+SLWmXeKFMZ/khgHSYwzhyswA9/6/54MAHjM/hnsXXfP1rlL4fG3no0DFeNArQpjYFtjkmbnZXkOhXNKVG7vuLir6/l62o2TTMsYfvmFzBolnSUr0V3G8Q+dh/0oWXRU8eGNT8D3Do7g9Jn0mgsvxLH5slwRE2G42n2Wi9vsr5Zpk2GT46oQQyjWIRVtn25hIroUwE8iXZ3tN9rKWZRyBdorWKB5mkyIebS2ubCbcrbrYcoo1y+3WC1tOh/vKf2o40jlGWOA1ILb+9tLcTReYykLAFx621lpJiIABypK53XydupCYFvJKiavNIO5nBpF2Nd0nLrMTXUKt+llwGqJo2jV6vt70S3nYBUlOHnNQ/4XMFHShPsrrNfp2PveKo1cMhVkoiKwvv9G5jSTsnXqqljbKFVXOb5zWW1ympRz7zC6BjQdJ6I7jO8nmPmE8f3fAfgtAI/pUsjilKsgCIKwbCyrKfrwADNfZZVL9DIA9zPzh4noBV0KWaRytblk21ivtnPbupk1VZHEZTlVrmpzu5abbI0tb7ujdNmrOLVcy/Vdc5RHhQLAitSsElA8/tZ0aEQx4XSyWXtVZ8xxnRZRN97ahMuYZ1eapuHoY/q4Bn0f00IApQiKFNYc4aJbzinUb07u4vM/+Ng0JiG7tgOjreTrKGO7H8jjH1D8DbZdLi7dXzqvhcVqk+Piwg3hUl4AzwHwU0R0PYCjAM4mov/EzD/nK2gxypWx7R1tmlqTbmunYDU+gU5lOVXBTqncjQygqEzTczds7av4IUfZ7Gui7U46TzauFXOUBrDoek1truPxD51X+H6gADN7kr63WqmW5yAC9qlQXRWg7/lNQUptjg99DXmwFCF3x+XTvyIgYsI6WzgiytIqlp/PlAKgzDm8CpvsSjrIzXwRKwzPbL0MG/EP2J56UeUG3t63XUeXoKWyHJssV2UYUrEOMkump0KY+TUAXgMAmeX6m20UK7Ag5SoIgiDsBnOY57oo5epivdq2uUzRcZFVFZykqUsUURXo5BKNbHcLbQdRaFnmurD62Nx61fn/swQLusyj77soP950tUbX3rtVdmh02eZ0mdNJdTSvLRFCndVq+96VEPJ8ZLhMxelSj4L1CuSWw0ESF9pJmsBke2GDsz5wYeH5DBE8p266DACwihQiYqyygboDIz2V6dUwvTfm72uz7EGx/VQN2ZiM5QauOrdNBihXWYPpO8YQU3HAzO8F8N625y9KuQLtFWy6zX1qjK8sLa9qgXSXrE42d13deIdLZ8tMUESAioBI5R2MLdPTJlI0QsIREhUh+e8/kNd9L06wFyXYjxPsR0lh/Fan2NvKHMWbZOemfP09HUPd3A/dkZv3qHZqSUXmoab7U5ca0yuKuOfxKB9XcNtrsipu2kRW215EyzLL0cZa8ennuSKVK0FAtxkutBtTpuIoHyc90JHhSYxDFeMwicFMiNfZqj+RQlySX74us26JStfxbRpLrBqyqbsPmjZK1VWW7fyQSno0V7DBHCxXydAkCIIgCIFZnOUKbN6izPerrhYnUB/9WxXklH6vlmWri0ugkymrK7m7T0VQNVlnzGxPiYqwTlILtnAPVpvjVyqzXCNVcB1uuWy1W7oUUKLzsup7QNpqLdeR66No69zBtfck4HEmbaMrbfmmm8p3iS5uohBJbVh6ZlRyVTYp27M28+1GSIco1hRhxSp33a4psj9rQ472cGhvy6GKcbBe4VBFqeWarcu7igkqovw3VmVlbwKUtr0cLveniT4Dl6rOc7VYbTInFcA0iUL9WIZyrbjRZRexzxSd9LiivDbjsKY8U2aVy7lcv7qkE2WZPpiRjlH2N40GpWx/qbNE0W2mXWeFYzidXJ9/p01nVe58y3LNiGXd+ZYjfyMwmDmftG/r9K1leHQwbRliOoJPGWR5kelC13tsKtWyuzh/rpwqQV1e2f1fkI3N9KrDzC2cKPtau4kxfUiXZZO3VVdutxarpm5YMOTYatX5bZW0j7wxdBxhHm7hZShXQRAEYTdgHiSgqSuLUa42CxLwm/8KbM+BTbdtzq2yXrfP3d6mZVa5dpuSTpjYgp6aqIo8TgNTNvP1bOWZqRSZybg32+5AyixYZVhPpnvP5qa1ueXK91C7r3Uu5bJ15jsvVF9LH3TNK+zLdkKSMOXXWcAugUzFudnFNrTZjvy5FgLqLH/LAUi2FJ9mvdOy0vJ0ykLVYPaYMm3X3qZfr2oPXdy2tvNDB0LVba+6DQtLKtGaxShXoJuCrdvuohD1uYBdyZrb69zOVTLLcqt+rFVutPJ15GXl8gA0dJbl6T0R2BLdnHZ4a0QwE5tXdWhV41xs+T+f7pDtikoiE+O4rj/w0MoxpBva+ow61req3TTdR/2SU9mpWxYBsCkEheLwgXbd1rWb8ksY17TLwgtkxcuYrf35KtY2a7fayqkrz0cB9iF3CkpV3MKCIAiCEBpRrsPja70CNtdt9yCnKrnmdlugk5apqbJibfKBZovWZoVvLWVXouot3pRbqEPZzUzYWAvYzntbZdFURWnXrQaUHtj/G3Tf81dDll03vzWX2cXyralPpYu41MbLz1Yh+10Z7aZcz7Jb2KSuXebfy5dR0f42+81z3e9Xm7noXS3WKtmuAUyVZdlPH9wVLJbrSFQqt+yvi5KtG4d1cefW1sNa3ma/i6It19FW10IdCxG3xbIKZdjO9XCNmYpVyyqMj9ZMPajqaJo6Yhuua7GGZCpjTbZx8NC4TvFpUljlbVXKryq6tjBWX1uP5rbuoozaun3ryrCV0yRvaPcyMB3FCgagpq9dl6NcLaZpnZK1NQfXcdguFqdPIFWd/HIZNplluUD91B5bFqjyMWVZ5cAqLaMcINU0nlt1HVWKoi4rTlUqurpy+qQvxV5lkYa8psqXNUdLp+r48nnWwCmHNlO2XG3WZdU0tuogI9fj3O5zk9L3kR1CqbaRbxNTef3T13uDsBzlKgiCIOwGM1Dgy1KuFSapdQw1++szDmvfbk/iEMJVrOVrqtzRJk2u47qpPW0iIZusV5OqdWfrqHJ91117H9bomGOsVXSpk8s4LNDuXjZZqT77mtpMk9Va/t91Cltba7V56lft7mAWa1VZvU61GVDhyZjrGHgo2JrDnV3E6bb0b52SBZqn7ZjbNVXu4s3+6vLKZdbJ7xLQUu68dD1tMl3GuaqoU54+ruDQTGWstUyVSxfo/2WhTRCPjaq1ijfllGQ3tGPfdt71Otoq07p9Q7iZp6xY0/Kmr12Xp1wFQRCERSOW61hU+HzbRhGXz/ENRNL4JIlw379dTlNy/7rpPW2oyoXsmqjclTqrfKrWIzB8QJNmCvek67W7RKq7Jr/vSleXb3pMu6GQpvJDWazADKzWmbBM5aqpcRED7kpWn1On/Nq4iutctk2KtuoYs+zisdvll3Edh3OpV1W6xhDeHNtC9C44zfecgEJype8pNqHLDOnJ61OhhnD3Fo9tUsrtXc++SrVuX6upNmMoVh6pXE+WrVwFQRCERUEASMZch6PKsqwzR9sEOQF+7lufaF/fbEx18zy3jy0fs31eW6vEJWiqizt4Oyl9/fFVS/D1ZZVO9XfeZinCtvfI9x6EyN0cKmiquRzX49zLnIq1mu/3OadDkFYw6jLITITFKFegXvk1uYjL55kSXMdim+rRZkzWVmZd2a7H+naGdR116HHdqsQXVfgq3yngq1zqFmSoYuz7EDIyt7Ycz/NC3JdQilTTZTy3rVIN6Qaewvj+1FiUchUEQRCWj7iFR6LSgq0zR1HvJq44xclVXOemdXXPulqzTXVxPX77/MZDMlnp37o38a5Rrr6W7RJYwjV2sW7c01m2LsIiqz/Lu2vQVJcgqbpb1MYVPEpu4enr1mUqV03lOCzgHUmsT4H9tFpXcf2+ohyf8ULf3LI+buQmfCKV0+M3/3eN8vWpr8uLQx8M3eHM7Tq7jon6KtGQz8NXls+19pF0AmjWRW0jgsdxB/P4Yx4OLFq5CoIgCMtDkkgMSFWO38YgJ8Ar0Mk8reLUVlHFm/3b8nyjX9uuluJj/fi6mX1eNF1yKJcZYoWYKTO162xjmXY1RkLdgykFVnVJOgG0dAE3nNjskq7dvTOMqlyJ6DcB/AGAC5j5ASIiAG8AcD2AUwB+gZnvdJVXpWDTfQ0uYsCuKRvOdXUVA355hIvHbW/zjd4F+sno4zu+23Ruen7Tedvb5pT8YVdp2+m2VXZ9vXD0/eLQZTw1P6ZLGa0jkxsKDckMNPhoypWILgNwHYAvG5t/AsAV2ecaAP939lcQBEEQAAZI5rnW8ocAfgvA241tLwfwZ8zMAG4lonOJ6GJm/pqr0Lq5pE0uWRc3cZMFW3G6mwyvuaz2MtpYtDZCLknWJZAqtItZCMuYrtyQlmkoz0efgVZtEkF4yZiDxTpqoX6MolyJ6KcAfIWZP0bF3u8SAPca309m27aUKxHdAOAGAFgdP8daTh9u4kYFXX+6swzbsW7H1+52VjZtOpsu456hFfDm/MbThZ4YS3H2OUTQpT31NbWnLxewy/ny+6qmN+VKRO8BcJFl12sB/A6AH7edZtlmfXzMfALACQA48v2XVD7iuvmkTZaky7zYfHfN1J0aEV7WatXxTecUz3c6zJDrfuwQgUddxnib5AjVjB0spAmtOPtQDm2u0cuCDSGv05Qfhwr0zRTq0EBvypWZX2zbTkQ/DOByANpqvRTAnUR0NVJL9TLj8EsBfLWvOgqCIAjzQzI0WWDmTwC4UH8noi8BuCqLFn4HgF8lorcgDWT6jut4a/N4Zr2LeFOf9hHFTTJcxmXLshrrVXOOrwy7XO9Tgoz79uFi9pUjuDMFi3So/jZE2/F2EYeS2cFiTfd3LD8Uoly9eRfSaTifQzoV51W+AmqnzdQEO5nnp/vbuYpNGXVybFt9FW5TGT4y6vBVyr6pEm00dbZ1AVehfuQhXMd9dzhzqCMwnXVgi3L7ve7WU4hClxFgPm2XzFHBYciqOC4w8xON/xnAr4xXG0EQBEHozujKtQ+cLMeaQCdXGY1hwdh+o6sNVrJsc3kf7Msl7CO/TXkuFkmItVldpxWVmYPreMw6TsWVO8Y9CBbg1We5PbqAW9UnEASWMdfBaBgHaJqq0RQR6+UqBpyVrZMCqtnn06y7/gj8XcPu5YWYShNqfq8QhikrzaEVQpdbEVKZ+sicqmI1Cu9NdJbg6M+QznZRAE4w8xt85SxDuQqCIAi7Q7+W6xrAv2bmO4noMQA+TEQ3MvOnfIQsR7kyVS6V4BJJrKmzgBot2PxA438Pl7GT7IpibIR8r+zylhoi0rlJTggLV6gnVH82lTmzoQjZzXtfU0BrNT0ugJy+n0vPAU3ZDJWvZf9/j4juRprMaEeVK7B5qA7TOEJEFDfJSg8ufXd0Qdvwds16HV1N15/KUNNnQi76LrRnjKkqrcrovQRLmV2vy6PSody/TrIm8qITGiJ6IoBnALjN99xlKVdBEARh8XQMaDpORHcY309kGf+KZRCdBeC/APiXzPxd30KWqVwbLFjAzcXr6y7eHBvOkq0rx0Yfaf26SgwV8Qy0Sw3pK0fYZoqu3Ck+waAWt+cF+i0CEFDeGFZrN+X6ADNfVXcAEe0hVaxvZua/alPIMpWrxkPJAm6K1mXMzisiuGp3y/Y6pQQKGldJIZRw3wk1fBhKgU9l7BHoXpexFOa4ka8dTvVeDCCgzNHuGfcdLUwA/gjA3cz8+rZyonBVEgRBEITZ8xwAPw/gRUT00exzva+QZVuuGgcLFvB3FafHNhTdNiLY5cWsxxfHkG/yztaki6wmGRNyDU/JouzCVF25s7m/I0ZZ+xh4k4gEdoHRq+XKzDcjQO+6GOXq5LY1G0bH8dhyuY1lW+SbtOrou7SvIVOBBhhPzWU1yXAqKZM1QRd6SEbJXNT1/Cl03mVGeLyt8xL3NU4bMPI4CJJbeHiclZ3neGwqr8nyLX73mVMZcgqOW4HhRW7hO2UvQMCWz2WF6sYnqRACMmlrc5rvNU6ECxLrqdxAc2P7QNIfCoIgCEJoRLmOS0hXcSqvmyXbWBfHcuuYlIuyBzMyVLQw0M7ombuN2nfrCGaZTqgZ+9D/MnZtzvF1IU3XYp0Ti1auGpfk/OmB7oo2lesx5caoi41QafnGdFF2UuyBArhCKl+r/E5nz5u5KM4luOm7Kq6+FGoq27MyoWEAauxKNLMTylUQBEFYCv3Ocw3FcpQroyFJfvF7o6XYwYpN5XvmAW44fA4J5/ta93VTQJMwBxGedZyUmz0wvVl4AW/ZEqxQG6F1Q6v71Je1OsRPRpTrwOj77dTJZof24C5O5dsFt+2s+2pLYyvtoCkdQ87T0SIX2rm3ZsQ5m2MwtT68e+L/GSvUQnkTezAWJEOTIAiCIARmWZarxnypacygZBzqa8UCzpbsprywFm1Xhn4B7CNaOphlmwt0PG5pTMFV2SgzuMhJEi54rN+sTuOs2wcJaJoELVzFgIcSsDXeFopyrBVvhqaPRc2DJ+BYwhydgZrKWIkQ5kqv7vAOsievUMsV4OmnaFq+chUEQRCWxQzexnZHuXq4ioFuqQy7uo4rxc4gCCSUdR0qerrvea+bgsKImRpzs0zn8BsJQtel/do8j4W28b5YiHIlMJP/2JuX69Eozbdd1/0QFuDqNel9Ok5eTp08HznhFhWYE6GVUB/Kc2cUZRUhVyPq8nz6yl3cFhlzFQRBEIQeELfwsPgsFZeeUPru7Grc3tZ6vqjPW97CrCdNaNdtKKu2KLPb23hIy3csiy5kf7azVulA1z2kpZqWN/DzFOU6Dq2zJXmOyxbL3N4WPEHDnDukDsolaLL+kTJhzUWZzEaBzuR+DkGQZzYHhbopWZTrVPC2aIHWVm2xXPv2sbMijUKPFnrIqTh9TBUamz77oaAdrCjMWoI+x5ay5vKSOAV2QrkKgiAIC4EBKJnnOinaLBG3Odn4v2t6T4+i52QhBSNgdHXwBBO53Nanzg6xTvtjsHY09BJ2fTODH+BylKvn9JpgijYX4ifCuajptyFngrwoVP3IA2fFWtrUmyaCd549dsZL+k30xhjKdMjnMoNGsBzlKgiCIOwALPNcR6GF+7brWqxb5ZaZmEdlLHoNFgqcqGMpySXmFLE7A2NkOkwpaYc8NyvLU64mLTIxAYGUra0edYgCzgk+TzWgK3lL9NTGovpgKpmCdokB7lP3tWHD1KNNuSyJ+wVBEAQhMOIWnggd56wGt2SthYQXORlCBpuGnDu8QzmfG5mji3fHHlEXgnhYpnS/Z+AC2Q3lWqbjtJqpLXg+eQZYpDx4wg7Xzmjqz7xnt3XwPm7it3PqhJ02FU7ULrKbylUQBEGYJ8ySRGIwuGWKw+zcAp2sp4HWDl0qTbdniikou1oKdW1iQsFSU0i9J6T0GwHe4dQh26u4hYenU3IIIGgmpi3RC5neMRpVtyWg0gUGzoq1VAUKiBL1YLSo8wDPaIy6s1iugiAIghASWRVndIJasbmg9vVxLnZC1kyfBLPQAz+nyVi1PdBLn9RzP7crv4dBmKmlOkcWrVxNOivaXFDFdmlv3vQ6Rt3Tc/JVTn0r40Ff4HsqSzrrHphSBqfQMGSe61TpZd6qpD8MTi9TngZ+OZqB92obUaLTZ9e9BZKhSRAEQRDCwQB4BpZrNFbBRPRrRPQPRPRJIvq3xvbXENHnsn0vcRbYMYTc/ASHO3ymVJex6ly+hNLzCvLcRrqWUenpmnt5Pk4FD/CZSp37qErI31KfMKeWa9uPA0T00kwHfY6IXt2mmqNYrkT0QgAvB/A0Zj5DRBdm258K4BUAfgjA4wG8h4iezMyJk+DyQ227yIOlcY02PWaOHbxvnQP1u708N5fTp+pBG+hFpz/h/YluzRTr1ALJ5FQNEcUA/j2A6wCcBPD3RPQOZv6Uj5yxLNdfBvB7zHwGAJj5/mz7ywG8hZnPMPMXAXwOwNUj1VEQBEGYIKy49ceBqwF8jpm/wMwHAN6CVDd5MZZyfTKA5xHRbUT0PiL6J9n2SwDcaxx3MtvWjoAulNFcXbtAj+6vQZ7bEO7IkV2Hvd3HEYcXdoFe2/+Yz6lft3AQPdSbW5iI3gPgIsuu12blngfgWQD+CYC3EtGTYHewWR8fEd0A4Ibs65l7XvmauzpXeliOA3hg7Ep4MLf6AlLnIZhbfQGp8xA8pS/B38O33/0eftvxDiKOEtEdxvcTzHzC+O6sh+roTbky84ur9hHRLwP4K2ZmALcTkULaeE4CuMw49FIAX62QfwLAiUzeHcx8Vai6D8Hc6jy3+gJS5yGYW30BqfMQlJRXUJj5pX3JznDWQ3WM5Rb+fwG8CACI6MkA9pG+lb0DwCuI6AgRXQ7gCgC3j1RHQRAEYff4ewBXENHlRLSPNMj2Hb5Cxprn+iYAbyKiuwAcAHhlZsV+kojeCuBTANYAfsU5UlgQBEEQOsLMayL6VQDvBhADeBMzf9JXzijKNYvA+rmKfb8L4Hc9RZ5oPmRyzK3Oc6svIHUeTs80IgAABuZJREFUgrnVF5A6D8Hc6luAmd8F4F1dZBDPMj+bIAiCIEyX0TI0CYIgCMJSmbVyJaJ/TES3EtFHiegOIro6205E9H9mqas+TkTPHLuumuBpHweCiH6TiJiIjmffp3yP/4CIPp3V678S0bnGvkne5xDp1vqGiC4jor8joruz9vvr2fbziehGIvps9ve8setqQkQxEX2EiN6Zfb88m2P/WSL6yyxoZTIQ0blE9LasDd9NRM+ewT3+V1mbuIuI/oKIjk79PvcOM8/2A+BvAPxE9v/1AN5r/P/fkc5XehaA28aua1avFwJ4D4Aj2fcLs79PBfAxAEcAXA7g8wDisetr1PsypIP79wA4PuV7nNXtxwGssv9/H8DvT/k+Iw2a+DyAJyGNnP8YgKeOXS9LPS8G8Mzs/8cA+Ex2T/8tgFdn21+t7/dUPgB+A8CfA3hn9v2tAF6R/f8fAfzy2HUs1fdPAfzz7P99AOdO+R4jTbDwRQCPMu7vL0z9Pvf9mbXlinRi79nZ/+dgMxfp5QD+jFNuBXAuEV08RgVLzDXt4x8C+C0UJ1JP9R6D+f9v725Cpa7COI5/f3gr34qwMgoXFxdtKkhLsJQybVEm2iIikBRctamsTZi0KCSKwmwXoRSlbXxBhaAEM6IoS03TsBcl8aVbukgTykB9WpwzOF3nKtx7ZuZc+H1guDP//7l/nvtwZ57/2zwntkTE2fzya9L31KDePBdpt9ZuEdEXEbvy89PAftIH6zxSQSD/fKQ7EV5M0gTgYWBlfi3S1wDX5SG1xXsNcC+wCtLNnxFxkopznPUAoyT1AKOBPirOcycM9+K6GHhd0hHgDWBJXl62jWI5nWn7WJCkucCxiNjTb1W1MfeziHSEDfXGXGtcA5LUC0wCtgM3RkQfpAIMjO9eZBdZQdoxbPS9uw442bTzVVuuJwIngHfzqeyVksZQcY4j4hjp8/cwqaieAnZSd57brvr5XC/TRnEW8GxErJf0GGlv7wEKta8ajHa3fWyHy8T8Auk060W/1mJZFTFHxKY8Zinp+9JrGr/WYnwNt8vXGldLksYC64HFEfFXOhisj6Q5wPGI2ClpRmNxi6E15boHmAw8FRHbJb1FOg1crXz9dx7pUstJYC3wUIuhNeW57aovrnHpNorvA8/kl2vJp34o1L5qMC4T75DbPrbDQDFLup30htmTP0AnALvyjWNVxtwgaSEwB5iV8w1djvkSao3rIpKuIBXWNRGxIS/+Q9JNEdGXLw0cH3gLHTUNmCtpNjCSdAlpBekSRk8+qqot10eBoxGxPb9eRyquteYY0gHNrxFxAkDSBuAe6s5z2w3308K/Affl5zOBX/LzzcCCfEfrVOBU45RKlw2rto8RsTcixkdEb0T0kt74kyPid+rNMZIeBJ4H5kbE302rqswzhdqttVu+XrkK2B8Ry5tWbQYW5ucLgU2djq2ViFgSERPy/+7jwKcRMR/YBjyah1UTL0B+bx2R1Gh8P4vUsa7KHGeHgamSRuf/kUbM1ea5I7p9R9VQHsB00rn9PaRrP3fm5SJNdnsQ2Avc1e1Yc1xXAquBfcAuYGbTuqU53p/Id0DX9gAOceFu4SpznGM7QLqGuTs/3q49z6S7r3/OsS3tdjwDxDiddGrv+6bcziZdx9xK2rndCozrdqwtYp/BhbuFJ5J2qg6Qznhd1e34+sV6B7Aj53kj6VJS1TkGXgJ+zJ9tH5DuyK86z+1+uEOTmZlZYcP9tLCZmVl1XFzNzMwKc3E1MzMrzMXVzMysMBdXMzOzwlxczczMCnNxNeswSaNyb+kRBbZ1g6SPS8RlZuW4uJp13iJSG8xzQ91QpJZzfZKmDT0sMyvFxdWsEElTlCZoHylpTJ48+rYWQ+eTW8FJmiHps6bJsdfkFnJIOiTpFUlfSdohabKkTyQdlPRk0/Y25m2aWSWqb9xvNlxExLeSNgPLgFHA6ojY1zwm9w6eGBGHmhZPAm4l9cr+ktRw/ou87khE3C3pTeC9vG4k8ANpAmpIrfKWteNvMrPBcXE1K+tlUiP+M8DTLdZfT5qWq9k3EXEUQNJuoJcLxbXRwH8vMDbSJOWnJZ2RdG2kibSPAzcX/SvMbEh8WtisrHHAWOBq0hFmf/+0WP5v0/Nz/H+nt7HufL9x55vGjczbNbNKuLialfUO8CJpgvbX+q+MiD+BEZJaFd7BuoU0G4mZVcLF1awQSQuAsxHxIfAqMEXSzBZDt5CmbyvlfuCjgtszsyHylHNmHSZpEvBcRDxRaHufA/PyUbGZVcBHrmYdFhHfAdtKNZEAlruwmtXFR65mZmaF+cjVzMysMBdXMzOzwlxczczMCnNxNTMzK8zF1czMrLD/AOPRvsLgnxBdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "aspect = (xy_dimer3.xlim[1] - xy_dimer3.xlim[0] )/ (xy_dimer3.ylim[1] - xy_dimer3.ylim[0])\n",
    "size = 6.0\n",
    "fig, ax = plt.subplots(1,1, figsize=(aspect*size, size));\n",
    "\n",
    "x = xy_dimer3.grid_x.reshape(xy_dimer3.ny, xy_dimer3.nx)\n",
    "y = xy_dimer3.grid_y.reshape(xy_dimer3.ny, xy_dimer3.nx)\n",
    "#z = xy_dimer3.Ecartint.reshape(xy_dimer3.ny, xy_dimer3.nx)\n",
    "#Try log scale\n",
    "z = np.array([np.log10(i) if i > 0.0 \n",
    "                          else i \n",
    "                          for i in xy_dimer3.Ecartint]).reshape(xy_dimer3.ny, \n",
    "                                                                 xy_dimer3.nx)\n",
    "pcm = ax.pcolormesh(x, y, z)\n",
    "#ax.axis('square')\n",
    "#plt.axis(np.r_[xy_dimer3.xlim, xy_dimer3.ylim])\n",
    "plt.xlabel('x (nm)')\n",
    "plt.ylabel('z (nm)')\n",
    "plt.colorbar(pcm)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "_draft": {
   "nbviewer_url": "https://gist.github.com/f08e978de2c5aae86677c0463b050321"
  },
  "celltoolbar": "Hide code",
  "gist": {
   "data": {
    "description": "Generalized_Mie_Near_Field_Plots",
    "public": true
   },
   "id": "f08e978de2c5aae86677c0463b050321"
  },
  "hide_code_all_hidden": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}