From c9709ee778f83c256589c10c5f057ed6a809eef4 Mon Sep 17 00:00:00 2001
From: Homer Reid <homereid@mit.edu>
Date: Thu, 28 Sep 2017 23:03:01 -0400
Subject: [PATCH] Arrayslice (#96)

* trying mkdocs stuff on master branch as readthedocs seems to have trouble with other branches

* synced with stevengj/master

* sync with stevenj master

* updates

* added src/array_slice.cpp and tests/array_slice_test.cpp

* added src/array_slice.cpp and tests/array_slice_test.cpp

* updates

* updates

* added real_valued and complex_valued array_slice entry points and added sum_to_all processing

* revamped tests/array_slice_test to eliminate reference to h5 files

* updates

* updates

* added libmeepgeom/cavity-arrayslice_ll.cpp

* added libmeepgeom/cavity-arrayslice_ll.cpp

* updates to libmeepgeom/cavity-arrayslice.cpp

* updates to libmeepgeom/array_slice_test.cpp

* added python/examples/cavity_arrayslice.py

* finalized python/examples/cavity_arrayslice.py and typemaps in python/meep.i

* updates to tests/array_slice_test.cpp

* added libmeepgeom/array-slice-ll.cpp test

* updates to libmeepgeom/cavity-arrayslice.cpp

* updates to Makefile.am

* updates to libmeepgeom/array-slice-ll.cpp test code

* added libmeepgeom/array-slice-ll.h5 binary reference data file for array-slice-ll unit test

* update .travis.yml

* try revising .travis.yml to avoid annoying link error that has reappeared

* fixed missing file

* fixed missing file

* update .travis.yml in attempt to decipher inscrutable travis test failures

* added array-slice-ll test to libmeepgeom tests

* removed array-slice-ll test from libmeepgeom

* removed array-slice-ll from libmeepgeom tests

* deleted .travis.yml and replaced it with version from master branch

* added support for derived_components in array_slice

* rebased; restored array-slice-ll test in libmeepgeom

* added array-slice-ll to CHECKS in libmeepgeom/Makefile.am

* pulled updated .travis.yml

* updated libmeepgeom/array-slice-ll to find data file correctly in out-of-tree builds
---
 libmeepgeom/Makefile.am              |  12 +-
 libmeepgeom/array-slice-ll.cpp       | 249 ++++++++++++++
 libmeepgeom/array-slice-ll.h5        | Bin 0 -> 42464 bytes
 python/examples/cavity_arrayslice.py |  92 +++++
 python/meep.i                        |   6 +
 src/Makefile.am                      |  18 +-
 src/array_slice.cpp                  | 482 +++++++++++++++++++++++++++
 src/h5file.cpp                       |   2 +-
 src/meep.hpp                         |  49 ++-
 tests/Makefile.am                    |   2 +-
 tests/h5test.cpp                     |   6 +-
 11 files changed, 899 insertions(+), 19 deletions(-)
 create mode 100644 libmeepgeom/array-slice-ll.cpp
 create mode 100644 libmeepgeom/array-slice-ll.h5
 create mode 100644 python/examples/cavity_arrayslice.py
 create mode 100644 src/array_slice.cpp

diff --git a/libmeepgeom/Makefile.am b/libmeepgeom/Makefile.am
index 819a774ab..7a67b7485 100644
--- a/libmeepgeom/Makefile.am
+++ b/libmeepgeom/Makefile.am
@@ -11,9 +11,9 @@ AM_CPPFLAGS = -I$(top_srcdir)/src -DDATADIR=\"$(srcdir)/\"
 libmeepgeom_la_SOURCES = \
  meepgeom.cpp meepgeom.hpp material_data.hpp
 
-libmeepgeom_la_LDFLAGS = -version-info @SHARED_VERSION_INFO@ 
+libmeepgeom_la_LDFLAGS = -version-info @SHARED_VERSION_INFO@
 
-check_PROGRAMS      = pw-source-ll ring-ll cyl-ellipsoid-ll absorber-1d-ll
+check_PROGRAMS      = pw-source-ll ring-ll cyl-ellipsoid-ll absorber-1d-ll array-slice-ll
 
 absorber_1d_ll_SOURCES = absorber-1d-ll.cpp
 absorber_1d_ll_LDADD   = libmeepgeom.la $(MEEPLIBS)
@@ -27,10 +27,14 @@ pw_source_ll_LDADD   = libmeepgeom.la $(MEEPLIBS)
 ring_ll_SOURCES = ring-ll.cpp
 ring_ll_LDADD   = libmeepgeom.la $(MEEPLIBS)
 
-TESTS           = cyl-ellipsoid-ll
+array_slice_ll_SOURCES = array-slice-ll.cpp
+array_slice_ll_LDADD   = libmeepgeom.la $(MEEPLIBS)
+
+TESTS           = cyl-ellipsoid-ll array-slice-ll
 
 noinst_PROGRAMS = bend-flux-ll
+
 bend_flux_ll_SOURCES = bend-flux-ll.cpp
 bend_flux_ll_LDADD   = libmeepgeom.la $(MEEPLIBS)
 
-noinst_DATA = cyl-ellipsoid-eps-ref.h5
+noinst_DATA = cyl-ellipsoid-eps-ref.h5 array-slice-ll.h5
diff --git a/libmeepgeom/array-slice-ll.cpp b/libmeepgeom/array-slice-ll.cpp
new file mode 100644
index 000000000..0040669ea
--- /dev/null
+++ b/libmeepgeom/array-slice-ll.cpp
@@ -0,0 +1,249 @@
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <complex>
+
+#include "meep.hpp"
+
+#include "ctl-math.h"
+#include "ctlgeom.h"
+
+#include "meepgeom.hpp"
+
+#ifndef DATADIR
+  #define DATADIR "./"
+#endif
+
+using namespace meep;
+
+typedef std::complex<double> cdouble;
+
+vector3 v3(double x, double y=0.0, double z=0.0)
+{
+  vector3 v;
+  v.x=x; v.y=y; v.z=z;
+  return v;
+}
+
+// passthrough field function
+cdouble default_field_function(const cdouble *fields,
+			const vec &loc, void *data_)
+{
+  (void) loc; // unused
+  (void) data_; // unused
+  return fields[0];
+}
+
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+#define RELTOL 1.0e-6
+double Compare(double *d1, double *d2, int N, const char *Name)
+{
+  double Norm1=0.0, Norm2=0.0, NormDelta=0.0;
+  for(int n=0; n<N; n++)
+   { Norm1     += d1[n]*d1[n];
+     Norm2     += d2[n]*d2[n];
+     NormDelta += (d1[n]-d2[n])*(d1[n]-d2[n]);
+   };
+  Norm1=sqrt(Norm1);
+  Norm2=sqrt(Norm2);
+  NormDelta=sqrt(NormDelta);
+  double RelErr = NormDelta / (0.5*(Norm1+Norm2));
+  if (RelErr > RELTOL)
+   abort("fail: rel error in %s data = %e\n",Name,RelErr);
+  return RelErr;
+}
+
+double Compare(cdouble *d1, cdouble *d2, int N, const char *Name)
+{
+  double Norm1=0.0, Norm2=0.0, NormDelta=0.0;
+  for(int n=0; n<N; n++)
+   { Norm1     += norm(d1[n]);
+     Norm2     += norm(d2[n]);
+     NormDelta += norm(d1[n]-d2[n]);
+   };
+  Norm1=sqrt(Norm1);
+  Norm2=sqrt(Norm2);
+  NormDelta=sqrt(NormDelta);
+  double RelErr = NormDelta / (0.5*(Norm1+Norm2));
+  if (RelErr > RELTOL)
+   abort("fail: rel error in %s data = %e\n",Name,RelErr);
+  return RelErr;
+}
+
+/***************************************************************/
+/* dummy material function needed to pass to structure( )      */
+/* constructor as a placeholder before we can call             */
+/* set_materials_from_geometry                                 */
+/***************************************************************/
+double dummy_eps(const vec &) { return 1.0; }
+
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+void usage(char *progname)
+{ master_printf("usage: %s [options]\n",progname);
+  master_printf("options: \n");
+  master_printf(" --use-symmetry    use geometric symmetries\n");
+  master_printf(" --write-files     write reference data files\n");
+  abort();
+}
+
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+int main(int argc, char *argv[], char **envp)
+{
+  initialize mpi(argc, argv);
+
+  /***************************************************************/
+  /* parse command-line options **********************************/
+  /***************************************************************/
+  bool use_symmetry=false;
+  bool write_files=false;
+  for(int narg=1; narg<argc; narg++)
+   { if ( argv[narg]==0 )
+      continue;
+     if (!strcasecmp(argv[narg],"--use-symmetry") )
+      { use_symmetry=true;
+        master_printf("Using symmetry.\n");
+      }
+     else if (!strcasecmp(argv[narg],"--write-files") )
+      { write_files=true;
+        master_printf("writing HDF5 data files");
+      }
+     else
+      { master_printf("unknown command-line option %s (aborting)",argv[narg]);
+        usage(argv[0]);
+      }; 
+   };
+
+  /***************************************************************/
+  /* initialize geometry, similar to holey_wvg_cavity   **********/
+  /***************************************************************/
+  double eps=13.0;   // dielectric constant of waveguide
+  double w=1.2;      // width of waveguide
+  double r=0.36;     // radius of holes
+  double d=1.4;      // defect spacing (ordinary spacing = 1)
+  int N=3;           // number of holes on either side of defect
+  double sy=6.0;     // size of cell in y direction (perpendicular to wvg.)
+  double pad=2.0;    // padding between last hole and PML edge
+  double dpml=1.0;   // PML thickness
+  double sx = 2.0*(pad + dpml + N) + d -1.0; // size of cell in x dir
+  double resolution=20.0;
+  geometry_lattice.size.x=sx;
+  geometry_lattice.size.y=sy;
+  geometry_lattice.size.z=0.0;
+  grid_volume gv = voltwo(sx, sy, resolution);
+  gv.center_origin();
+  symmetry sym = use_symmetry ? -mirror(Y,gv) : identity();
+  structure the_structure(gv, dummy_eps, pml(dpml), sym);
+  material_type vacuum     = meep_geom::vacuum;
+  material_type dielectric = meep_geom::make_dielectric(eps);
+  geometric_object objects[7];
+  vector3 origin = v3(0.0,   0.0,  0.0);
+  vector3 xhat   = v3(1.0,   0.0,  0.0);
+  vector3 yhat   = v3(0.0,   1.0,  0.0);
+  vector3 zhat   = v3(0.0,   0.0,  1.0);
+  vector3 size   = v3(ENORMOUS, w, ENORMOUS);
+  double x0      = 0.5*d;
+  double deltax  = 1.0;
+  double height  = ENORMOUS;
+  objects[0] = make_block(dielectric, origin, xhat, yhat, zhat, size);
+  int no=1;
+  for(int n=0; n<N; n++)
+   { vector3 center=v3(x0 + n*deltax,  0.0, 0.0);
+     objects[no++] = make_cylinder(vacuum, center, r, height, zhat);
+   };
+  for(int n=0; n<N; n++)
+   { vector3 center=v3(-x0 - n*deltax,  0.0, 0.0);
+     objects[no++] = make_cylinder(vacuum, center, r, height, zhat);
+   };
+  geometric_object_list g={ no, objects };
+  meep_geom::set_materials_from_geometry(&the_structure, g);
+  fields f(&the_structure);
+
+  /***************************************************************/
+  /* add source and timestep until source has finished (no later)*/
+  /***************************************************************/
+  double fcen = 0.25;  // pulse center frequency
+  double df   = 0.2;   // pulse width (in frequency)
+  gaussian_src_time src(fcen, df);
+  component src_cmpt = Hz;
+  f.add_point_source(src_cmpt, src, vec(0.0, 0.0));
+  while( f.round_time() < f.last_source_time())
+   f.step();
+
+  /***************************************************************/
+  /***************************************************************/
+  /***************************************************************/
+  double xMin = -0.25*sx, xMax=+0.25*sx;
+  double yMin = -0.15*sy, yMax=+0.15*sy;
+  volume v1d( vec(xMin, 0.0),  vec(xMax, 0.0) );
+  volume v2d( vec(xMin, yMin), vec(xMax, yMax) );
+
+  int rank, dims1D[1], dims2D[2];
+  cdouble *file_slice1d=0;
+  double *file_slice2d=0;
+
+#define H5FILENAME DATADIR"array-slice-ll"
+#define NX 126
+#define NY 38
+  if (write_files)
+   { 
+     h5file *file = f.open_h5file(H5FILENAME);
+     f.output_hdf5(Hz, v1d, file);
+     f.output_hdf5(Sy, v2d, file);
+     master_printf("Wrote binary data to file %s.h5\n",H5FILENAME);
+     delete file;
+     exit(0);
+   }
+  else
+   { 
+     //
+     // read 1D and 2D array-slice data from HDF5 file
+     //
+     h5file *file = f.open_h5file(H5FILENAME, h5file::READONLY);
+     double *rdata = file->read("hz.r", &rank, dims1D, 1);
+     if (rank!=1 || dims1D[0]!=NX)
+      abort("failed to read 1D data(hz.r) from file %s.h5",H5FILENAME);
+     double *idata = file->read("hz.i", &rank, dims1D, 1);
+     if (rank!=1 || dims1D[0]!=NX)
+      abort("failed to read 1D data(hz.i) from file %s.h5",H5FILENAME);
+     file_slice1d = new cdouble[dims1D[0]];
+     for(int n=0; n<dims1D[0]; n++)
+      file_slice1d[n] = cdouble(rdata[n], idata[n]);
+     delete[] rdata;
+     delete[] idata;
+  
+     file_slice2d = file->read("sy", &rank, dims2D, 2);
+     if (rank!=2 || dims2D[0]!=NX || dims2D[1]!=NY)
+      abort("failed to read 2D reference data from file %s.h5",H5FILENAME);
+     delete file;
+
+     //
+     // generate 1D and 2D array slices and compare to 
+     // data read from file
+     //
+     rank=f.get_array_slice_dimensions(v1d, dims1D);
+     if (rank!=1 || dims1D[0]!=NX)
+      abort("incorrect dimensions for 1D slice");
+     cdouble *slice1d=f.get_complex_array_slice(v1d, Hz);
+     double RelErr1D=Compare(slice1d, file_slice1d, NX, "Hz_1d");
+     master_printf("1D: rel error %e\n",RelErr1D);
+
+     rank=f.get_array_slice_dimensions(v2d, dims2D);
+     if (rank!=2 || dims2D[0]!=NX || dims2D[1]!=NY)
+      abort("incorrect dimensions for 2D slice");
+     double *slice2d=f.get_array_slice(v2d, Sy);
+     double RelErr2D=Compare(slice2d, file_slice2d, NX*NY, "Sy_2d");
+     master_printf("2D: rel error %e\n",RelErr2D);
+
+   }; // if (write_files) ... else ...
+
+  return 0;
+}
diff --git a/libmeepgeom/array-slice-ll.h5 b/libmeepgeom/array-slice-ll.h5
new file mode 100644
index 0000000000000000000000000000000000000000..9d359db83c4e00f0f26c222cc99daff40af22f63
GIT binary patch
literal 42464
zcmeFZc~s5s7cYKF6hh{Z9K<Jc=8XFwgd~I{IZ{+A&GT?HiBhLQ8qkpt!a>Nmk2!PZ
zgdBtrG9~jp_I`ivy7zw9y6dj{$Nl}*?{|IH`s{jop1ohM*Lj}zdG_9?R$E$3YSyG<
z6Vcy)dU_%qkpX@A@BjV?|I}M~<lpcA{kc^5LQAj-f-ThMFSJB|{}YJ}h4Gex|Nn8G
z?W~CtMVgua`v2Do%S2{6{MEmH{vH3n_WvIWSWX;g^OvceKvyi-QemG7--=vA`}_X=
zT4UY*TECFL|NpzQ_@B98|5uHcs1dIf9B2tUSxbutA_YhP+6J1R|NTizF^Kf&XA;qn
zrlKK#8t}x1n$f>M*VfV|2YSN%23k^{E~*t~^88kUm%o$urC*&d5(V=Q{vZ1PZ0i{o
zTA~L3^(hhD@P7_OASDw2cdY(%&A;b5q0xW-?Ei}+|L69VjXvx5|L%SDpL&=QY~h%Q
z{ykqkh0X9ke*eFJ?j*t+o?_j9yU+gf+`0a@KmNa8|1T<Vw_4P+>GC^Jb1KaJQi&RZ
z+b&d2-1QW+X8iDN|M@BCeSq`zeV@VJh5b7<o&5}I`l${Wxjcj0n+<gpPS0TX{0q8+
zr#*uOxvJtI<7eRQnYy;w!>5owuRbt4<|#Zi+F<;w%~No^|2u4e{0Tg(=-z9@S2gJ6
zRb)EOP(z2y(j~JC%3;58ve6fv&G4e>^QX;Zx3ELTf#1tj&v5&JN3rQ^Yp_lGM+<eI
ze87YiHx`|~_7SZmu5rcBK4Qm#p(CbNeZ<D4gUYX7|A?Q&d$zRN_7RmYMuoIk_7VNI
zHtBw#)km}&*x0JF@B<!re&>N#uMc?9EUJEPMlDWSxA{Wut@pTKV_tgpoOigk@znA8
z#Wkq37!qYIdV@Nz@^ZzUp5V2DxNQApw;*idwBn2oAK*~a-LFdz{Qy<$`{$uzZHb9x
zTFJ0J4JE@q7fiT+s-Z*~d9w7Rs-fiGhnq8ewl<XPZ7=FQ#=D^;Ej_klc-MxK)qVSX
z@HnC)>9;-Ttckgf<UrD?@tGSNNH)KH?XmThw&a6!P0C~=ZONnIhp(O)swFXSc%)-H
zTqNn`ZkHD`<u{n^`g!e7i=VJ=>#47;<K97dT2g(7R-Z6WZ1(%PPGi;Jxkr9EbZM=c
zx_$0yi}{ACgMH^;xcs$~YUskB7t}X8t4sx3ER0VQ<_!?m6$pM73w~XMIJ1Pf@7x$i
zYXMi4fb+G0d%eH|1zrOMp4kHL6hVg~L6@t7PXFpA5p?V;=z3Dnd5fU?I^lIVBfKun
zgx9H5c-;zw*HJ3Gt{}Y5H-*<-DeQy!!oKhj_Q_6R-&6_vsDZGrx(NGhrm*j#gnigt
z*q4WeecC|Sw`25q9~XAveQhl4^ADSN-(L~V!9w9&+!xNtM&aD75YExR=W3sD&aMjQ
z?wxQB4+`h<XDy%8=EAv+7tZl#;aryp=lqA*f7TS^=dj@J({h~|PvBUWxHrd6SAxgY
zy(I^qJpxJl0OkIQhu|`A`lKOyAHdp>F%8|FE8uatzB2bmIcU2s4wnSngBazs!Erb4
zz}DOg-wO<G!^aitwoVv)6FO?$j`=;Z46d*4bKGv;6-b%CNOi}k6kfFauyS4P87Q!x
z->&P@O#C>s{foQqmDnUsdQCg<4JIyVWM$?55lsfIemOS%3ue{ntxWz_hb}#59t<w8
z!?st!G;@C)*7O^bkg=i;@Be!9{&34WJa!?;GBoE4HZJ$MdaubB9M*F0uHYq~@t~~z
z@Y@BSuu`XG>Kx5;VMJI|;i_jJ(6x<o{-UR~m@w|d%3Hbb@wvy#t95hU;^a;KjyD>=
z!PaLYte)9D#DM#LUeP<ALF1fuF`tjugR5)uzRS_t60f)8a=&)am5klKUm5?cv1FuT
z?7cpAO(Z@2?;Kh)v5BO|aI?N0S~Zc(`jfnCDmIqnwCZ<ckac6put`s5B%z*UMVqF@
zXIkq?E>G?CtLZdd$(Ad#_C~lgl4Q@8#gALkQ1WK-jwfmh9mx&5&667IH;}aMXFo6M
zgqCDSa3>vq7m*|^`_A=V{vTkRdD6lN|F8H&G(4uRLvs~wxfs5>vYpCgLG0s0x}8<Z
z0qZ;8C=jdKUC26Ob*!t(P_PGuaj`J(zOe3+;Kx_+>n+6DF2vOraJ&$3RS7sV1>6||
zkFx?VFM;P^fp?Oi!$Lt9LqVsZf^LHa9VZI91_?SZ5_E4SybeLa>tZUrPXE4c=Y-dB
zjqtkq3a|4h;dOr{?1Ow^UyKy?Ni$*J92WM`O<`X>7xtN0*mv&2K1>$&Wd~uO#tZwl
zfv}IS3;S9&o%eaMu<xCPb8t^M7h>U@{CjS)gmV-roGX!V&MJg+mm{3RWx}}}Eu7ON
z;oO!9=h#&^*YAXLK4?yQ*_-)}*lSC0)AcJBW6QVNbM`fHMw#O}xpA{)*#B;eF(VS)
z(Lgk-I4R!?t8esEYpZ<l^Xgw4WAg$q?_5@TPpvR~cSRgjSt!FM+iy+`T^Wn*iVoaJ
zh1J+|#Fu(?RWh!B6DfWam4?nmJ`b-ttiv@WLGyzRvv9`^t%8V;S=dS*l4<%r6?acL
zX1@A;Ahux(?RXrH@I2Vh>mZ)}KqmVIDaV0pyC^O!U>xv*alul?3DJxj4lxgS{AcE0
zUSPpIA(MH7DeC|d)&-ncC%7@I?QypZM>H+HcJ|kS#b|qCe8TBiXM8n(qDxuXGQ85;
zX0c!Kay)ipMx%tj-gsB{nC$UyUldK-(rf4PARMx;&Zt@)fe&L_thzWd1_#-REEjcK
zjjekR^*fd>$LQH#Tm9Iz4!3K~INxVtF4k_|yz<43JS>fEYvM1?#i@gj9r~gbiZyfJ
zO>!_@2(#DBURiQRiY5tU;r)CX$1t7;P2*@C^v@?hpruWIVJF9dmK+y8&m$agfpNh#
z#tF@*5^i|GJm4ntf-2?-X7<Dz{8<NZU|rxR>jbJZ-|#|XM{HIT>~0aW7}Lhx8#~L-
z8GGGYcSQ7V8Cn~+Db8)@fhnQYcelyCvFAG{QMbK**lf~yw~=YVXfb-$WCvXtexLEO
ze|=OOIt=bPUU76Sx<42azvW#9UI<9e5zW|$S10879oFB1w!OcOEH2xEe(?hx2YlFy
zZK7{&JM}6H1DYBS>Fco)dd|`b8gL{E%x%fy2YnjH=#?}NrtmtjVn1+~{X(NaiUaF8
zE{GThY+_u{l5s)_<A#yU12UKwtYMxojCsRD)&YjIF3|lf=>%tY4tICk;ec+wCMRy(
zkm5tN$4}AuB`8`r+;hhCW%%h`QKRaV<+!W3Y*o3lH=fkq6jZd!4>wh8`?C3AFqZ9@
zzNJsP43nM3UbA&q;rsWwOEZh)7<q2bwh=F~(8=OML{F<N*g0rfgAHC=(M<cp{ad0G
zY!&x0Ru+;CGhfGI%CdA&yR>c;-YO8Yr(85R@og!7vLFkm7t=V{zM*-zl-FSp`vEid
z3kx|8WO7{i$vB{T3E=`W#tBdF6K*)eJm41d0;l)H6JnV+SU({hpn`P)*IlF&+`Ox2
znc;4aNe#ExwivVsXI^TrKI-d)x7tpBVExzyO+jzVdOde^ygB@=Q@<4`X?^a(+Izlu
z@>q7qZpOiw`L@}Qan3Tdk6Yz5f7&Ya{O7Rk{jsU|!(@y0jR`qu=^^gwvsi)Nmov_G
zs*S=Q%fy!LE+xZ}<)3G6s96p6X~zuvxkNybo9m@eg$<b1|Ej5V8h~PVvS1KE;~1Du
z^H7)9p=v4l!68@3FD&CY;Ms%XLLuV-8RLRv#t9Wi2{+g=50Ei0h+&=}X-B-_JnI0#
zjYt>BzehTOYnD}ivGqJOs~UdSaNh!a?mshTXOI*x*-pJ4mA3>BRxCb}7Nj|sR*kC*
zM|oiW!p}eFn)sl!-KWZpV**k2qQg2{RV0q}3~OoGG#)>k>FYFmf*j567yQV6lY<d;
z7e`vXUx|;0I*b~SD2HnMm+@nFBmzWTIuU*^2D&veUF`<JP%8^NAEo6CM*Z4%TX=aq
z3>icgH^tL9zSgCA7#T_HFrNM3Gxm$N90!~^uIBs@4#;6#pnOR<!LTjihAid*&dduc
znJ4s{L%d;30O<gx&q)_B$|jxQ&7}*K9*bt<`mUGFyY-)k@i)JuiH<D5^;Qw`1({OZ
zwff<Z^yf=Z+Iv%3%Fkukxm$R2Zf{S#WPM9#(P%$B`|hj5oK|7j&0|B8Q0ExDtEk-Z
zZcievY2^68d|5hL&3W_Mx=948CboNGrn?dn<eSXZs+CZv`)h@|O&F|++2Z*8fDcG3
zbxPZ{c7nbh)sGUoPK47-$YR+e8prQE4_zD4I<#Xyc$59&Xj_T{P8=7OZYCTM%ecV*
zE#ZW48^R5PnFqKsFL=W|p@ezE#9O2T7<VIGpx`p;1ftx?CZ}AdW2dJdM?E_`6Aj+p
zb4zuZgUT;WLpm4FM}=#<=<0`#n18NO^{zB0+_rbat!I%-v3CExleo|eH;nk$qx4At
zezw<cvweRQ7NtcD{@OJGP171qx0oM+!7Y-a6Fhw&cmn$D`X>y;o%0v1o#zkUhSiGB
zbG<<CTo2D|FK6%<|1-Bdd@huId(gwR@o1>|L>6O>XdDyH&^$cep4Q<`_Jgh#<QIz!
zC=QmTQd|(F6ApODxS-aYa6%U2h8X4nyO|f<W}Yy$j(9@?>i~MsNEaBmfpmiXz3rE5
zEFX^|yYS>rqbK7Y?N3&6TPrj>kSw+7G!sK>zVsci(iSW8tICcC+F_~3Y@@i|3$dg}
zs9|=XGn%b=bE3AR7dk(>|9e8RKUN(Gc1ci$q2){G0V}%rp=?s4n}<#Mp}MZ;ydSNn
zK<%*|1v)b(gKSC1zIlr#f@$4byOE8?fw*mt9XAFJ1=;QwHkLDtKv_c;Ra!KT<}Ya;
zs@-TEDvpvLlvR^olue;H;LLGB#W=v6aY1<^;e;XG2{*VdCmvwJykIo*gaynS7B(Oq
z;3ew<>3v8iuy31|;O5m22cH>X9kOpQj;sIjU}5hOC=af>d-T&-G+!Xmp0RX1Rt<Lh
zy1f5HEDFEcb7JOX%<AQF+Tr9Zv=NW*R-W&G@=iOare0l$vYtsp#FJ*B*`^uN8f`C(
zoI5qREGiU@<5~|K2LWg%x-)L_El;eO_qdnGwk23Q(YojI9rI8eJ}TbA+!AZz9?bmp
zPe1%@K^AIaXdJ36X&yw|(>lm%Kz^Vmp8SIIJc<MIJc<jV9)tr#-3S+W_aL0`(}Zxt
z$S~poa|RJFh_ol3FffRC!xq*7Vptb&YEL@BNtc#A4nFM$hr4ZSV0A(a1>fuHmQ)%;
z^3I><_vna0av|~f?FJ@r^v9ud6<hm3^X~hnq<1z4kzd0_X3YnJYC^y2W#$9H+jMs7
zuI6@d`%q|DYK9+(I&ZR_cxEU*SXwyq=*u{qmeqS?!r^#K+P1f$$Mz^RzV~3`eft26
zs~qHHG{zktsUs7gKC{PbH=4}+*;0bZm&rmEk3-gDng{Xav<^NLk{?)?O@2Y_N^w9W
zrnpdDLO4*?k8q)47~w>TJ>kZ{9>fDvwGRK~h3d=16L0My-dO*WbU<wr(ghuzNhhpX
zY+^gT**HjxS?sahZ9F`X7JXdeVgV~+{aiaNpA63o#>_a-&I(dw`kh~Fn+c<=BJ92n
zo(H;*)5ni?um|^fxA!@vEQU6@U#y+R`@oAht+wy4%b>>h$pniJi6EbD{JQ?;EKoJC
zi7tp;gHE2`=Ol;3qj=BD#U@Rnut(P8(W!fV(Wl*{{(1Wsqd|p_<Ls0P=v_$`zVJ90
z^E`0aOY7h!`+-OU@(VHxih~{fDJ~9+Cmd*<L%7g>4&lTu4+uA|Wgb|~yl_ZM;)w;!
z8}G9Y=*qg_UDgSm8cM%To;e+^WuLBHGRy`>XI&1yHf9za-_Uc>js|n!Qj-D6ziaGZ
zbQ_(wBT^i}{pi7!<6V|Of^jR4FCUk}smjHtzdiMW-|g4F8$Krt&g`vyxFujU6qSaA
ziQ8m@vw5S_gA)SbtkWI4oiEby=$JmuCuk?*L95gW9d5*6KYwqFDa`|Me_YbkU&W62
zWk6q_o=b;eGeffQCYHwGPyo#X_vy3_MzSBMX20MxoZ?_I$HjJx17{f$F1%hrIPodt
z#!JitFEB6c*OYkTXXcIGtOLGeT`-Gv!b4}It~v@^Xrh(X;7ILU=)AY2-16>x&@SvT
zyG`CgP&6vXxJ8R$#X|kV3qCl*-fv%Rt##dDv`$!rL8cd^XnFN-+}R(#8E;>?^o<O<
z4eUBPGA<d0OkDlr<F*`_zgO#ehSz3jFm>AS?%QOrvz1Zv)Fa8L>;GfR&`zt+Cd1xu
z$>~5$GVea2_nhI7;_y;ib;S~9dy_?jfi#Zw=`;^|^E!xOKTybiu@T3?^A;2rJsAgH
zVO;1cBAj?!n{eY9=7E093pX%NtlLh!@#7xS0eiD9c#w5M)yVRc-m~qYbli}Yk1Q5~
zt@4D;26rix-b|9F?{otDZs}9XH@SfEul%Mn^4%cPuIG`lumWV;B7PNb_knM&yI+rf
z7yyU0mUlh1dL_hsD&0HvSu#{6tf@81&w)70<{OTV-v%MpDHD!8ECjDD&g-7oX29XM
zJ_Dmn*5l!vr_Xi{NJWQ>p*LLC%0Oh><Da#aUf?jFELQM1*77_k=XFq-PkvC$esMR)
z!KWM-3mFH-FfI&qBb->hop9q(=7H^*7mjA0xQ%(^RMr7qSQk9WI$_K*-|Ce8i(o_Q
z%{KE>7eh+o`t^YsOW^gc^={sJOQD_fH>U?Z-Jsc{8O9EiJ;7|q%8&uUKHzp>pJT*y
zUuf~fVbiT<{xE9QwnVGuk<hGZm+6L2R>Nx7b5VP4u7i6vS9`Vpr2xI??WKdXw!@s7
zPb*T@+d=hfb6LxoJHgkY_MgFv70`Ez+;#6)U##t@>#lDRhU43lg*BEm4&8YkM)Epz
zXFu4J{o*r@gNYm$Com55k0)HH^OSJnOU8{N=7Gl^5-&Wlgm_{I^Tu%20jIJqXvaEX
z&eLzkK@LvPd~6fvV;0VUk3u7LmoJ55{~R_Cedh{=<?9oF{aFrsrg!l<_+|y%3Y3X|
zr29grp=QIhM*2Wb<Jyj3?+Mil11`7g<_8_-=2vx!mO;5;+jf`NCPCuI-<GSEW<vP7
zR<kDb%7f~Ak(XxF<iqNzbpihvE8y13n7r5=8O$Bkc+sI$3GSIFKQ*X!87#9Pi)LeK
z9A5K0RPZ{iVL!N({bC%)!SNj^F8VPJ)ORFYsNIxsVi4oTE6f99m=|U+Pu$GBv6ywh
zeyj`bVV%%E`DcjEL1)NnpJ6%f%u?`d6{BbV!xhd9+B_=6XE`jeh_&8l=>-G&PPlyZ
zp*K|9_V3-Rg&%m_`r;~f^@PGG*`o{c`JgnjbGp!RI;duiIMtz@6b$dIx%=dvA0)Z-
z>tSCO4XxKSKYOjmTKM8!yDC$b4y%42-!{Q83mUpS?`qRI2SV!f-Ui-CgRGI0A6K3Z
z0__oG(TT_LkTcE03SNhA*bf%2CchZYac~XC#qv3X17-aQ7b;^2Cl+2O+_;f>U|k>L
zg@c$Uu43Lem36>q)&*@?CoKN)&$ouvOCjZ4e!Ol6H<&&z-*>zHau_zn-9O`iCphb7
znr`Xf4Y74y&KwEz0ofDPo*o_jG`g70bIo1~vc&X%3<p`E)Wc}s_01ld`+wa#ryl-T
z=67?H`7eKT&RE{}!B!u<6uK%e=HgORW@YFaJfDX3c`4I}=G#L3pU<OCKU)fxW~vPr
zIy=La#KE?Oi>86qTC#Y%KaHb^=b<#8)?wi~@&oz7<QK}*C=S^9Qe3ECPB=iqxIhv?
zIKgxd;l|#~1Cy8+j$oc>!@Tj{T+#t^Sr>F-op6QGnWq;u`)SAZ{4@G<mP7u*7p9>F
zo{&~_?MjGx1=L(M2`Ul!!0mlW587(`f^wdg&DJ_!P+uN(_+#5eAgvJ{Yg+D(cD2sR
zsCF`x+)6+4wrv9HOforl?|d>!`hE;u)HMx%8O=SsY;XqlZiZdHFI<Nkzcw82YO@Z}
zJJonzi*!sJcX5SQk2N@>$asJCTOTxF3p;rnQg|NJ@j7^-Pkx{=`vq%`1M-O!7it*?
z*t{oPpngR-L3W97qndeOE%QRjHR6dcm^W%4A{{V=b-`%X2`?2-%eDIA4wdB%Mqh~2
z?5Dd4tzLd!0T*08=09oe1GOHTPvCDK*yg0$?Z-A>(D6)qlK06M6q7ctNK(!Q)eDoY
zXEz0);qj0s%MY!>A^A6wP7hB*^D#G^71>!BemE&&&+AR-9J}Pt^jY~B-g<4!0-pkG
z*t*DAba^YvPW0PzvF}#Q_+q_vuggZ%M|g2c6^(sWWTBMDA%^F{yH>Og8nPdl&3++?
z<A7>3#f9UH0~E&y7Z|)DoS>LMxG|V{U_JA~j-JF5%a}KMvkthKbwNAU2}h)~h_$`r
z0g-_(pR_x;0tWlnc6~6=2TXH3(jVIUf<f*%hfH%nFbMaWAS&>K&&Q*e%~<9Kl6$ST
z4Xv35N|(o<_nizyWv7IDpPQ^jLyy)|s=j8Tp3b1ei@g-+6{s^je$-YB$6$l+&3EE}
zzFnWp9k3gVV$2K0p@kTmHFm_f<8ln1*kk?9Deic%eT8L^TL8Knkp&wGje|SSgGam$
zT-Xn^V82itKykp4<AR+h;Q#~11!6bC36dv-8x2kp53FHcSj;@Jka^=C)&Z+n7xZSG
zaMaMZ|3t6zf(y!i$G(2>hM(W|C>uxlf_}S}bxZE}LG;=y`t7g!gXyl(JJ0qDfQO!$
zzXope2i1-}v+lI90q?4@s#EX6u=c~n#?$;#P#*TC$>HbOD3t~_jWF4QCr2lrDaqZ5
z19w1!D`%9rR#CHX(ZGFpWpw`dYhw<e<KPH~eqVQBLg)P<%DR4NVmEF=t6QO9J%=o2
zn9?{L9YXV<hS$MN_5)MdFXV6>NZ_~-%s4>&iEx4Z1>poy3E@WP6~qH?GcW9xOgwQc
z^TuA^Ne8^ax?l&^3AK&lttwl2!~ULsBwZT#!ttEUyk2TQcyYDjK)=-i;EPr3<<|ls
z{NQ50+wMVN*z%Ij<+?ynZEy4~QnCQlOHbd9C=W)_tUqH_8ObP)==ZTrb`B=g^qbT6
zO#v!RAGUa7ro^S<n<sGce!STZM=u|rkKaANVE+7dFm=OeLz5nI&>vNOdYAiJ?Au_)
zSljU-Sf@i4WacyuJ9B9s>UhyQe9wNsmi@vhjsv9}7i_`_2YA~OE)Xv!oS-Tr+$d^I
zJh06n;)U8O;)&mxH@3|s9k7UX!EUS*4w!xL%gBS?VDF@B{NRr-1l<^WdXi-T$UbDe
z>f#v$b2r%S?l3w8R92qdLft|^fAp7u?V>|Ldg12zL*>gr6&G(&waEv|=3L8(oUjJH
zdmsA#cySi$t@$>pZ_5IlJ1TcZS=A04;@0g+=9g3yMO@s|%}0}mcWs&)q=*6at(x*C
zk366>G28w~%u-C)v!ijveg_m^BMas}G!C*dnulr=T8H)F<Oe^ylV7MBNpZlN<ARED
zfQWH{ESzwH`V!&BvQfkXlb9DyWuEBIywQ(!z}c(|c43{+YI~7)4`UyA7?Yl6e8Uea
zjwN}fYz%~e?gclSEDeE3N7v*hKf=J|Wv8Bf+C_p#j%!e_2~l7(==GAFo%}#{JvQIm
zWdVAhS&|=E7LBf>a*u~jNXN?A&aKw{laG?yrh)Ue#Nu)HcmJ4<&w`SNY1{P2q-yF3
zK9|zJ#zFbEv5UmUK_J=RclZFX2T?OOErml5kZ6%b@gW+=;5?d#Hj`-`O4Z~C&AiAj
z>fNC@U};Nnp^R~W@-*Sn)FBBcnDrvusLCcDc>D$NLTl!U6PY*eU>&fMb-{tG6Q*2V
za_q<uABg{~NO08DiHF5E(Cc(WbMDgqbeKLd3~U3Y`#--I39%XH>k>Pzgy?ZDtr|u~
zgDf+8i2d6z5Jjga*XG!RYF@}O^P4^>{xMV9VC*V%Jr%Ze;IA~48R?z*K5RaEH#tz!
zy*?Gn6<Ws47pFjJQQr+|edD2Q`_6w#qQby+<=0OqJ1+;(T9*XBRx`kcExLNpI4ZW$
zJe2LFbtvCTey~=b{GzBo#R1DR6c?oRgagzI2p7mw2q!912{-m<9tg|}^_Zuoo=v>b
zlyyKS)&(<IC(LhYpua<NZ~xLxe}=b*KP2|CowCn081y^jR`k6T4n{+lyy}87&|8q*
zwZB;m*x!EGqHKK}=(Rrbb;N^c5C!cScKNF2{61g0v+2m0Aa0&*HP_u2WkFqHl+!(+
zWLU4Gdpo6rMJL_7Pp;{pXKquF-D@qB73^yhGcyj#bQhd9OAH6YzAIcobC0RL+djm9
zmMut2$zrvb#<8>o%|pW-v<_t>$PemwB)_N{L~&5!L~&8vig19ao^YYei*TYP<Hpao
zhzI7{6E7?uLOe0cka%MN>wqfO1v6PEoPOP~XP=KN;Po*t^B*t#;L45cM&-W2aH;8_
z?IGJD;G<EGP2w{v;lUPb#j!?l5c9a3wo66=WVCGmI{)e_P<_vc+5SEnq&>?!<qnMm
z@x+3ZG1r4ZcBNwa{IT(1xBIQ*BpW$+Z(j4utkqg5H2n~qH$MT+?rYU|r-|nIb$jv4
zlzEzceCznLEt9<=D=oX@%e8Yre3&eHT&HmqU!i#zzMIyeWe@U$W_jcnEsH1)mR+Q{
zC>u*SQ1TDqLiJa|iDm5xH-<0|JjlH881qDD=8b`@1D3EZ_z&xZhhA00>{olht=FNq
z=jQoBsPEt){f$9Tbz$EJ|M&<fC@=UpcuX`XSDib5*eV_(MWbh$HAn=9cH@@swMYVQ
zo%KzA?^y%VPNuWFcSr*H4&$ny@8w_+WxRhwPBJJPTJNuPj05GM_3u6<#(?>k&e&&2
z474eKo3<ib2IhOMTI{h4frLF7dhUl;K%vMRO*7}f=x;0jTGZQ0<ESX2d052j@XAN>
zqdAA<7bV&h2PHC!izSQ$y*m;vO<$dGVsJ~sjW7O;`pW~SFfVjuo|whF(U*0=E36Bi
zXPvOCb=!+aq^>aj&_9oz1H9pJpX<+Z4g|o6{!!Vsox;H_ysgCY(@J<f;F8Zjd*Z>M
zVqTQX@+7eO?GW%MF9jS{wVu8<PYz}0M`<m;oDTAi<1Az!azHP6Xv>w?GC-L$wf!zh
zBuEbX$2aX00tWjFl6!uR1oNq~TTN#~Kv|4#e7jdcP<m{0iadS=NLnW5y>2%bl7^DS
zhkIxoZNg|CDkso7blpsTP;`L&qCqOfL1_xbMK$9<@dd)Asp}9<JpYq$WAuIEfyvAZ
zdoxeG-H&)<A?tvnSr;^BoiKHIWQRHzCwMgN;P&Oa-QiyAYj?V)`NGvvp-US669VcL
zrIBXEGN^Z6``lL@2TkAlKk;=;g16i1bZ+Kq?oBHU8n!-_2Hw3~qtlkJ2SxMW(;BHZ
zf#vIY-So~DK(()={8&~V$o$`|ed3k?S8mEK_}7O+f%xRk#N}aNv$v(E=er<qMQz_`
z9UpKw9%u6Op&iWnNEVY0(>ONTLi13?>#+0^`9TX$@{9G}6bG$1F1j)f)L~rMj&W-0
zr-U1?)DRCm#k_D2^TcD!8xvRubYfj_8ta5-hkaZ#MRzV_e;!gbF~JGu+3F4GGtUDG
zOlNI&TjdWE$_6*CJ{b-&^|*=h-7%2T-7o9;nl(`GllfjYMh<U%y9_AyN{6}b?sZQU
z*--3}dhp`O&EOoEc5P4JT_AP7x6a|nMo`c3Pqdgc9&>thHW*sy1M{Ef?ta%R1X`?K
zJL<xj5V+NN@afxgHSZgv>N=!nNx{;BEZ)(hakPk_d03~Wb=a{R`9Zmi{DO`P#X)rd
z#l`cUghO*bBV72yhj8K*#*Kl@1Jjuox-m~&$GmYc>wt4v7qn%au(tQb0~?bEgU`&K
z&rI!R!j>gtch2j&7|a}W`*n-*fcxJayO-PzfXc+669Yd-Y0mA?VmKNPHkVzx_PCV<
zvnSlMTBer<R`X6<cJG-DFTUh;KVqK`idiixM(^GU^3`GMj34ZUkD_5l+jeaR)j#4o
zOA8rn2tNirFchk5EV`agjD+vm&nqKN$iP6ad+n+4<*-zXEDp@2aZGwc^ECM$t;5LP
z<OdzJ$*-nwOmVOy$3@3F!hvV!5H7sQII%tB#x&-E)=9(*ML&oqPG#Qsly$&b)&&h%
zC)B(DbX@l{(@<{GCilGI5EOl?NhuivAWA&8@W!|KkW_QA&SZxhxQ$t@|1~E7y7p_E
zXjv(PUJI_qow^hUu8GDOdmbjjow->aL(XNujgv#Z+!(5Wf`<*4+3qWZzFT584<E1}
zCOM}*8rbmw7>u)4G}^dL^E~_2@$>{owAr$>$!-%@bl(4I(6lddQ2&}Ve$3NQ`22}1
z8Xll=lx5O9&Ao-zK}R4zn8kk4gX5t7BZ`Zei~}p@5H1{@L^#onabqv$fhU+3_Iyh`
z@i+5E**DSwZF-U}P{ule(Or+3(>6;rH>_o=T30$?b(7nf+IDl$AizLx?1%BFetod>
z#wrP@F6x!PShg4rn8q%iHpv&NsxMb(W&}dePDAS-qr%|+@4&t#g{z?DgF}hwYt})V
zF>fXguh|SIPaJ>LHF+o4{;v6xe?$qp<Bnu*0440PUZgYEbsyASpWN`$twJ!e&tABC
zSSXaM&Caary$plP$U^WU8pn4$4;8!)onpxkUT42J={Uu~#T*y)DhLPaGA_KvII%h7
zMknThQs#xCYT}7D!NeO>tOJ;_E}&$cz@n?f>+&BrtiAcX@zc5PXxjhCgpB=eDC!e7
zz#w}WCPe+7`ov)g%ERl|F4wa|QR)1hQF<RviiY-nG40AsRQ?|GUiZ2yD)!`^DfgR$
z8R;^eQC@zae|~&|m0L3Oxv{;nrvC<*)XVL}sig&)x|FH+zU77R)oarq{}u(1ayMS5
zrZ^s2>m*BW4O<AY?SJMxF7}4uhGa3|&)UD^*q-O1lGkAp`@t0Ui;6OegAN=Q{TT;J
zju0*|V4NTtO}IhiN<7A@Lx>m1B8Vr{_afd<&pJT;Ez$*IyOB;%dL$=gO+O#JvO2p{
z%pPA1UcI&D<);1^+|FxsxHteyKGn{+eJ}uP{aTlu2n|5-;j>GA)cK*@{J>Sm!z(et
zF+WK9VJ*t<k6yPdIu(Qc=RA9GFa?$K8*Qyih{Yn;lfI677ozmZrWBi?NU;0$OZ?kP
z4#f-iZJTm;J!CYl@@~E)7xr|%J!H}LbTI6&{8apsU{G`<i-E>8j-n<s50e~d9qRog
zKN$Rs{Gxa^#et%Q6c-8u2nSeI5H7G;OgN!@B;kf&<^ko*3(P(fPbg*HP`HM4fZ#Qx
z3wR$TouIGn__9Ay{+QLe<AR>kgHSK%oy8TqP_%S0C_9}Kf#R19AZeuxgU1cNHfdBe
zR$KbGKdy>F=^?Y3^DE?NFz@5XyDM|C{*mhecj;ynmyd9Lq_qX*FM<y(JeH3)U0Y{8
zd8fc{jq}Ic9KQkeU#`}<n!66|K1@q)*?BFNFLThj7%oGDWvwQbW_zP5w)xqh0dr9-
zB@1$A8V8jd%>(I2S_k`IkRMR8U#R0apm<4f!GUo=LM7n>wJqTU3FC%x<^d|^1w)u8
zh?zGiSO<`{Azd_e4$?_8UiH{#YarTX78MQ42}N^HUx$JLQMlpCdr^n6(HJ~s(9ug)
zap)SSnhCcOusRTDi=QW9waQ{phw)ihR_O6LI%^9`U#&AL_1J-;GW`dA?RKMR;>A)o
zx4o!0F>>9<XM0iAAMz@i?!~(H#gT40yHUM*=H&*rcHsAvPOFl?6=2<XnMLUETs-)?
z?XEnRXbi6*3#Aeo2ft*R2WJy#9o#b{Kj1rn{DSEuiUV037nIit2beQ15bq$IV9B^a
z!aSflmUw}hd4kbZ;tgsq(gBp~Nf#&&CY?a(e|vYKLlBy`xp05&$Z)j!{``Pn<CR#O
z(yiI&q*x3tefH^Yspj)<@%@U!s@I}Yw<O|V>oioXm#i8*c@rv@Wo!?fQHbi6r_rTq
zFP8Wm9J=wwek{rDzr9UaF_x9+{U~jC2n$l3S1s^AglIc|LWe&Gv7$|GK;?~M%sU_S
z_{pUGcs*ydjeeU#v^u{g(O_*dMxP-I;XDpocpe0`rFBrre!!6ZLeUY517$fB7c3bE
z$XXFD&~qo8AX!Jap_X}oI*xdOz!S{m#2e~a2QX_(x<DE01PN!qy=@;6hza@kH*I_$
zhRT(W-v`ZJiRRr`Z2lY+hvn{@3twfe!NRsiXFr9@(dN7&H8?y2ZLEtn+xOUv(vtnz
zK~B5S+&wfgwDEov>1TB?jXQ|Sos0G3N0*?{N&B8pERSO8n}5C(`XAH$UPSjv8QRBi
zg7KMO`>YOP+4A^iiz+swd2Xi%qaGyU+a4CnKdp?zoxRCIRwj+ZS)K<nUI#Do$qy8;
zU%0|?ppfH2E#rVfH^K!5j1$ZmH&mY^9$?43K+QZs%Dh1yO*(+yR?-EO=Se3h_g?mB
z!7o3oUtc6|tO&v2O_iAw7Ru06KXYaC@>o=@c<(iQ<r*~jHplPOJvqvahR%rfT!%`%
z_`}QUH>2Wbv#^foMHrl3Vle5hMz8gGjl=Q|VY#u_>D@1mphsHTuqmI9qka9>uANVw
zL`4(*B{r4E@otlICQ;Uij}GhioA>lEcG`M6<XuuGN`rnBn;n@5-u2)9T1?(R<B)ui
z=E1ebv<_~vA2`l_p*WS|K%o)E1yjZW;xxhq(r<(l)VhQlWSxiy*whm*sAZl|&b*<7
zbpUl1>7vPJNGB-!^3VJs?^j^3V@`6$!T?mXYtXlgZ8+9%kPU}z(P+?Z{>*641XQXX
zy1agxf~sL&IrbMau&jQgH>7Vy?-i0uYpQpl-uD+97A@Y7;^gq$(f=I6QbS{hyAO|G
z)mWVRx-DX^-p&b;tyG9-`YiRGb`aHV`k3GAy$0oNC3syKgWm%>PR#F@k4e3@JB(YI
zh@S_O1&JGtgH(s+flXIh2hGouA6UbFp>zzz0X4@3OU3~PKL{7d87IiU5pD=(9-ur$
zyg<b~L9~&0gN-5S0IEf#3y8mxP9Q&*{77q_3o0+Sa2{Lch2GzPnI$R%(B_TJW&58I
zC_4VwYnXp5D$cHIR-B)R!Ov5^pG-`{>ioowXY>?k)+jDiD|-hj_x389e@%(iEw79h
zcPmCov98JN_`~R*d*=D%8AnicwU?Q>$sUwH_o^^jxE^KCrgc4@x(@XsD)L2-($Qw|
zn+svX($Lb^ykQ?%xaNDZP~t@6;3}tiptwitKv_n9!2CM-g|i$7lpGgKTN4f_QxYyP
z(<Pi>RzSGH>?QF4Df0r=QQ`@bSHv64SO-lVk8}a)dD00K$<Kl$UuUAEaM;K-3#2Hu
zRz1Iz<c>-=*E=3#{88GsxKa7;2(;0+xUgnW9E$H{&kH!V7Nyt6|L|>`i>k?Y9;E%*
zilWVz)Q^uBp-p+L=bo!dl-Z6PTr{8<MFTh8+OerX^S*g?=;Fp^C}{|5T_ka+j+Dnt
z3`;=obrsbgrmRBm311}>k4K~T$cRCT3J;WrkOlR48VBzTng^m4v<_sA$&cobVwC;$
z3-SRJ2Wn?iT&TZCI6$$HaDkk0g3VRJ4K{a)2h=ttULfm3JT><g;te7z(gDP*3#d{^
zCy*>SxpafA1Vn3ABphwg3PmQZO|3dwqqI{a=Sy1_qiERn5z*S-C~f$vUZ+(k+SH^d
zDt^VFY?|$Z@gW%~@wEB&BYZQ;j;?6a<MK9?Nv4b}-k{lU+mo&Oh<2jcpzuLjQx&Lc
zv#sdUl30}OT+yX+jW4R2P44JA%pWDMx{9j%`lI;4x%qE(y)~aFfbda=>`?WDtfmf6
z<C@QH&^(ZRn)G)a$O_00sLaSOC<am-5N)Nnpd3RuKyr+5f%kI438E>48|2IbRLl$1
z%o8L=#2aL+1Bf<}E})u2I-#PVu5+IZDX1!JtDbdr1j)DwY4KI_K(sT==Z*1H5G|~4
zasH{9=5yig4##*oqDbATTh1^|Uqf5HV(z}xC<#7emv}A<MZ4E-aG9pi+{4!T6nx!`
zsuRuw4WFf;Y(xD$^V^zw(3s0Td&~?0(TV!l^C4~^>i$*#a=Io??6c0j&kqBUQJ=W(
zqF$Q!*S$Zp=XOVts`#(f)ZYjE9S702s=xD4gtQJt7swCFdXZn0DJc$0N+>Rf83#yr
z5iU>#5>8O*5^ha>ig<u*4)N05M~Nqhl*Ajw7Ni5J9+EB?%sOFVk#_WimdnBB+N-&y
zS3NZK)0NH}E_s6Lc2Auna!-)y&2)4OS`Lz&$OzjU&F67ib<P=-y#Peh<F}`r@I}$P
zJ*q7S5>Y019O^kv(+^s=MfdEH6cm}J{b;E9T)AZUg6L^E@gQ>vJUnA!94Pj!-t=Pk
zN{}c<?P-}125JKnd+SM>&lfm8E`Dgg6eLwYlX~cBJ{R?qEUFsQILcmi|2q#=IkXN{
zNPbZ5L4HwUMsYNGbL!u?Aju#cAYxn~dPq1yWl6aI<$;O~#0zDOi6=^#H_BKClwBuX
zP{ulGJ~wyn(LGIF(>X9TT<+utrIpcLMpp!Yjl;*?2WM&WIMtC8##=+7Y@x-xxa_}m
z!|XSEpGSb~l;hfF*W5w0<K*HMwaZX6e|Y%Kp*|?`b}_YZl7g(X^x)q57*M2*H?y`)
z1XcZlyTf}Vfx^<$(mq7<{p9(}R?b@mC9US`7kRD(vpDDX&t?RGs;gA{UWNlmrY!nv
zHTM}BM?Eu|hmvu$4we1L52{)?{Pl~H%@hYE92Z4Vgabv6gbPF$2`7k-6K)iEppto^
z!71X2)eDF>o@E_yH0y%zSSK8MWalDp&F7T56<YV+ofrU?k(oy~Yz={`n<J*BFN%bs
zmJLt;c(xLP?S@YoSr-ecQ&SD(UI`$3apCpVm$4x8XcVyJwx$oc_EJ|f&F4?WlWt9$
zJtGyQr#~fjH_iaF$rhol3e%v>VMg+Zq;#;f{waS2>Ck3*`s7qEIq1m!M{1u;0M~lc
z&r?Gq!LZ=i#0K%sn)?%3EZ<4vXu6T+VXQZ;!`NBm2NOKWFN!TG4vJ1tT$D^D94K2)
zxKQLnI8pJIaAP_1!0LI#3k{hkN|`raiy$5F;dasmo3T!)+p&k&L``32QrwirpLG1e
z$o%!h&D}y^caMQ7-kwol8u0#|MP>}tc<O2!E=&N)v;&2<6O%RdzUQION2h>#R`5jJ
zy%xkh+IO{nkP2#(qgB~UvO#h;*yGi^El^e}E~x6iO{1UFlCFN6VMA45!PlqRu>6F*
z$sDV7Aa9~>vAHY-Y`$#~<;r3}-|A}TcLUu)r!iSf`#|IPvjfe;=I3c0o_<MwF#Ig}
zMae~qg9djfE{ggR4pbNsE|iWToM>>6aAWyY;(?XS3tKWz^kCj7WgYN(Ea`&7SSPfw
zEG_Wv<qQ?;t><`+_JCG7OYXVN_k*Ep*4Ovy9txAk)aok3R%+@ZrAbr9u7ZyDo4tO0
zF&Xr#cShu$)a;{7UC)l;n!Zrr#|C3#vp_z0-sYGb1*q%RPo40y5G=cy6#7MM0dbhp
z^78TxVCeo}%xU+H@N>>4`)8@yP?TNK#6Yv(t^^p@r+Kb|voTxd=UHj`YWvCJQZE|E
zA*nPE&)U&CH1#1rD2pP$SjKU%ww~f@-j@jnnne;WEL%@F(Wr`WYw9S(19R^aFMK?j
zc;b)M#2Z&-lMc9&b;0A4Nhj<mUfuAFt`VGZEDE#?od&JK@7krUb%Y3=MHlV<xWj{P
z2k#Ess_8?Fn@|`WCWCL&tGljSkN_XWCX;?NO9Atj!^4j0rbDoQ=8a@^HY8~4-d7pt
zgLHR=IDS(hDEsxCT@bblf)jpu$2TklgS%}42VdC&=Gr?vJEm`f+MN6@ZMJ5@8n0$@
zla`5K(Q$aURquVEW-3`cy_Ck$HI?R}s41;O)gAJK$_Vm{%03haRYepR%d`jw)~_I3
znDmNpVo_(pjat^k1Le#Mr*<cv_=S1nSk?iLEFxV{ZA?0$dQrwqGfjQ<iuLp4+nrpo
zBJai$14n06cF6N+cWME;`rAI%%@~Qx%N7T846}#6>l-+Bu=j(0<Cbi=@H`TVH{P__
zcOV|leTy?$RGJ1?=6Ti+x84NB*BeCb4cw;rqa9Z+hPK%Q%CB2p@{4wXY-h^1{u{SL
zty%Q<d84*~f#s8pEwZw}rrYtW=VFsV+`hF|m!AG$(}FC@`_VWWSkOE)=u7KRT10*{
zeG2l6X6Y0M!*5btESgI=F#0>;!b=+nCpKW**qM3YA?Ahke#BGrdxpdtY*`1eWnG{r
z>jX_I9(^>r9Dy|!-aTJcM`5_n%g^5{Wa#m3apHm%D^ani@u^N4o2eHtppTCXV^<gt
z(o=+C)&665`n(81quY;sXBh>c=)Co;WfO*KK4)L%Y4I=)(t3?F%HFQ&7j>MsF`z)x
zhnySJc-r+sFsmC<<W^h&-m|`y-m%LEvoixfH|(nE8*M#vxNZ76&F3e#u6P?40pbj@
z7`%(d(cl8jL&Y3gr}=yy`9X(y<QIn+Q5-zZaWN-~aA^9qgbVF%5>7BGCfs1bJfM<!
zK{4}$j?5bzSqDgBT_FBG=>&SI#b0NoCSd)+!14RD5^=-E?jthaCu3EQNp<#<Qt@+U
zp=WKMbS%8}tjoz3>(C*7TJNi8GO=vT>-*;Vn!4&A7x`|-JT!gva74&{O+73su+y?{
z>v6}Guw9{|Rj7WIpc^^F7s^5(rapL|uBi*JYY}6*5lVXN2R^u%2@Y*e1}whg4~kDc
zG9MmY3VIvPyIg9&1Z*~uMZ2Fgj!tc99_sa_byz4RKPZnPzbJc1aj?q`iVL!pgaZ=Z
z6D|-lPN-(wP|7^uIP(I@dg2LvdlGM`-9S2MKL1R*z<4?71ZF*Djv+y**nfle@PIBE
zXts3xWy{g)QNKRJed?1O&F6_+e(H74!z;N@+Ep~kN4aSImbT9eu=@7<8J)d%qO{xL
zfO1bI+C3~-azRt~)6<jtR?ORv<$K;;8WW?$qKz0Zw$o1ho^I{tGe(m)nLZDFJxz|H
zdeP!{u|q*r(OvKN-8hggJ$AOBZ#1Yj6!hq}-Wk-p$f9fojiclV%|nrr)}f*w`GJge
z<QI%KQ5-PcOK~-MCE);5#swu^2`8A|BHUo1MLfWA1Mvbi^MtCe#2e~eNC#-ny1>UC
zq!Z|TeD`E<K_*_&t^U<^Weys&^4Zv~x90C@EDrVmvMe8E+Xh|KF5ilu4-Cj}tG^Qq
z&E7f|4Bd+cBeuuf_bEoD%ixt!bC04;yXxN$E}q2tSv%?$=c-WJc%g25feK3>Hc_6@
zLbU0UacSuBqj>vnBdy^^2T^R)Y~$pz-J0{(!f~ThJ{lP0hi4Y0Yd$BYGTe|FqUmpt
z1%u`^4#7n<53Y2fb)d(7z`>aOf{f#UqzT0ZgV}@w%-Rz!5G^H~U}sLaLCQR!hIxU}
zcH#+D%o|Kt2TlG&x<DJ&2_9x!U)5f&`MV>(J9N%8S770h<LAyy--2h??RWcPv>i=K
z;$Pmb*oER*1FTNJP-5A%UaJo8DAx2t#x9Oh9>c7>yVKi`IgRo@BTH(Yo<pVU!o+g-
zQuK~|=AoZ*5l!zui+<VvA_iwvwdgYE0+!C}9x-j$S?vANc++oF73PG8AK5kPFg~>$
zGgnW&3;QWwAHKLT4V{X}LhH#i4nIO^p5~rR>p*EoejtJUf~1_{Ky4Jo1slczO2!4Q
z3c?As%Lq4^SrHE~VO~&rn0P`d^M(S}0qR*7P_a%>e#P;dvs)g{oxd>p?%{m2I=D+~
zlK(c8qVK`s;$3JsrlI-FXw7-EzkGY~_+qqM@S(H3%@Gv)Zy7$V`6(>Jx*;PkpGULd
z!`gR0e_8W;@W)m?OS+D-PJ<<Kbs1*GcCCEh@&<~Y)DN67qYS<K2P)1_x{7vwH^&dL
zFU29_hrUW1dKzon^tkdl{1AGCwdlIgeG_{7lGXgZJ{kvu2{aGNd}$rX`;i|AR+C@w
z<~UHzaY1o_aDbF?f%p{R1eFir25;s8vCIpGFi$XK-Vn<=KyTIsj<Zf+_PDOsDwEA<
z+SeoTKwJTi?pW!ZcA^l6eAxefO0zv^a(~*r;-&jhF>5-M9Xo^uN7nj8*&j#M-L4kV
z0cTLwJ8f?Gs7si2<dS=(aT%6+Ef~-`@+L~%d~Zi;-$v8JQ~bWlZewAnj(d&zmZl#P
zck^eLn`pRe*Y3)KYiO{!q?yII3zz^G?+(^GfzBSU=QTUA1ucEZLK%-k;U<~~%AK?h
zB);Sas@X4?9ilj(Tt;z0%s3#Jae)Ejgu?rT8>Gwwf|(bn%ZVqPttH-2(1CP->OQ0k
z)K`*DFiZRB{`#uTs1EmNRui-pOTL&qob_P`mZm-54q2LhSo60pA1*zBCN0d8Z#6B^
z<kjsPcp08V%S@eW)ycDH_RwJ1y*5{{G@x6TZ?kS_^45lXzI?fb@^@F-DU0u5*`Yxz
z$1l8#lKR*;oy+c^P2TKrdm7xv^3Q)>8a29scG06}&-#8*lgG*0dreZIxOxB1KYMLM
z)5c_>d_IkXO$^NgLtY1#><3D|l3!4B9PoZiaY1yNa6lR30!0PkgpzfH8*1%{2beK0
zFlC;gV%{KU9iW(XfsU*bI0<^0MzCJD1--n}S+7_@FR7rHQqZfK<6>F4pjRO4B@y%z
zaQ}y1dJ5L7<O}On8_#+rED-c+!+M!DWWBNky=1Yhmz>vW`Vr&@^aQ=U1-)toy=qmg
zSAw9InV^?g(94i{fQO)$fuNUM(Cg=C+Bc9P=vAsk`wGOudF(5kM=RkxW(nspTsV*Q
z!g;JN;qz#EkI$o6IFEWq`8>*m^C%Y1qi{}86~gDyTsV)xe6Fxe!RJvvmd~U2GCq$M
zu6!P)qi7t=c%J6-uzVgxetaG+!^kg~3+K^HIFIJSc~lDL(Yp^lH^7pgBcNs;AQ#T#
z5atOU%o~hY2PhTJqlk5a(O3C-UnD&5%Z2B?OnBZ~2+w<G;dx&uJns{H`FU?}lAreq
z;dw6?p7-S)`FU>;$Itsh;dyT(Jnv2Y`FStl=Qv9ExsEbFe%{}X;OD(*JAU3<2+#Xk
zexCqK_5(`cd9M(j_ohY^7i7i!yjKmQ_aPAT`x40ceF{p1=e_Vg24dlP?;$+zGlb{;
zEa7=C6`uE!9KKKJ{Nek=PvJgcmdW=Csc@fA2=@sKejmmxSH4f^@%yyq_xAWc;k}IS
z6JmZ}$J%>*pRh^i`$Ua!pExeuCnksSePYK1zE2qJ;QNHuUcOIUo67fz(|qp$M_vae
z!hJ%xr+}ewpO6Um2{q$@GR6f4!hJ%@_aLa@dl6U&_lY*l6N-iVgj~2!lnM8VA*>S^
z3HQ}%;l3&n?yEJzeO0)}p+Pa<SCztj)pQu&S7pw8Ulkkv_xDwW9^Y3>D*3)zFXsEI
z`8~d``V03}vmJb2Eo;d4Rf7*m_f_d@vQXBT@2lQ?uZ3z}2gd9NOojWZ`Vqwe8@?xl
zIp3Q#`8(fNl}q@(YB`+mtL4Id)s%U`ZN8U7tZ-i~ZAtfbs2A?5w?ELm9(oJ;P_U2>
zSqk}(k&q7wc>>l5`A~t756SOxJ`^nELwZ6!6fERJdYqS_cOK_Mf$KRRGUzYlLuzfx
zhg4NUKIFmqkR*-sA;oddhor8Y56O6*rrxdfcO5FsI3J46;C#q#1Ls3ZAs<pV=6onY
z$cGYye5m{m=R*qS0XjlHR3_v@Hq09&tOJ+``A|pJ2|_uK#P34BtuN%;jX3Ybr$WB1
z67p>eA>USx<a}H9i}UT^8qT++k(_Veyv+Hw<R0hSrqP^lTQ1^!yS^9a+o$z8-?sa~
z`F7TN&bJHqa=u+y%lWpdIp^D^rkrorFW`K;UdXo<m7H%ky3YBwVL0d8-mf{|wiNPh
zvvGtQ>{e0T0B#HUwvLc*f95)arq9j!_E{m{-X!GP9zs3GYa-WkjtljiRzf{TB;|U}
zaiN~$>c#b(hX=WyqjQJrIpP*GsGc)ksOMPr;d;)(PF&9^x8-_H?Psp%D7h{J2A8>>
zV_3}f994l(&pF8T9MLwe=cvwcJx4v4>p5AUxSnGx)N_7y;ChZb$HkG11H~0w&oQ~j
z^&G2xT+eyPJm9QQ&(UU{;3(8{;)QxnnNZKsW}V=(P%ktV>V-~1z0g9a7b@JjUYIY`
z3&&=1z0h<l*9#Bia=q{}*X1CZ>vZ6MgX@K*LcOq->v~YVgX@K|-dry%o51x#@8Mi8
zw6WuQVQmMl7b<xiy+?7qP|=_3h3Y**z3?a33&k7<uW(%KJBRCqSwg)~k;nByrG#+f
zSfO6{>pImzQO$Kx&Ap!Mg|$Mx@T^cTEEej8vpRA;db&`Lwu|6;bVx(4M{9GP7Dfs6
z=<Py1y40BK(Q2U{o#4y$Xj85OgNW<GAnV2T=&~rTN7r8EdbH|CA=RTTHgP>#9Ln`*
z^I=?%4!7cZ^vOuBM;9k^Jv!?xtwROZwNWpL>(OQdgnBgB!O=#<_2^PRu15zmPV8yK
z_2_%d1GftGXgi@E9m>3MBI|&sg?e<O4qT6p7Wy5AS={e<%XNP67y2E;gnma2_XWU{
z`vl<niu)akEbezmLb>1JJ&^kyBJMjtLRap0sKdG6QSHtB4p#^6cUYuxzeDcN{f<dR
z-0!Hn$^8y@U+#BIisybu!4&RysDrrQQ7ZI1?7C1Klneb11ID5G{b=rYsI$4>@%aV!
zJDPLf28VJV2kn_BdNFT2DfBxetP4t5C+y#e`)OU{xu14I=%;P$$^EoD{@hRV^5TBl
zXrZ4b68dS9ncPoPcjbPXxIg#P%($-winH8Lb5(Iat-gl)Y2lZ-pH|n3`)Q3{3H`L$
z+)oSa%Kfy0I_{?x{CWJh?}jld?x*E(Uk?4*ujcas+)oqT<9?c=ANSM5-MOEpvgLl7
z`W*Mu%B;Abc8Pgmw9rq}XWp1E^wZ3^&xk%kKds+U?$?go!Ts8-DDKym+~R(%Q7i7(
zDusTn>uI51`<(l=qA2dys@8G8R>gf%keuRvt#~x|YZI1pzqW85^;tEa=jMKG(t7UK
zmYm{#tyIPRT7$FPuQmL_{n~<9?$;*G=YFk>`@WjbCvv~`zjSrx;W6Ls8vj`iDcYu%
zh^0Ezl3LS>iZ(oBPt_6xp`?S5Ac)vv9~wjjBUED_Ya>)c9gX#wdeOY98hezN;kD$g
zrBkI+>lxqsKIghl|C;NXYv%WTp6@@}%zb}8`R(z0twZs9tp~pwt8+eouN{ZqYaO@n
zdu<{1(9FL0z4ifruf2o4^$dQuRvYd`zvl0?o&@9pawj4Wkn2Srz&!wY0Q({20Rm%@
z2e4m79>6-tJiu1u0i4~L2k47DfaehM0G@@&1Grx?zuWu!9n1r?K^`E~hCD$2BIE&F
zn@kSCpL0eYz_%EA03SZ5_D0MDG(#T1Q4Dzi&l2PTES&5AI|}ju&Th;D{D?e&^91q$
z-Wz5QwP7#yU{5VT4nmJ#G54TTnFqMWJV5Ra<V`|PBX8on%)H43<V|d^B5&eZ!@NmZ
z<V`#Y%$wXt-o&vKc@z6$=1u-IxexWci@b@IjJ%2MDDx)skvFkDH#riqO<>-n>NAr!
zv8I|=wbe!5#8%SeQj|Fre;x;U6X!tYO<FQ<vY&aACCHoD)6HCBJ!IabAM++bW)F3)
zN8ZF;4tWz>f!SN#%$r2!BX1J(wYewF+l)Mu?FRBp&M4-Y4l>Ua!8}tJ<e5B|nP)1x
zV)9ImQOGlSvXEzT7b4H({0@00D~oxi<iF*amN3tB^A_?<c@2?gYWp4XOb50i&lF+a
z7t3=Md8QELnVk8^Gg)<+XUbrnX%F&D&O2rf@gQd<wsFWaIbN8##m+pF6~#Q$DCU`V
zAkXCOfIO4;U&u4rm}jyNW?rr_@^a3K$je!kke9PJMPANc40$=<ZRF*=-H?}ak3wF~
zz8`rxi+MR4^Ky3N*whn@yqq<RdAVCA_omJ~<mK!~k(aZD*-c)~Ru_3W_a)@zy!()s
zv+|Iavt}bNXEj1z&e<G!IXlk%zr38a8+kb!^KwqiEgt6OtP*A~v6`7Z{rOkq<*d5q
z9{l+<<mIeB$je#tk;imxMjq4o5%QScc;qo{?;wxqMs86<kz-VQKJu8(cFbe0VIC7X
zNX=&+)7^}D%-|CykLmoQ74nz^k;k-6Kpr!+CGwcQACSj%{xS7sE>pbikjL~aMjq2K
zlzGg1<~jb~%}mZye8_!@qZjg+)-T9o+SVhF8Ol7SJr8+I`&r~MxtCbXV>%O%$8>i_
z9@EV{rrm`+rVBY%bs^WP4(5Gr%=-p1?;CmqdEelh$oqOWAnzN69IZN-_l-f$R=v#o
z#;r!)*LNIwU(Ye*eM9q^_nnTsZ^97deM1wG_f2vl?;Fq;dEZEUPOQ_&`}$5G@0)i5
zdEZ?)N79)04K0tnuNQNPrx^0SnVXRJEyNz;!(RU1CXn~d$KK+4Vsg(C_!IKJwjYr9
zbrc~_UcfxLw=(kNIn0yW-$tI?aR_;GS99da3m+p-?oLLY+(K?!GZ!&WUIlsb2grHr
z0OrZd97dizr!w;7?y|^}2aZIZJii(8<mD<MPo9>5Jh|(Iyv&`8rxNnyui!Z{ka_aj
zf;_no=Zc+q@}$SelSg4r$$Nr4c>wnC|D7Lsau@a#7xUyc=E;4`lUvM_+s2`<;8=;i
zLdaV56+G-KxY<{5BiFA!<own53i=ApXXq;wv9A!rzJmQt^c7s_5$N_u=qp6Buke_C
zg<5;jS6H)=eTAjyE7(7HYWfOAtI=1Oor}IgZ+uUNv#)Rz&y`i^MTi|e32`kmy$Q+0
zT$1wweFcks1sD4Yk?bqvv#*fDzJiZ^h2Z7rD?}!vuaJVChsL1yp?T<mXaRa5nmGb}
znB2ehVG6pT4`XK^rVzao4QC%FU@Q7C7W*)V*@wyfTOTGT5q+5C4D?~T)<Pd9tpWNl
z%T}Qelf^#FNcLgAVjqS*8_7LjdN-2$w?0fDdOG6SV|qK1hPlPVK1>q(Fdp_{La?`l
zvkz0S5q+4fBk03qPC(x(AH5>YXWz=hzEvRmRyOvn*i%w#Df(7E^qBsyZ{<ir-zpNl
zCw<DkRT}$Nwb-}1z`oT1_N@Y!qi@x2Kl)ZBJm^~$!~60(`&PZ!w<==ass{U3SvW^3
zvTs$yzLg6-FmbbQ<#wTO6~ew%F8fvi*i#Cyx13?$Dib|4af~<jBq5W~=L^n6pU)bN
zKA)R?zCiZ*!qI!veD?WlhtTIM*o;13NFn-sMd;CKXb1H9($TZis_5Nml`81-z4O0)
zzI62Y+O<ZX@9JIj`8wcL?x5!<aqRQC@jc1GbEG!=d`ax{IoRjRM6Xah?DM(5LZ2^B
z&D`Q1Z+eLm!aiSNr0Fe61ooCJ^cp1#Jx57H?@^Yb2dR<hMXHs8zNjY-eNoR%^hHD2
z7foPaG>U!EoZje*MxvLgo!J-dUx>bF$zkY=mI^>$GzvXW&G`*|QO8j9MQy{-7Y&_>
zzNlv$`l8<9=!@ozKwmV=hrZ}4_C@=!FWLj=NHorsboND)*cZ*koDz_0dae?~zG$x9
z^k5}`ebEr?Em;%HJxC~ew0`ar`l7k$*=izsxBea{`p91PkptOB&Phfe*}V&W<h=LM
zM-FzNkDS9ka$XAh$cg>XM^5dBK62ho^pWlCBcEX(Igx$jcIX{Tr)uaUA44x$t+&ue
z&cWyW{b%%%t3O8{IgNef{y0bO;9SW?FIt@JBWI#FEjIM1#j)G$A%*NChhtCS-m)L}
zAc5>7??-Q2N}|WDY3Oxpek1hV?e)=j4?quG-ROm@70JGPB>L{2r|7%e_Mz_{oPoak
zv?J)de~%uz1hMZvfqnPj(dfHBJ&V5kpT*I4-%$m9_bBw<b-`ow-8-}IK97C(e(bxi
zWZ(Tu_T97Dcb~k}%%P55=(~HZRxf+`;{AYqcVEAkz5P~}6Mgp_^!g<fJ%3^EU$WSD
zcYzBaoyk+kB2OWNJOx`acna1X@DyCfz*Bfgp2E3X;3+iP3ZB9g@)Rb4V;~F2Q<x6U
zfixyhA&oqRUF0cDCr_a;7CZ&t3Gftb7r|2q17|_g@;kodF6jH+42MBW;W@g5Jcae-
zDU>5mVG(%>%g9ss6LYGaJcVlHDLBbfaNh(^!9kuvcv-`hkfc_IGa)we6fE*Oyfwk=
z2z>%xM@~9;9S6zlSlkl4j-w61>*(7LypBHLV8}1vVn~33*HN0hj`ieqM3L9wX$4+~
z^$fg@z&+q~oXP>OW84GqIyOE5uVXtn9{%6u!0Yh*0A5ELc^#$5>)1?Q$8(&kYsl-U
z4vvU62UkQ}lGjm{ypC?zOGjf*&HZ0qM-g}(f#8^ko4k&ICg8yg>j56j4f0@mlLvF3
zJeWJ=!89TdX7_RMV8*s152iPGFayYgnMEGV9B^9X^J3t^^aueDX5(1!V2)e?59X&#
z@L&!V2M;C~To}E-3OtxM?t%w%Y65sLWyyoF$b+%SgIPfyOe6AOGBKC7CJ&}H=GGzP
z!MsZzOd0ZE+K~q{nLL=MQRZGW>K$`W8cE*J0`i9XkvH@$c|+vz$U5?d4w5%?Fd4j|
z666iNHXOX6lHmOK|9t^(s5W^+?)Kmf-LZi;RO%FXL;b)Jl8?z7N^S<;(D?4)4F#10
zZzy~Uctfw90dHtZ4e*A%v%njwJ{Y{Ab~sntkT>MSTzU<2YGcf;+sPYBBX6i5c|)_w
z8)`=0P+RhbUL|j6B6(Kdk!O`fp4DdZtPYZAb%;ExLU5pDI=E0$cNBP5tH`qoC(mjs
zc~)L<rX*qvcvdMtf@fu`0G`#(jo?{*4sMmqSq7d}z+LdHO8z%^Rt?vIXSEC8)9^g-
ztPbP38a@U*tJ>sQ1(IjwsBbu1dY?S2$KY`3Hu9|2kZ1KRc~&XpSxq3%%HN%S$@$Vj
z<XJuF0$$qt<fV}lCda`Ilh)vfNj7<DyTKWg`sAhEAulZwTr$ZbFU<pPnK+h$mlnJS
zytMor@Y1~P!Apyt0A5<{ufR(y`2BOkOPl%ucxmZ5;H4cN3tn1a8}QP6Tfs|nlb7aM
z1YTOv0`Srb$xAB&$Ib8i123%{d1(pcrL7?^tq%6o?c}AcZ(z7_8bn^&b=;F)`WugL
zBYAuo<nbku#}`E&A31gsSP?wF{9@qoxr4#ub7z6a=a>l|pS=Wld`a(u#~1w=Jia>@
z!Q+eB1Rh`bD)9Jn|6@3Oa=fhJ@@e2s@c3Nh@r9f<-_xlNz~f750v=ysG4S{T$>U4G
z9GX)ZJiY)6JidbR;PEXdkM9lg_?nW(7fT*rSMvB)kjFQjJU-Xoc#qlSJ<cKTaW{F7
ziR3*#yA0lAb@Cp)Mc_T!$a{3Q2Jg{x0=&m0a2O?Rwc#>K0XU7ak-W#Oli)oj90%{w
zz5u+({2kyu+Va7B3^@(nV>o$_IpjTBo6U1IYCCz4;7qEgC3uhC5bz#DF{eH#58mTJ
z?4j}GJ+2||u@`xdS>!#QBJZ&$d5@>a6Ah{ip6FuoL_OeY%29AO<#TX1{r~+2o@huA
zc%rs(;E7rfz!S9s!4tKQ0Z%kM3_Q_QRlpOqfdeX$o4^wdd=)%Vt1ft=xmn<eT4%u%
z4Q~vdXnr7gqMps*iF)1#Pc)1?QCoBHL_O2N6SZRwah?WG)IJ?NQBQd@w>rrajVDiZ
zFL|Qt$rBwzp6F@vM6<~gJxgBg*W}g4l2==TyxNW9z^g4!UhQQ&c(qXr$g3SqUhPow
zYX1VSHkZ6w>se{Ts|^5G)}M1OW_YzeaA(E+HF&kwIPhxo2ZLAZTWel5bTW9gh4`En
z=73l0TWy}BS>WKRvjlmyhrp|KJqEAV<^ZqOngU)eIlCH2UTtaeYPVreJxyNiquGYr
ztJBD<T|i!K9C^r+3LdhxA3WrA@{mW7hwK73*xy$K57}KGJY;JDc*t(@kc(P?hx`k9
z$id(gt9Jx=$m1@8hn(9UJY?tgR)&Y{%>obEnVevF$kxy1Rqfw^hwS^td`^8IfrlK7
z=V<CM@Q~f)A(O)_!Q>%(Jmewo1P|HQ2t4Gc<RNb(4>^oH<iX@2Um*{<7VbqI<RNz=
zZ##~>?fK+wFC%ZeD|y=~<ZVZQQ!S1Q;B7l5gSYML0p50D1Ms#}!M)Zy;9#p`1bEw-
z;ABfmUGTO8kAk-yGZ4IO`{&?o`%=K$wvPmF+jRlFZCeO<+cDsFYv)Apw%-HS+s}nC
zoNvj;9OCO^<`Tyw@V1>lfwvtFZnzeZx4n?O?W5#v|3KdM9P+k@lecXr&%Jm(@Z2k`
z1kc@d0zCHz<hiGVdoCNxgXeC~2G89w2t0Q$dG5L3sB0E^?wRDd4=2yPVrB5$H<ITb
zvl=}2Uz>sFJ|h4;_sPCdFS+irvpIO~qcXvBFGrqx1@hc`k>|dOJoj+&++9<^bNBWF
z&)q$QJola8x!bXahLh(WMxOgD?5!Pf5Biuq_r~P8527!i)DHLp=7N(i3G@X-f}<~i
z+u#cbPJ%DM`8|9A*01mdIOq#VnG9dRQu+dlw}daC7kvRy^aW(o7qH?z_yWRb<``eV
z75V}?#KRYGCeFMsZRiX5lfHn(^aVuI7mz|<Kn#5WNti=i_uvaix(8o?{p<hl3*cTl
zoW6i8`T{o57f>1Zq6g>;SVSL1@I?41T=Y?_p^xGOeH0C+!AFq}kAt7zP{H^p9HZc)
z@KlG7B9A_bVEQN=^id3@kK#T0D7Mf?apOz)C>pkdkHTFAK8h{D@KL-Y@KL-?A4N8O
z6zAxp=uRKSuk=yyTqzhqA4L^2m){cuAB7cU=GN*8A4LrI(sT4t?5B@n1$`72?nT4t
zqu5E`$r<`ivgkXp4TbNd6FeYNd@y_`5%isSTfldcPv41mFMKCC^qqv#cM?qBNq_oI
zs?m4y0NxTEN#98peJ2G6;X4VV@1z%fCrjx&X#r1)%%JaNE`2A?-NvgT?dUs6T4%f~
z5_kZ<lh7ORomkD_JF)0H@zHnkK7A*T=sOui-^ot;PClXU<XJWNbW-WlSw^2uI(<5q
z=+o&$pU&mp^yz&6q4DW>ufwO4a~D3Hz%%gaMAD}dK%dTK`gE$$r!#^+oe=tT4%p$-
z*+rjDCHi!t>C>4?pN@w<owj(6tfx=sTm$3nk%l-|TxZQ3;-*i>eFZ)p?>;lP?xj!1
z+G#vN8cm;0$*+w^$p8OA@afE<Pv=wmbZXJp6iZ*z&+rt<<w*FN%F@@A0I!i`d*N$x
zbcL_UK8L=hq3|`OjfAf$9Nr`ips%U;6!@Bg*TC15^c8$f?Mgf}zNYjh@HJIj2wziC
zDtt{Y`kId8bDFsgz9uUjzNT%H;A=|y0=}k(<KSyb!W<GBXXX-5EPPEdz2R#LrLSog
z_R^a4HPKV0|DmsG8huS+^fmpm06ws5LGXe3-hdD6WBS1U(-uCkd+=sSSTFd%ww;F$
zEEt|G$*Bt;SX@u|z#Q~}&7lwMV|cry+!&wnfmMG2AJ_|czGMM?V6U}+4{Q{?V7iDt
zupRV)<<bY%70*$6#+2SMjinDPp*MVBp-tce3#1P${33i{ov?>Cgy&5A&<8f3KCmGA
zz*^D=R)#*X`t;4MZ3*ApYxK=^rEjh(eRJjLn_EWT+<~U>&GmrCP1^Q{Z>|%4b8o@>
zCe!Jg>p|aKsRH=sR?#>2)-3qu{%k?t96WRKV+Z)=UWbQH?$9?klfJo?^v$)QZ_Y*E
z+{ZXqYtlFOJ$-Y&((uhK4Ks7=9{T3`(Ki=M-&{O>bBXlLUBSKRNc!f^4~5S#h(5#G
z^cl9M&oCLDKj{tcpX4-!&+y+n;WNyKCs4kCH&7bWXSj(z!w2*kwxQ3kHGPJY=`*y`
zXZUkR_zWGv@EL|yfzPlPJcklSpWzAm46D*-czLGrB<f}Q45!d%ScX2sL$l#C+)STg
zPx=h+U=KY;pWzSm8P226up@njedsf+MxS9l`ZC|7FZ0DW@MVsnFSDp7e3@sT!k0OK
zzRW}LPRcoWC}k&onfvI=+(%z#E`6D?^kr71FLM@snP=$Bbkmnvw?BND75l)KIsBOM
zXlm1y#;Ym4?itUf9HuWbg1*d^^ko*%mpO;N%pCeMAJUhZkOg1n5&AOg(w7-dU*@;;
zWroq0IgY-}b@XK}rjNBGeXJ$uW35FWtA{?;_4Kj+K_6>L`dEk4$9j%FR(B|TtQLK&
zE$Cw%L?3Hd1^8IU)5p3Ko>lcmz{lEmIDD+*@v5G-@UaHL+bUD(V|_#)>zR*?=T&X=
zv0kN*wLE>SKhwwBl|I%{m|LqqHeOj>NgwM@`dBMsZ+!*#pjGH&T}dD7D*ApO(Dyrw
zzTX<F;rks153amZ3BKP7`{4VHw88iL_9gg!>D86O?LOoCJx<?m20Xm7o4(%+`hI_b
zw^ysf<E#DY``uFmzTbcb@cp*>48GrA;02bPzxjUC==+Tfhwt}1eZNoX`+chle7~_5
z;rq>~X6Dw(^!?7D@Am<`#(J5)-){8%7Si{dLf>!7Ao!#kq{1hCav*%t>t@3zJ+(Z1
z(nDv$Cw-kh>Fo3HNhe3cC!M(*KIz@5@JSbpy=r{Yp04mouYVmr>E-lEkD*U`3w_e{
zs=z1x`!)Ea-@6K*bVIx^<>-@6{u(~%lJrUM#&gBn5kBc%I9Kn(YpuKKlRn%VKIxrj
z;gi1fvGHQ-2KuC1(kFe6KIwOH54wXs>EGy+-Zd1y@-ubdD=#em$oR^$cEDFY^=J6X
zZ<I7%ZwaagUwNnlzH)DC_{wt{!&jaKPq^Q&4qy2-`pS>dSDx_|eC1K}m2bNSU-_OJ
z@Re6MWIW{>ea5^>3A`^hc+KU~S@_D`@Se+FJXiL|!&h#v1z-6t`pP{G=qqmwU%4wD
zzVeOmuIooHjE7xo9fPktgue1q^p!86ul!4R-F149@w{u>LoK>a+R;?IKYDm^K>ZdP
zHQspPb$4ImiPwtn8gIP5KiPQXwPHQvl~?Oc<C)jMw~TjQ^WQcediAU^UV81+%XsRw
zcB=8#>%!+}T6a7>QnNoXUVDw~Vm$YL{+sdMYvMZ>yFaVfNFA-rt6ZO9-k0O;&F7>;
zwE3Ro;5l-5j(M(R#+Y*?cda>BqC(6Z;;3Wh60gn7DQ?Uy2d9`lr008PFUjs__LQ6n
zW^X?~(A<MW{KMRf9II^ZNm3g08@#%pjrNS0y0BVAM?FydkIt^>&U$#~LPu1!?t1R_
z!-h>}M(a=OuIKNc5aa(lR?q6Q>krV}+LMEqPl(fk^xFq~RY$4)Ul*RYT{uq1eVZ~Q
z%`-(aZ^bnk{nl8uDwaC+U2VT_{zctBPbz%m_b_eQSh8D3KM%2L!?mBbm#i0eV?HR=
zTAU$uww!#onIs*%85P;Bx&$`4_43v4mox8c=pW{D>K$Ufr@n*cIchIi@8xs-o?3H`
z1QwceC1#hIL&{+;Nx+=4s*agkBCv-%$6j&~drIYR&EB#M_aOanFVYJ4B!8`(UgrCR
zF1qR3{SW_rxVyf0_x!NMLwf1Whh-<dxvQ^+CIrnYycDZlkCZChnlVVjzgczq-_wU`
zT+aD}I$)GWhu^Q0cXhl*_A6fTtDQ4+R>riZk7}l=XYqe^Y1em}?tdeFR@%rhlGjeh
zY^^s%?$=oV=>or}&HZMf`_0~?Brvbe$mOHrB%)@$-UU;8OX6=A>m<AqA)RlTSGCnJ
z@9R{2&d=8~-_s;K$Im%6&((+u<{XW~xqknfnL`RNmv}Iz<Zm}~OK<EUv#^&;z@Bmn
zdrLI#LB7Df$VuFj%>VI+%ZZ=%)Sy~FhSzG`N2{KGm|^YtMBSaIxjtFu(C_A53R`kx
zuvV?VCV$b>c)eD0>a3OZN9v&Nr%Qes;pcSUebm2K`6+4}xvI?MxX<*G`|_51+vn@Z
zQn44q<Ckb=&)dItE#v2lbqKl|T6(6uws=jr*E3UUyz*?cx8fApo4O(P#<yeS+@+Fv
zK@|o`z~sVzrPXXNxn0exrs91qoNYd*;bqMCv=GnH>K)B<)zjUaBZF|R+`}BwBhkzy
zt1zeZ$K3J<_K<SeOEzLpNx<H6?3B3&SzpuKi-h5xq|Wlrjcr}~=(Fq1&js)JMBg6Y
z{j=<e4t){b>E{`12m5;=1NUT1jMw1KuM}o|n5d2c51#g}Jw}Tv{Be18%L(ecUge8u
z&vZSQIOUDVvU9Y};f8~HwsmRxg^f>^-(I2Z|NQ3%oqwCDe>Us?VUziOPuM57=Cs-}
zU-r*9c3|rH&;0)Gka@!-MJ~;qa&$(qIEiXpxz>Om+lZ%`dDUE}d0)dyo6qSrd{58y
zGtbfIc&;|bIWhp}%CG-0bI72PW-e)X)66M>mCW4o9`+FZ*z6^rU{6_vz2)3Ea}RzV
zthpDNm1*ut#-0y|Tspj;#+LfZ_xgLE>Y8Jfy?c`fX}O>-7rR~=qG1xS;md&G+UivE
z!-E?6e@m@!<c9_4lXd9+<yC9dNYNGdUVYL!&Z)yw=ML)o@!Xfb@8}Dy{9epD8;YN7
zyF$Ap+G>s&=+?c*mmaR>o1wdh+9HFFEtJ1TJm|ja&OCX(bYF0zA=4ycQo^yx?kMq$
zes;$<sEMCDY+l8eY~I(x4D&fXJKTIvYsQ=BXyYHvbEOr|kwBa)O)-bu!CX=ubIQ%V
zW^Q>>-s~a8v6nn;W%iWQmCWAK0rw!M{>$8pTt8>-Ny=<(K781g{+f|o`@qSQ4y|xC
z@u_umuvQt+d}OyS@ml=1ZKVnqjnKKf^5+)BC8@Lcnu*)Cj@3m4n;)*bl%l7yx@>5$
z#HodYZKHaYo~xyOD~B%I=F-$7+c&qknW+=n{XVU3=M5S$F}}pAQJb_<@!v0;wifA>
z+4obOla|YZ+o`%Jeu*p{SfSL^@)Pvaj{gii@YhiN{3G)!!9&dZvh`Q<IlY`}zNf$9
zIZ^}9l{lOuAvjka{buHn`<P2AZZUJp&*#kC@)7p%b3M#na_^qmQ;K14v2YJE9rq#y
z+2)?)$UpwlA6|*o*YBhZZ}-_iUD5h{s~%(HH2Cxvd#ro$8d{^`u9ZQFI(~BDxy`O5
z4ZJtVs`brSo!$QEun#g)bmWIUI`%K^)SKOIH7=SlTO*p*ZR6g%K*KJE&Dn5oseV3w
z;R>tUYCYC=#n|u4ZO{?*AAPVZe}hh#^KNXD-#2OFnWr~@|DQ}P_iBx>hKV!fRQHk-
zo9vk(J$sl}?OxWruhuN{Ir$9Vli_%d?8S3sBhC>!&J{c6kk_KjTv9y3%qbTzw-m!3
zl8C*eKK7Jnr_A0m6Zas$-!}ImYjIDqphu<s!<NSCTX`oAKl^r|h8$><b|^JYkN0gq
zBtJi1)9XhK9J_0TZaUemR9eF%-B)^_^UCvNJ=kf(JJ};9YSl0QtX}Ayrt4<EKKGv=
ze5MgGR^+F(=IOijn~mMK#HDk4J}UoY>lb>kShu4M+OJe={koYqPp;KpQY+OSXIrVG
zo9^5Fu-IhH?tQOA|BwL^Kd<M^+<!Q<OuBiMVZF@zS_hw#tN5PWUuT{p+wffZ4ChE)
zoGTfaL%zgZV#A!0_pzB<u3`__kG=dnMYE@z!`{*X_aJZMUSuimNuEb9$-mt<Rs%Mq
zl>Kt_KrOgayw{F0aoR5S3s;R-<Mqnf*ITxjKU_Dz=)Nm`{U}|2vXQsAXN>Nh_iT0Q
z+VR@wXy%f-#inTCaofbpSx)Wv)v5NAXQui+&EK9}Eb~-aJS_5FU#R_WPq+PRMutw_
zHs;b#{lC!txfA=Q1ZQgY$zmI)-Otp!izY1#`E9Yj|IxnidB082>`3z}*YUpWz~^Ki
zz9&g|j&#Ozr4G)K7dTg<F^Am8Typ5JnNysfo4Mr__K+voOWr?c_LR}sTVA_s?m^1q
LUL>@axhMG_fgiwq

literal 0
HcmV?d00001

diff --git a/python/examples/cavity_arrayslice.py b/python/examples/cavity_arrayslice.py
new file mode 100644
index 000000000..bc97e9686
--- /dev/null
+++ b/python/examples/cavity_arrayslice.py
@@ -0,0 +1,92 @@
+import unittest
+import meep as mp
+import numpy as np
+import matplotlib.pyplot as plt
+
+##################################################
+# set up the geometry
+##################################################
+eps = 13
+w = 1.2
+r = 0.36
+d = 1.4
+N = 3
+sy = 6
+pad = 2
+dpml = 1
+sx = (2 * (pad + dpml + N)) + d - 1
+fcen = 0.25
+df = 0.2
+nfreq = 500
+
+cell = mp.Vector3(sx, sy, 0)
+
+blk = mp.Block(size=mp.Vector3(1e20, w, 1e20),
+               material=mp.Medium(epsilon=eps))
+
+geometry = [blk]
+
+for i in range(3):
+    geometry.append(mp.Cylinder(r, center=mp.Vector3(d / 2 + i)))
+
+for i in range(3):
+    geometry.append(mp.Cylinder(r, center=mp.Vector3(d / -2 - i)))
+
+sim = mp.Simulation(cell_size=cell,
+                         geometry=geometry,
+                         sources=[],
+                         boundary_layers=[mp.pml(dpml)],
+                         resolution=20)
+
+##################################################
+# add sources 
+##################################################
+sim.sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
+                    mp.Hz, mp.Vector3())]
+
+#sim.symmetries = [mp.Mirror(mp.Y, phase=-1), mp.Mirror(mp.X, phase=-1)]
+
+##################################################
+# run until sources are finished (and no later)
+##################################################
+sim._run_sources_until(0,[])
+
+##################################################
+# get 1D and 2D array slices
+##################################################
+xMin = -0.25*sx;
+xMax = +0.25*sx;
+yMin = -0.15*sy;
+yMax = +0.15*sy;
+
+# 1D slice of Hz data
+dims1d=np.zeros(3,dtype=np.int32)
+v1d=mp.volume( mp.vec(xMin, 0.0), mp.vec(xMax, 0.0) )
+rank1d  = sim.fields.get_array_slice_dimensions(v1d, dims1d);
+NX1 = dims1d[0];
+slice1d = np.zeros(NX1, dtype=np.double);
+sim.fields.get_array_slice(v1d, mp.Hz, slice1d);
+
+# 2D slice of Hz data
+dims2d=np.zeros(3,dtype=np.int32)
+v2d=mp.volume( mp.vec(xMin, yMin), mp.vec(xMax, yMax) )
+rank2d  = sim.fields.get_array_slice_dimensions(v2d, dims2d);
+NX2 = dims2d[0];
+NY2 = dims2d[1];
+N2 = NX2*NY2;
+slice2d = np.zeros(N2, dtype=np.double);
+sim.fields.get_array_slice(v2d, mp.Hz, slice2d);
+
+# plot 1D slice
+plt.subplot(1,2,1);
+x1d=np.linspace(xMin, xMax, dims1d[0]);
+plt.plot(x1d, slice1d);
+
+# plot 2D slice
+plt.subplot(1,2,2);
+dy = (yMax-yMin) / dims2d[1];
+dx = (xMax-xMin) / dims2d[0];
+(x2d,y2d)=np.mgrid[ slice(xMin, xMax, dx), slice(yMin, yMax, dy)];
+plt.contourf( x2d, y2d, np.reshape(slice2d, (dims2d[0], dims2d[1])) )
+plt.colorbar();
+plt.show()
diff --git a/python/meep.i b/python/meep.i
index 55627b0c8..a2dd0230d 100644
--- a/python/meep.i
+++ b/python/meep.i
@@ -261,6 +261,7 @@ PyObject *py_do_harminv(PyObject *vals, double dt, double f_min, double f_max, i
         SWIG_exception_fail(SWIG_ArgError(tmp_res), "Couldn't convert Python object to meep::src_time");
     }
     $1 = reinterpret_cast<meep::src_time *>(tmp_ptr);
+
 }
 
 // Typemap suite for boundary_region
@@ -295,6 +296,11 @@ PyObject *py_do_harminv(PyObject *vals, double dt, double f_min, double f_max, i
     }
 }
 
+// Typemap suite for array_slice
+// TODO: add (cdouble *, int) version
+%apply (double* INPLACE_ARRAY1, int DIM1) {(double *slice, int slice_length)};
+%apply int INPLACE_ARRAY1[ANY] { int [3] };
+
 // Rename python builtins
 %rename(br_apply) meep::boundary_region::apply;
 %rename(_is) meep::dft_chunk::is;
diff --git a/src/Makefile.am b/src/Makefile.am
index 5be09c990..3f5deb098 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -9,15 +9,15 @@ BUILT_SOURCES = sphere-quad.h step_generic_stride1.cpp
 HDRS = meep.hpp meep_internals.hpp meep/mympi.hpp meep/vec.hpp	\
 bicgstab.hpp
 
-libmeep@MEEP_SUFFIX@_la_SOURCES = anisotropic_averaging.cpp bands.cpp	\
-boundaries.cpp bicgstab.cpp casimir.cpp control_c.cpp cw_fields.cpp	\
-dft.cpp dft_ldos.cpp energy_and_flux.cpp fields.cpp loop_in_chunks.cpp	\
-h5fields.cpp h5file.cpp initialize.cpp integrate.cpp		\
-integrate2.cpp monitor.cpp mympi.cpp multilevel-atom.cpp near2far.cpp	\
-output_directory.cpp random.cpp sources.cpp step.cpp step_db.cpp	\
-stress.cpp structure.cpp susceptibility.cpp time.cpp update_eh.cpp	\
-mpb.cpp update_pols.cpp vec.cpp step_generic.cpp $(HDRS)		\
-$(BUILT_SOURCES)
+libmeep@MEEP_SUFFIX@_la_SOURCES = array_slice.cpp anisotropic_averaging.cpp 		\
+bands.cpp boundaries.cpp bicgstab.cpp casimir.cpp 	\
+control_c.cpp cw_fields.cpp dft.cpp dft_ldos.cpp energy_and_flux.cpp 	\
+fields.cpp loop_in_chunks.cpp h5fields.cpp h5file.cpp 	\
+initialize.cpp integrate.cpp integrate2.cpp monitor.cpp mympi.cpp 	\
+multilevel-atom.cpp near2far.cpp output_directory.cpp random.cpp 	\
+sources.cpp step.cpp step_db.cpp stress.cpp structure.cpp 		\
+susceptibility.cpp time.cpp update_eh.cpp mpb.cpp update_pols.cpp 	\
+vec.cpp step_generic.cpp $(HDRS) $(BUILT_SOURCES)			
 
 libmeep@MEEP_SUFFIX@_la_LDFLAGS = -version-info @SHARED_VERSION_INFO@
 
diff --git a/src/array_slice.cpp b/src/array_slice.cpp
new file mode 100644
index 000000000..b2656ec2e
--- /dev/null
+++ b/src/array_slice.cpp
@@ -0,0 +1,482 @@
+/* Copyright (C) 2005-2015 Massachusetts Institute of Technology
+%
+%  This program is free software; you can redistribute it and/or modify
+%  it under the terms of the GNU General Public License as published by
+%  the Free Software Foundation; either version 2, or (at your option)
+%  any later version.
+%
+%  This program is distributed in the hope that it will be useful,
+%  but WITHOUT ANY WARRANTY; without even the implied warranty of
+%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%  GNU General Public License for more details.
+%
+%  You should have received a copy of the GNU General Public License
+%  along with this program; if not, write to the Free Software Foundation,
+%  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+*/
+
+/* HDF5 output of fields and arbitrary functions thereof.  Works
+   very similarly to integrate.cpp (using fields::loop_in_chunks). */
+
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+
+#include "meep_internals.hpp"
+
+using namespace std;
+
+typedef complex<double> cdouble;
+
+namespace meep {
+
+/***************************************************************************/
+
+typedef struct {
+
+  // information related to the volume covered by the
+  // array slice (its size, etcetera) 
+  // these fields are filled in by get_array_slice_dimensions
+  // if the data parameter is non-null
+  ivec min_corner, max_corner;
+  int num_chunks;
+  int rank;
+  direction ds[3];
+  int slice_size;
+
+  // the function to output and related info (offsets for averaging, etc.)
+  // note: either fun *or* rfun should be non-NULL (not both)
+  field_function fun;
+  field_rfunction rfun;
+  void *fun_data;
+  std::vector<component> components;
+
+  void *vslice;
+
+  // temporary internal storage buffers
+  component *cS;
+  cdouble *ph;
+  cdouble *fields;
+  int *offsets;
+
+  int ninveps;
+  component inveps_cs[3];
+  direction inveps_ds[3];
+
+  int ninvmu;
+  component invmu_cs[3];
+  direction invmu_ds[3];
+
+} array_slice_data;
+
+#define UNUSED(x) (void) x // silence compiler warnings
+
+/* passthrough field function equivalent to component_fun in h5fields.cpp */
+static cdouble default_field_func(const cdouble *fields,
+				  const vec &loc, void *data_)
+{
+  (void) loc; // unused
+  (void) data_; // unused
+  return fields[0];
+}
+
+static double default_field_rfunc(const cdouble *fields,
+				  const vec &loc, void *data_)
+{
+  (void) loc; // unused
+  (void) data_; // unused
+  return real(fields[0]);
+}
+
+/***************************************************************/
+/* callback function passed to loop_in_chunks to compute       */
+/* dimensions of array slice                                   */
+/***************************************************************/
+static void get_array_slice_dimensions_chunkloop(fields_chunk *fc, int ichnk, component cgrid,
+				  ivec is, ivec ie,
+				  vec s0, vec s1, vec e0, vec e1,
+				  double dV0, double dV1,
+				  ivec shift, complex<double> shift_phase,
+				  const symmetry &S, int sn,
+				  void *data_)
+{
+  UNUSED(ichnk);UNUSED(cgrid);UNUSED(s0);UNUSED(s1);UNUSED(e0);UNUSED(e1);
+  UNUSED(dV0);UNUSED(dV1);UNUSED(shift_phase); UNUSED(fc);
+  array_slice_data *data = (array_slice_data *) data_;
+  ivec isS = S.transform(is, sn) + shift;
+  ivec ieS = S.transform(ie, sn) + shift;
+  data->min_corner = min(data->min_corner, min(isS, ieS));
+  data->max_corner = max(data->max_corner, max(isS, ieS));
+  data->num_chunks++;
+
+}
+
+/***************************************************************/
+/* callback function passed to loop_in_chunks to fill array slice */
+/***************************************************************/
+static void get_array_slice_chunkloop(fields_chunk *fc, int ichnk, component cgrid,
+	      			  ivec is, ivec ie, vec s0, vec s1, vec e0, vec e1,
+				  double dV0, double dV1,
+	  			  ivec shift, complex<double> shift_phase,
+				  const symmetry &S, int sn, void *data_)
+{
+  UNUSED(ichnk);UNUSED(cgrid);UNUSED(s0);UNUSED(s1);UNUSED(e0);UNUSED(e1);
+  UNUSED(dV0);UNUSED(dV1);
+  array_slice_data *data = (array_slice_data *) data_;
+
+  //-----------------------------------------------------------------------//
+  // Find output chunk dimensions and strides, etc.
+
+  int count[3]={1,1,1}, offset[3]={0,0,0}, stride[3]={1,1,1};
+
+  ivec isS = S.transform(is, sn) + shift;
+  ivec ieS = S.transform(ie, sn) + shift;
+  
+  // figure out what yucky_directions (in LOOP_OVER_IVECS)
+  // correspond to what directions in the transformed vectors (in output).
+  ivec permute(zero_ivec(fc->gv.dim));
+  for (int i = 0; i < 3; ++i) 
+    permute.set_direction(fc->gv.yucky_direction(i), i);
+  permute = S.transform_unshifted(permute, sn);
+  LOOP_OVER_DIRECTIONS(permute.dim, d)
+    permute.set_direction(d, abs(permute.in_direction(d)));
+  
+  // compute the size of the chunk to output, and its strides etc.
+  for (int i = 0; i < data->rank; ++i) {
+    direction d = data->ds[i];
+    int isd = isS.in_direction(d), ied = ieS.in_direction(d);
+    count[i] = abs(ied - isd) / 2 + 1;
+    if (ied < isd) offset[permute.in_direction(d)] = count[i] - 1;
+  }
+  for (int i = 0; i < data->rank; ++i) {
+    direction d = data->ds[i];
+    int j = permute.in_direction(d);
+    for (int k = i + 1; k < data->rank; ++k) stride[j] *= count[k];
+    offset[j] *= stride[j];
+    if (offset[j]) stride[j] *= -1;
+  }
+  
+  //-----------------------------------------------------------------------//
+  // Compute the function to output, exactly as in fields::integrate.
+  int *off = data->offsets;
+  component *cS = data->cS;
+  complex<double> *fields = data->fields, *ph = data->ph;
+  const component *iecs = data->inveps_cs;
+  const direction *ieds = data->inveps_ds;
+  int ieos[6];
+  const component *imcs = data->invmu_cs;
+  const direction *imds = data->invmu_ds;
+  int imos[6];
+  int num_components=data->components.size();
+  
+  double *slice=0;
+  cdouble *zslice=0;
+  bool complex_data = (data->rfun==0);
+  if (complex_data)
+   zslice = (cdouble *)data->vslice;
+  else
+   slice = (double *)data->vslice;
+
+  for (int i = 0; i < num_components; ++i) {
+    cS[i] = S.transform(data->components[i], -sn);
+    if (cS[i] == Dielectric || cS[i] == Permeability)
+      ph[i] = 1.0;
+    else {
+      fc->gv.yee2cent_offsets(cS[i], off[2*i], off[2*i+1]);
+      ph[i] = shift_phase * S.phase_shift(cS[i], sn);
+    }
+  }
+  for (int k = 0; k < data->ninveps; ++k)
+    fc->gv.yee2cent_offsets(iecs[k], ieos[2*k], ieos[2*k+1]);
+  for (int k = 0; k < data->ninvmu; ++k)
+    fc->gv.yee2cent_offsets(imcs[k], imos[2*k], imos[2*k+1]);
+
+  vec rshift(shift * (0.5*fc->gv.inva));
+  LOOP_OVER_IVECS(fc->gv, is, ie, idx) {
+    IVEC_LOOP_LOC(fc->gv, loc);
+    loc = S.transform(loc, sn) + rshift;
+
+    for (int i = 0; i < num_components; ++i) {
+      if (cS[i] == Dielectric) {
+	double tr = 0.0;
+	for (int k = 0; k < data->ninveps; ++k) {
+	  const realnum *ie = fc->s->chi1inv[iecs[k]][ieds[k]];
+	  if (ie) tr += (ie[idx] + ie[idx+ieos[2*k]] + ie[idx+ieos[1+2*k]]
+			 + ie[idx+ieos[2*k]+ieos[1+2*k]]);
+	  else tr += 4; // default inveps == 1
+	}
+	fields[i] = (4 * data->ninveps) / tr;
+      }
+      else if (cS[i] == Permeability) {
+	double tr = 0.0;
+	for (int k = 0; k < data->ninvmu; ++k) {
+	  const realnum *im = fc->s->chi1inv[imcs[k]][imds[k]];
+	  if (im) tr += (im[idx] + im[idx+imos[2*k]] + im[idx+imos[1+2*k]]
+			 + im[idx+imos[2*k]+imos[1+2*k]]);
+	  else tr += 4; // default invmu == 1
+	}
+	fields[i] = (4 * data->ninvmu) / tr;
+      }
+      else {
+	double f[2];
+	for (int k = 0; k < 2; ++k)
+	  if (fc->f[cS[i]][k])
+	    f[k] = 0.25 * (fc->f[cS[i]][k][idx]
+			   + fc->f[cS[i]][k][idx+off[2*i]]
+			   + fc->f[cS[i]][k][idx+off[2*i+1]]
+			   + fc->f[cS[i]][k][idx+off[2*i]+off[2*i+1]]);
+	  else
+	    f[k] = 0;
+	fields[i] = complex<double>(f[0], f[1]) * ph[i];
+      }
+    }
+    int idx2 = ((((offset[0] + offset[1] + offset[2])
+                   + loop_i1 * stride[0])
+                   + loop_i2 * stride[1])
+                   + loop_i3 * stride[2]);
+
+    if (complex_data)
+     zslice[idx2] = data->fun(fields, loc, data->fun_data);
+    else
+     slice[idx2]  = data->rfun(fields, loc, data->fun_data);
+
+  };
+
+}
+
+/***************************************************************/
+/* given a volume, fill in the dims[] and directions[] arrays  */
+/* describing the array slice needed to store field data for   */
+/* all grid points in the volume.                              */
+/*                                                             */
+/* return value is rank of array slice.                        */
+/*                                                             */
+/* if caller_data is non-NULL, it should point to a            */
+/* caller-allocated array_slice_data structure which will be   */
+/* initialized appopriately for subsequent use in              */
+/* get_array_slice.                                            */
+/***************************************************************/
+int fields::get_array_slice_dimensions(const volume &where, int dims[3], void *caller_data)
+{
+  am_now_working_on(ArraySlice);
+
+  // use a local data structure if the caller didn't provide one
+  array_slice_data local_data;
+  array_slice_data *data=(array_slice_data *)caller_data;
+  if (data==0)
+   data=&local_data;
+
+  data->min_corner = gv.round_vec(where.get_max_corner()) + one_ivec(gv.dim);
+  data->max_corner = gv.round_vec(where.get_min_corner()) - one_ivec(gv.dim);
+  data->num_chunks = 0;
+
+  loop_in_chunks(get_array_slice_dimensions_chunkloop,
+                 (void *) data, where, Centered, true, true);
+
+  data->max_corner = max_to_all(data->max_corner);
+  data->min_corner = -max_to_all(-data->min_corner); // i.e., min_to_all
+  data->num_chunks = sum_to_all(data->num_chunks);
+  if (data->num_chunks == 0 || !(data->min_corner <= data->max_corner))
+    return 0; // no data to write;
+
+  int rank=0, slice_size=1;
+  LOOP_OVER_DIRECTIONS(gv.dim, d) {
+    if (rank >= 3) abort("too many dimensions in array_slice");
+    int n = (data->max_corner.in_direction(d)
+	     - data->min_corner.in_direction(d)) / 2 + 1;
+    if (n > 1) {
+      data->ds[rank] = d;
+      dims[rank++] = n;
+      slice_size *= n;
+    }
+  }
+  data->rank=rank;
+  data->slice_size=slice_size;
+  finished_working();
+
+  return rank;
+}
+
+/***************************************************************/
+/* precisely one of fun, rfun should be non-NULL               */
+/***************************************************************/
+void *fields::do_get_array_slice(const volume &where,
+                                 std::vector<component> components,
+                                 field_function fun,
+                                 field_rfunction rfun,
+                                 void *fun_data,
+                                 void *vslice) {
+
+  am_now_working_on(ArraySlice);
+
+  /***************************************************************/
+  /* call get_array_slice_dimensions to get slice dimensions and */
+  /* partially initialze an array_slice_data struct              */
+  /***************************************************************/
+  int dims[3];
+  array_slice_data data;
+  int rank=get_array_slice_dimensions(where, dims, &data);
+  int slice_size=data.slice_size;
+  if (rank==0 || slice_size==0) return 0; // no data to write
+
+  bool complex_data = (rfun==0);
+  cdouble *zslice;
+  double *slice;
+  if (vslice==0)
+   { if (complex_data)
+      { zslice = new cdouble[slice_size];
+        vslice = (void *)zslice;
+      } 
+     else
+      { slice  = new double[slice_size];
+        vslice = (void *)slice;
+      };
+   };
+   
+  data.vslice     = vslice;
+  data.fun        = fun;
+  data.rfun       = rfun;
+  data.fun_data   = fun_data;
+  data.components = components;
+
+  int num_components = components.size();
+
+  data.cS      = new component[num_components];
+  data.ph      = new cdouble[num_components];
+  data.fields  = new cdouble[num_components];
+  
+  data.offsets = new int[2 * num_components];
+  for (int i = 0; i < 2 * num_components; ++i)
+    data.offsets[i] = 0;
+
+  /* compute inverse-epsilon directions for computing Dielectric fields */
+  data.ninveps = 0;
+  bool needs_dielectric = false;
+  for (int i = 0; i < num_components; ++i)
+    if (components[i] == Dielectric) { needs_dielectric = true; break; }
+  if (needs_dielectric) 
+    FOR_ELECTRIC_COMPONENTS(c) if (gv.has_field(c)) {
+      if (data.ninveps == 3) abort("more than 3 field components??");
+      data.inveps_cs[data.ninveps] = c;
+      data.inveps_ds[data.ninveps] = component_direction(c);
+      ++data.ninveps;
+    }
+  
+  /* compute inverse-mu directions for computing Permeability fields */
+  data.ninvmu = 0;
+  bool needs_permeability = false;
+  for (int i = 0; i < num_components; ++i)
+    if (components[i] == Permeability) { needs_permeability = true; break; }
+  if (needs_permeability) 
+    FOR_MAGNETIC_COMPONENTS(c) if (gv.has_field(c)) {
+      if (data.ninvmu == 3) abort("more than 3 field components??");
+      data.invmu_cs[data.ninvmu] = c;
+      data.invmu_ds[data.ninvmu] = component_direction(c);
+      ++data.ninvmu;
+    }
+ 
+  loop_in_chunks(get_array_slice_chunkloop, (void *) &data,
+		 where, Centered, true, true);
+
+  /***************************************************************/
+  /* repeatedly call sum_to_all to consolidate full array slice  */
+  /* on all cores                                                */
+  /***************************************************************/
+#define BUFSIZE 1<<16 // use 64k buffer
+  if (complex_data)
+   { cdouble *buffer = new cdouble[BUFSIZE];
+     cdouble *slice = (cdouble *)vslice;
+     int offset=0, remaining=slice_size;
+     while(remaining!=0)
+      { 
+        int size = (remaining > BUFSIZE ? BUFSIZE : remaining);
+        sum_to_all(slice + offset, buffer, size);
+        memcpy(slice+offset, buffer, size*sizeof(cdouble));
+        remaining-=size;
+        offset+=size;
+      };
+     delete[] buffer;
+   }
+  else
+   { double *buffer = new double[BUFSIZE];
+     double *slice = (double *)vslice;
+     int offset=0, remaining=slice_size;
+     while(remaining!=0)
+      { int size = (remaining > BUFSIZE ? BUFSIZE : remaining);
+        sum_to_all(slice + offset, buffer, size);
+        memcpy(slice+offset, buffer, size*sizeof(double));
+        remaining-=size;
+        offset+=size;
+      };
+     delete[] buffer;
+   };
+
+  delete[] data.offsets;
+  delete[] data.fields;
+  delete[] data.ph;
+  delete[] data.cS;
+  finished_working();
+
+  return vslice;
+}
+
+/***************************************************************/
+/* entry points to get_array_slice                             */
+/***************************************************************/
+double *fields::get_array_slice(const volume &where,
+                                std::vector<component> components,
+                                field_rfunction rfun, void *fun_data,
+                                double *slice)
+{
+  return (double *)do_get_array_slice(where, components,
+                                      0, rfun, fun_data,
+                                      (void *)slice);
+} 
+
+cdouble *fields::get_complex_array_slice(const volume &where,
+                                         std::vector<component> components,
+                                         field_function fun, void *fun_data,
+                                         cdouble *slice)
+{
+  return (cdouble *)do_get_array_slice(where, components,
+                                       fun, 0, fun_data,
+                                       (void *)slice);
+}
+
+double *fields::get_array_slice(const volume &where, component c,
+                                double *slice, int slice_length)
+{
+  (void) slice_length;
+  std::vector<component> components(1);
+  components[0]=c;
+  return (double *)do_get_array_slice(where, components,
+                                      0, default_field_rfunc, 0,
+                                      (void *)slice);
+}
+
+double *fields::get_array_slice(const volume &where,
+                                derived_component c,
+                                double *slice, int slice_length)
+{
+  (void) slice_length;
+  int nfields;
+  component carray[12];
+  field_rfunction rfun = derived_component_func(c, gv, nfields, carray);
+  std::vector<component> cs(carray, carray+nfields);
+  return (double *)do_get_array_slice(where, cs,
+                                      0, rfun, &nfields,
+                                      (void *)slice);
+}
+
+cdouble *fields::get_complex_array_slice(const volume &where, component c,
+                                         cdouble *slice, int slice_length)
+{
+  (void) slice_length;
+  std::vector<component> components(1);
+  components[0]=c;
+  return (cdouble *)do_get_array_slice(where, components,
+                                       default_field_func, 0, 0,
+                                       (void *)slice);
+}
+
+} // namespace meep
diff --git a/src/h5file.cpp b/src/h5file.cpp
index 4c3ac1dcd..5f7b89b9d 100644
--- a/src/h5file.cpp
+++ b/src/h5file.cpp
@@ -309,7 +309,7 @@ void h5file::read_size(const char *dataname, int *rank, int *dims, int maxrank)
 #define REALNUM_H5T (sizeof(realnum) == sizeof(double) ? H5T_NATIVE_DOUBLE : H5T_NATIVE_FLOAT)
 
 realnum *h5file::read(const char *dataname,
-		     int *rank, int *dims, int maxrank)
+	 	     int *rank, int *dims, int maxrank)
 {
 #ifdef HAVE_HDF5
   realnum *data = 0;
diff --git a/src/meep.hpp b/src/meep.hpp
index 9c9d4da16..6774568b1 100644
--- a/src/meep.hpp
+++ b/src/meep.hpp
@@ -23,6 +23,8 @@
 #include "meep/vec.hpp"
 #include "meep/mympi.hpp"
 
+#include <vector>
+
 namespace meep {
 
 /* We use the type realnum for large arrays, e.g. the fields.
@@ -1165,7 +1167,8 @@ class fields_chunk {
 
 enum boundary_condition { Periodic=0, Metallic, Magnetic, None };
 enum time_sink { Connecting, Stepping, Boundaries, MpiTime,
-                 FieldOutput, FourierTransforming, Other };
+                 FieldOutput, FourierTransforming, 
+                 ArraySlice, Other };
 
 typedef void (*field_chunkloop)(fields_chunk *fc, int ichunk, component cgrid,
 				ivec is, ivec ie,
@@ -1285,6 +1288,50 @@ class fields {
   const char *h5file_name(const char *name,
 			  const char *prefix = NULL, bool timestamp = false);
 
+  // array_slice.cpp methods
+
+  // given a subvolume, compute the dimensions of the array slice
+  // needed to store field data for that subvolume.
+  // the data parameter is used internally in get_array_slice
+  // and should be ignored by external callers.
+  int get_array_slice_dimensions(const volume &where, int dims[3], void *data=0);
+
+  // given a subvolume, return a column-major array containing
+  // the given function of the field components in that subvolume
+  // if slice is non-null, it must be a user-allocated buffer
+  // of the correct size.
+  // otherwise, a new buffer is allocated and returned; it
+  // must eventually be caller-deallocated via delete[].
+  double *get_array_slice(const volume &where,
+                          std::vector<component> components,
+                          field_rfunction rfun, void *fun_data,
+                          double *slice=0);
+
+  std::complex<double> *get_complex_array_slice(const volume &where,
+                                   std::vector<component> components,
+                                   field_function fun,
+                                   void *fun_data,
+                                   std::complex<double> *slice=0);
+
+  // alternative entry points for when you have no field
+  // function, i.e. you want just a single component or 
+  // derived component.) (The slice_length parameter is only
+  // present to facilitate SWIG-generated python wrappers;
+  // it is ignored by the C code).
+  double *get_array_slice(const volume &where, component c, double *slice=0, int slice_length=0);
+  double *get_array_slice(const volume &where, derived_component c, double *slice=0, int slice_length=0);
+  std::complex<double> *get_complex_array_slice(const volume &where,
+                                                component c,
+                                                std::complex<double> *slice=0, int slice_length=0);
+
+  // master routine for all above entry points
+  void *do_get_array_slice(const volume &where,
+                           std::vector<component> components,
+                           field_function fun,
+                           field_rfunction rfun,
+                           void *fun_data,
+                           void *vslice);
+
   // step.cpp methods:
   double last_step_output_wall_time;
   int last_step_output_t;
diff --git a/tests/Makefile.am b/tests/Makefile.am
index e03ce8fb2..cbd1bcd74 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -74,7 +74,7 @@ h5test_LDADD = $(LIBMEEP)
 pml_SOURCES = pml.cpp
 pml_LDADD = $(LIBMEEP)
 
-TESTS = aniso_disp bench bragg_transmission convergence_cyl_waveguide cylindrical flux harmonics integrate known_results near2far one_dimensional physical stress_tensor symmetry three_d two_dimensional 2D_convergence h5test pml
+TESTS = aniso_disp bench bragg_transmission convergence_cyl_waveguide cylindrical flux harmonics integrate known_results near2far one_dimensional physical stress_tensor symmetry three_d two_dimensional 2D_convergence h5test pml 
 
 LOG_COMPILER = $(RUNCODE)
 
diff --git a/tests/h5test.cpp b/tests/h5test.cpp
index 247f38d5b..00bb79c4a 100644
--- a/tests/h5test.cpp
+++ b/tests/h5test.cpp
@@ -185,7 +185,7 @@ bool check_2d(double eps(const vec &), double a, int splitting, symfunc Sf,
     delete[] h5data;
   }
 
-  file->remove();
+  //file->remove();
   delete file;
 
   master_printf("Passed %s (%g..%g), err=%g\n", name,
@@ -300,7 +300,7 @@ bool check_3d(double eps(const vec &), double a, int splitting, symfunc Sf,
     delete[] h5data;
   }
 
-  file->remove();
+  //file->remove();
   delete file;
 
   master_printf("Passed %s (%g..%g), err=%g\n", name,
@@ -384,7 +384,7 @@ bool check_2d_monitor(double eps(const vec &),
 
   delete[] mon;
 
-  file->remove();
+  //file->remove();
   delete file;
 
   master_printf("Passed %s (%g..%g), err=%g\n", name,