diff --git a/visualisation/Visualise5D.cxx b/visualisation/Visualise5D.cxx
new file mode 100755
index 00000000..80825cbb
--- /dev/null
+++ b/visualisation/Visualise5D.cxx
@@ -0,0 +1,729 @@
+// Derived from VTK/Examples/Cxx/Medical2.cxx
+// The example reads a volume dataset, extracts two isosurfaces that
+// represent the skin and bone, and then displays them.
+//
+// Modified heavily by Peter Boyle to display lattice field theory data as movies and compare multiple files
+
+#include <vtkActor.h>
+#include <vtkCamera.h>
+#include <vtkMetaImageReader.h>
+#include <vtkNamedColors.h>
+#include <vtkNew.h>
+#include <vtkOutlineFilter.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkProperty.h>
+#include <vtkRenderWindow.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkRenderer.h>
+#include <vtkStripper.h>
+#include <vtkImageData.h>
+#include <vtkVersion.h>
+#include <vtkCallbackCommand.h>
+#include <vtkTextActor.h>
+#include <vtkTextProperty.h>
+
+#define MPEG
+#ifdef MPEG
+#include <vtkFFMPEGWriter.h>
+#endif
+
+#include <vtkProperty2D.h>
+#include <vtkSliderWidget.h>
+#include <vtkSliderRepresentation2D.h>
+#include <vtkWindowToImageFilter.h>
+
+#include <array>
+#include <string>
+
+#include <Grid/Grid.h>
+
+#define USE_FLYING_EDGES
+#ifdef USE_FLYING_EDGES
+#include <vtkFlyingEdges3D.h>
+typedef vtkFlyingEdges3D isosurface;
+#else
+#include <vtkMarchingCubes.h>
+typedef vtkMarchingCubes isosurface;
+#endif
+
+int mpeg = 0 ;
+int framerate = 10;
+
+template <class T> void readFile(T& out, std::string const fname){
+  Grid::emptyUserRecord record;
+  Grid::ScidacReader RD;
+  RD.open(fname);
+  RD.readScidacFieldRecord(out,record);
+  RD.close();
+}
+using namespace Grid;
+
+class FrameUpdater : public vtkCallbackCommand
+{
+public:
+
+  FrameUpdater() {
+    ffile=0;
+    TimerCount = 0;
+    xoff       = 0;
+    t          = 0;
+    imageData = nullptr;
+    timerId = 0;
+    maxCount = -1;
+    old_file=-1;
+  }
+  
+  static FrameUpdater* New()
+  {
+    FrameUpdater* cb = new FrameUpdater;
+    cb->TimerCount = 0;
+    return cb;
+  }
+
+  //
+  // Must map a x,y,z + frame index into
+  // i)  a d-dimensional site Coordinate
+  // ii) a file name
+  // Need a:
+  //     loop_ranges
+  //     sum_ranges
+  //     loop_vol  -- map loop_idx -> loop_coor
+  //     sum_vol   -- map sum_idx -> sum_coor with Lexicographic
+  //
+  /*
+   * Just set this up
+   */
+  int old_file ; // Cache, avoid reread
+  
+  Coordinate          latt;
+  Coordinate          xyz_dims    ; // List lattice dimensions corresponding to xyz_dims displayed 
+  Coordinate          xyz_ranges  ; // 3-vector
+  Coordinate          g_xyz_ranges; // Nd-vector
+  uint64_t            xyz_vol     ;
+  
+  Coordinate          loop_dims;    // List lattice dimensions put into movie time
+  Coordinate          loop_ranges;  // movie time ranges
+  uint64_t            loop_vol;
+
+  Coordinate          sum_dims;      // List lattice dimensions summed
+  Coordinate          sum_ranges;    // summation ranges
+  uint64_t            sum_vol;
+
+  Coordinate          slice_dims;      // List slice dimensions 
+
+  Coordinate          Slice;
+  
+  std::vector<std::string> files;        // file list that is looped over
+
+  int Nd;
+  GridBase *grid;
+  Grid::LatticeComplexD *grid_data;
+
+  void SetGrid(GridBase *_grid)
+  {
+    grid = _grid;
+    Nd=grid->Nd();
+    latt = grid->GlobalDimensions();
+    grid_data = new Grid::LatticeComplexD(grid);
+  }
+
+  void SetFiles(std::vector<std::string> list)      { files = list; old_file = -1; }
+  void SetSlice(Coordinate _Slice)             { Slice = _Slice;} // Offset / skew for lattice coords
+  void SetSumDimensions (Coordinate _SumDims ) {
+    sum_ranges=Coordinate(Nd);
+    sum_dims = _SumDims; // 1 hot for dimensions summed
+    sum_vol = 1;
+    for(int d=0;d<sum_dims.size();d++){
+      if ( sum_dims[d] == 1 ) sum_ranges[d] = latt[d];
+      else                    sum_ranges[d] = 1;
+      sum_vol*=sum_ranges[d];
+    }
+  }
+  void SetLoopDimensions(Coordinate _LoopDims) {
+    loop_ranges=Coordinate(Nd);
+    loop_dims= _LoopDims;
+    loop_vol = 1;
+    for(int d=0;d<loop_dims.size();d++){
+      if ( loop_dims[d] == 1 ) loop_ranges[d] = latt[d];
+      else                     loop_ranges[d] = 1;
+      loop_vol*=loop_ranges[d];
+    }
+  } // 
+  void SetDisplayDimensions(Coordinate _xyz_dims   ) {
+    g_xyz_ranges=Coordinate(Nd);
+    xyz_ranges=Coordinate(3);
+    xyz_dims = _xyz_dims;
+    xyz_vol  = 1;
+    for(int d=0;d<3;d++){
+      xyz_ranges[d] = latt[xyz_dims[d]];
+      xyz_vol *= xyz_ranges[d];
+    }
+    // Find dim extents for grid
+    int dd=0;
+    for(int d=0;d<Nd;d++){
+      g_xyz_ranges[d] = 1;
+      for(int dd=0;dd<3;dd++) {
+	if ( xyz_dims[dd]==d ) {
+	  g_xyz_ranges[d] = latt[d];
+	}
+      }
+    }
+  }
+  void SetSliceDimensions(void) {
+    Coordinate _slice_dims;
+    for ( int d=0;d<Nd;d++){
+      int is_slice = 1;
+      if(g_xyz_ranges[d]>1) is_slice = 0;
+      if(loop_dims[d]) is_slice = 0;
+      if(sum_dims[d] ) is_slice = 0;
+      if(is_slice) _slice_dims.push_back(d);
+    }
+    slice_dims = _slice_dims;
+    std::cout << " Setting Slice Dimensions to "<<slice_dims<<std::endl;
+  }
+  virtual void Execute(vtkObject* caller, unsigned long eventId,void* vtkNotUsed(callData))
+  {
+    const int max=256;
+    char text_string[max];
+
+    auto latt_size = grid->GlobalDimensions();
+    
+    if ( vtkCommand::KeyPressEvent == eventId ) {
+
+      vtkRenderWindowInteractor* iren = static_cast<vtkRenderWindowInteractor*>(caller);
+      std::string key = iren->GetKeySym();
+      std::cout << "Pressed: " << key << std::endl;
+      if (slice_dims.size()>0) {
+	int vert = slice_dims[slice_dims.size()-1];
+	int horz = slice_dims[0];
+	if ( key == "Up" ) {
+	  Slice[vert] = (Slice[vert]+1)%latt[vert];
+	}
+	if ( key == "Down" ) {
+	  Slice[vert] = (Slice[vert]+latt[vert]-1)%latt[vert];
+	}
+	if ( key == "Right" ) {
+	  Slice[horz] = (Slice[horz]+1)%latt[horz];
+	}
+	if ( key == "Left" ) {
+	  Slice[horz] = (Slice[horz]+latt[horz]-1)%latt[horz];
+	}
+      }
+      if ( key == "greater" ) {
+	ffile = (ffile + 1) % files.size();
+      }
+      if ( key == "less" ) {
+	ffile = (ffile - 1 + files.size()) % files.size();
+      }
+      std::cout <<"Slice " <<Slice <<std::endl;
+      std::cout <<"File  " <<ffile <<std::endl;
+    }
+      
+    // Make a new frame for frame index TimerCount
+    if ( vtkCommand::TimerEvent == eventId || vtkCommand::KeyPressEvent == eventId)
+      {
+	int file     = ((this->TimerCount / loop_vol) + ffile )%files.size();
+
+	if ( file != old_file ) {
+	  readFile(*grid_data,files[file]);
+	  old_file = file;
+	}
+
+	RealD max, min, max_abs,min_abs;
+	Coordinate max_site;
+	Coordinate min_site;
+	Coordinate max_abs_site;
+	Coordinate min_abs_site;
+	for(int idx=0;idx<grid->gSites();idx++){
+	  Coordinate site;
+	  Lexicographic::CoorFromIndex (site,idx,latt);
+	  RealD val=real(peekSite(*grid_data,site));
+	  if (idx==0){
+	    max = min = val;
+	    max_abs = min_abs = fabs(val);
+	    max_site = site;
+	    min_site = site;
+	    min_abs_site = site;
+	    max_abs_site = site;
+	  } else {
+	    if ( val > max ) {
+	      max=val;
+	      max_site = site;
+	    }
+	    if ( fabs(val) > max_abs ) {
+	      max_abs=fabs(val);
+	      max_abs_site = site;
+	    }
+	    if ( val < min ) {
+	      min=val;
+	      min_site = site;
+	    }	    
+	    if ( fabs(val) < min_abs ) {
+	      min_abs=fabs(val);
+	      min_abs_site = site;
+	    }	    
+	  }
+	}
+	std::cout << " abs_max "<<max_abs<<" at " << max_abs_site<<std::endl;
+	std::cout << " abs_min "<<min_abs<<" at " << min_abs_site<<std::endl;
+	std::cout << " max "<<max<<" at " << max_site<<std::endl;
+	std::cout << " min "<<min<<" at " << min_site<<std::endl;
+	// Looped dimensions, map index to coordinate
+	int loop_idx = this->TimerCount % loop_vol;
+	Coordinate loop_coor;
+	Lexicographic::CoorFromIndex (loop_coor,loop_idx,loop_ranges);
+	
+	// Loop over xyz sites
+	Coordinate xyz_coor(3);
+	Coordinate g_xyz_coor(Nd);
+	Coordinate sum_coor(Nd);
+	for(uint64_t xyz = 0 ; xyz< xyz_vol; xyz++){
+	  Lexicographic::CoorFromIndex (xyz_coor,xyz,xyz_ranges);
+	  Lexicographic::CoorFromIndex (g_xyz_coor,xyz,g_xyz_ranges);
+	  
+	  RealD sum_value = 0.0;
+	  for(uint64_t sum_idx = 0 ; sum_idx< sum_vol; sum_idx++){
+	    Lexicographic::CoorFromIndex (sum_coor,sum_idx,sum_ranges);
+	    Coordinate site(Nd);
+	    for(int d=0;d<Nd;d++){
+	      site[d] = (sum_coor[d] + loop_coor[d] + g_xyz_coor[d] + Slice[d])%latt[d];
+	    }
+	    sum_value+= real(peekSite(*grid_data,site));
+	    if(xyz==0) std::cout << "sum "<<sum_idx<<" "<<sum_value<<std::endl;
+	  }
+	  imageData->SetScalarComponentFromDouble(xyz_coor[0],xyz_coor[1],xyz_coor[2],0,sum_value);
+	}
+	imageData->Modified();
+
+	std::stringstream ss;
+	ss<< files[file] <<"\nSlice "<<Slice << "\nLoop  " <<loop_coor<<"\nSummed "<<sum_dims;
+
+	text->SetInput(ss.str().c_str());
+	
+
+	vtkRenderWindowInteractor* iren = dynamic_cast<vtkRenderWindowInteractor*>(caller);
+	iren->GetRenderWindow()->Render();
+	  
+      }
+
+    if ( vtkCommand::TimerEvent == eventId ) {
+      ++this->TimerCount;
+      std::cout << " This was a timer event count "<<this->TimerCount << std::endl;
+    }
+  
+    
+    if (this->TimerCount >= this->maxCount) {
+      vtkRenderWindowInteractor* iren = dynamic_cast<vtkRenderWindowInteractor*>(caller);
+      if (this->timerId > -1)
+	{
+        iren->DestroyTimer(this->timerId);
+	}
+    }
+  }
+  
+
+private:
+  int TimerCount;
+  int ffile;
+  int xoff;
+  int t;
+public:
+  vtkImageData* imageData = nullptr;
+  vtkTextActor* text = nullptr;
+  vtkFFMPEGWriter *writer = nullptr;
+  int timerId ;
+  int maxCount ;
+  double rms;
+  
+  isosurface * posExtractor;
+  isosurface * negExtractor;
+};
+class SliderCallback : public vtkCommand
+{
+public:
+    static SliderCallback* New()
+    {
+        return new SliderCallback;
+    }
+    virtual void Execute(vtkObject* caller, unsigned long eventId, void* callData)
+    {
+        vtkSliderWidget *sliderWidget = vtkSliderWidget::SafeDownCast(caller);
+        if (sliderWidget)
+        {
+	  contour = ((vtkSliderRepresentation *)sliderWidget->GetRepresentation())->GetValue();
+        }
+	fu->posExtractor->SetValue(0,  SliderCallback::contour*fu->rms);
+	fu->negExtractor->SetValue(0, -SliderCallback::contour*fu->rms);
+	fu->posExtractor->Modified();
+	fu->negExtractor->Modified();
+    }    
+public:
+  static double contour;
+  FrameUpdater * fu;
+};
+
+FrameUpdater * KBfu;
+void KeypressCallbackFunction(vtkObject* caller, long unsigned int eventId,
+                              void* clientData, void* callData)
+{
+  std::cout << "Keypress callback" << std::endl;
+  vtkRenderWindowInteractor* iren = static_cast<vtkRenderWindowInteractor*>(caller);
+  std::cout << "Pressed: " << iren->GetKeySym() << std::endl;
+  //  imageData->Modified();
+}
+
+double SliderCallback::contour;
+
+int main(int argc, char* argv[])
+{
+  using namespace Grid;
+
+  Grid_init(&argc, &argv);
+  GridLogLayout();
+
+  auto latt_size   = GridDefaultLatt();
+  auto simd_layout = GridDefaultSimd(latt_size.size(), vComplex::Nsimd());
+  auto mpi_layout  = GridDefaultMpi();
+  GridCartesian    Grid(latt_size, simd_layout, mpi_layout);
+
+  double default_contour = 1.0;
+  std::string arg;
+ 
+  std::cout << argc << " command Line arguments "<<std::endl;
+  for(int c=0;c<argc;c++) {
+    std::cout << " - "<<argv[c]<<std::endl;
+  }
+
+  std::vector<std::string> file_list({
+      "file1",
+      "file2",
+      "file3",
+      "file4",
+      "file5",
+      "file6",
+      "file7",
+      "file8"
+    });
+
+  if( GridCmdOptionExists(argv,argv+argc,"--files") ){
+    arg=GridCmdOptionPayload(argv,argv+argc,"--files");
+    GridCmdOptionCSL(arg, file_list);
+  }
+
+#ifdef MPEG
+  if( GridCmdOptionExists(argv,argv+argc,"--mpeg") ){
+    mpeg = 1;
+  }
+#endif
+
+  if( GridCmdOptionExists(argv,argv+argc,"--fps") ){
+    arg=GridCmdOptionPayload(argv,argv+argc,"--fps");
+    GridCmdOptionInt(arg,framerate);
+  }
+  if( GridCmdOptionExists(argv,argv+argc,"--isosurface") ){
+    arg=GridCmdOptionPayload(argv,argv+argc,"--isosurface");
+    GridCmdOptionFloat(arg,default_contour);
+  }
+  for(int c=0;c<file_list.size();c++) {
+    std::cout << " file: "<<file_list[c]<<std::endl;
+  }
+
+  int NoTime = 0;
+  int Nd; Nd = Grid.Nd();
+  Coordinate    Slice(Nd,0);
+  Coordinate  SumDims(Nd,0);
+  Coordinate LoopDims(Nd,0);
+  Coordinate  XYZDims({0,1,2});
+  if( GridCmdOptionExists(argv,argv+argc,"--slice") ){
+    arg=GridCmdOptionPayload(argv,argv+argc,"--slice");
+    GridCmdOptionIntVector(arg,Slice);
+  }
+  if( GridCmdOptionExists(argv,argv+argc,"--sum") ){
+    arg=GridCmdOptionPayload(argv,argv+argc,"--sum");
+    GridCmdOptionIntVector(arg,SumDims);
+  }
+  if( GridCmdOptionExists(argv,argv+argc,"--loop") ){
+    arg=GridCmdOptionPayload(argv,argv+argc,"--loop");
+    GridCmdOptionIntVector(arg,LoopDims);
+  }
+  if( GridCmdOptionExists(argv,argv+argc,"--xyz") ){
+    arg=GridCmdOptionPayload(argv,argv+argc,"--xyz");
+    GridCmdOptionIntVector(arg,XYZDims);
+    std::cout << "xyz : "<<XYZDims<<std::endl;
+  }
+  if( GridCmdOptionExists(argv,argv+argc,"--notime") ){
+    NoTime = 1;
+    std::cout << "Suppressing time loop"<<std::endl;
+  }
+
+  // Common things:
+  vtkNew<vtkNamedColors> colors;
+  std::array<unsigned char, 4> posColor{{240, 184, 160, 255}};  colors->SetColor("posColor", posColor.data());
+  std::array<unsigned char, 4> bkg{{51, 77, 102, 255}};         colors->SetColor("BkgColor", bkg.data());
+
+  // Create the renderer, the render window, and the interactor. The renderer
+  // draws into the render window, the interactor enables mouse- and
+  // keyboard-based interaction with the data within the render window.
+  //
+  vtkNew<vtkRenderWindow> renWin;
+  vtkNew<vtkRenderWindowInteractor> iren;
+  iren->SetRenderWindow(renWin);
+
+  
+  //  std::vector<LatticeComplexD> data(file_list.size(),&Grid);
+  //  FieldMetaData header;
+
+  
+
+  int frameCount = file_list.size();
+  for(int d=0;d<Grid.Nd();d++) {
+    if ( LoopDims[d] ) frameCount*= latt_size[d];
+
+  }
+  // It is convenient to create an initial view of the data. The FocalPoint
+  // and Position form a vector direction. Later on (ResetCamera() method)
+  // this vector is used to position the camera to look at the data in
+  // this direction.
+  vtkNew<vtkCamera> aCamera;
+  aCamera->SetViewUp(0, 0, -1);
+  aCamera->SetPosition(0, -1000, 0);
+  aCamera->SetFocalPoint(0, 0, 0);
+  aCamera->ComputeViewPlaneNormal();
+  aCamera->Azimuth(30.0);
+  aCamera->Elevation(30.0);
+
+    
+  vtkNew<vtkRenderer> aRenderer;
+  renWin->AddRenderer(aRenderer);
+    
+  double vol = Grid.gSites();
+  std::cout << "Reading "<<file_list[0]<<std::endl;
+  double nrm, nrmbar,rms, contour;
+  {
+    LatticeComplexD data(&Grid);
+    readFile(data,file_list[0]);
+    nrm    = norm2(data);
+  }
+  nrmbar = nrm/vol;
+  rms    = sqrt(nrmbar);
+  
+  contour = default_contour * rms; // default to 1 x RMS
+
+  // The following reader is used to read a series of 2D slices (images)
+  // that compose the volume. The slice dimensions are set, and the
+  // pixel spacing. The data Endianness must also be specified. The reader
+  // uses the FilePrefix in combination with the slice number to construct
+  // filenames using the format FilePrefix.%d. (In this case the FilePrefix
+  // is the root name of the file: quarter.)
+  vtkNew<vtkImageData> imageData;
+  imageData->SetDimensions(latt_size[0],latt_size[1],latt_size[2]);
+  imageData->AllocateScalars(VTK_DOUBLE, 1);
+  for(int xx=0;xx<latt_size[0];xx++){
+    for(int yy=0;yy<latt_size[1];yy++){
+      for(int zz=0;zz<latt_size[2];zz++){
+	Coordinate site({xx,yy,zz,0});
+	RealD value = 0;
+	imageData->SetScalarComponentFromDouble(xx,yy,zz,0,value);
+   }}}
+
+  
+  vtkNew<isosurface> posExtractor;
+  posExtractor->SetInputData(imageData);
+  posExtractor->SetValue(0, contour);
+  
+  vtkNew<vtkStripper> posStripper;
+  posStripper->SetInputConnection(posExtractor->GetOutputPort());
+
+  vtkNew<vtkPolyDataMapper> posMapper;
+  posMapper->SetInputConnection(posStripper->GetOutputPort());
+  posMapper->ScalarVisibilityOff();
+
+  vtkNew<vtkActor> pos;
+  pos->SetMapper(posMapper);
+  pos->GetProperty()->SetDiffuseColor(colors->GetColor3d("posColor").GetData());
+  pos->GetProperty()->SetSpecular(0.3);
+  pos->GetProperty()->SetSpecularPower(20);
+  pos->GetProperty()->SetOpacity(0.5);
+
+  // An isosurface, or contour value is set
+  // The triangle stripper is used to create triangle strips from the
+  // isosurface; these render much faster on may systems.
+  vtkNew<isosurface> negExtractor;
+  negExtractor->SetInputData(imageData);
+  negExtractor->SetValue(0, -contour);
+
+  vtkNew<vtkStripper> negStripper;
+  negStripper->SetInputConnection(negExtractor->GetOutputPort());
+
+  vtkNew<vtkPolyDataMapper> negMapper;
+  negMapper->SetInputConnection(negStripper->GetOutputPort());
+  negMapper->ScalarVisibilityOff();
+    
+  vtkNew<vtkActor> neg;
+  neg->SetMapper(negMapper);
+  neg->GetProperty()->SetDiffuseColor(colors->GetColor3d("Ivory").GetData());
+
+  // An outline provides context around the data.
+  vtkNew<vtkOutlineFilter> outlineData;
+  outlineData->SetInputData(imageData);
+
+  vtkNew<vtkPolyDataMapper> mapOutline;
+  mapOutline->SetInputConnection(outlineData->GetOutputPort());
+  
+  vtkNew<vtkActor> outline;
+  outline->SetMapper(mapOutline);
+  outline->GetProperty()->SetColor(colors->GetColor3d("Black").GetData());
+
+  vtkNew<vtkTextActor> Text;
+  //  Text->SetInput(file_list[f].c_str());
+  Text->SetPosition2(0,0);
+  Text->GetTextProperty()->SetFontSize(24);
+  Text->GetTextProperty()->SetColor(colors->GetColor3d("Gold").GetData());
+
+  vtkNew<vtkTextActor> TextT;
+  TextT->SetInput("T=0");
+  TextT->SetPosition(0,.7*1025);
+  TextT->GetTextProperty()->SetFontSize(24);
+  TextT->GetTextProperty()->SetColor(colors->GetColor3d("Gold").GetData());
+    
+  
+  // Actors are added to the renderer. An initial camera view is created.
+  // The Dolly() method moves the camera towards the FocalPoint,
+  // thereby enlarging the image.
+  //    aRenderer->AddActor(Text);
+  aRenderer->AddActor(TextT);
+  aRenderer->AddActor(outline);
+  aRenderer->AddActor(pos);
+  aRenderer->AddActor(neg);
+
+  // Sign up to receive TimerEvent
+  vtkNew<FrameUpdater> fu;
+
+  fu->SetGrid(&Grid);
+  fu->SetFiles(file_list);
+  fu->SetSlice(Slice);
+  fu->SetSumDimensions (SumDims);
+  fu->SetLoopDimensions(LoopDims);
+  fu->SetDisplayDimensions(XYZDims);
+  fu->SetSliceDimensions();
+  
+
+  fu->imageData = imageData;
+  //    fu->grid_data = &data[f];
+  fu->text      = TextT;
+  fu->maxCount = frameCount;
+  fu->posExtractor = posExtractor;
+  fu->negExtractor = negExtractor;
+  fu->rms = rms;
+      
+  iren->AddObserver(vtkCommand::TimerEvent, fu);
+  iren->AddObserver(vtkCommand::KeyPressEvent, fu);
+
+  aRenderer->SetActiveCamera(aCamera);
+  aRenderer->ResetCamera();
+  aRenderer->SetBackground(colors->GetColor3d("BkgColor").GetData());
+  aCamera->Dolly(1.0);
+
+  //    double nf = file_list.size();
+  //    std::cout << " Adding renderer " <<f<<" of "<<nf<<std::endl;
+  aRenderer->SetViewport(0.0, 0.0,1.0 , 1.0);
+
+  // Note that when camera movement occurs (as it does in the Dolly()
+  // method), the clipping planes often need adjusting. Clipping planes
+  // consist of two planes: near and far along the view direction. The
+  // near plane clips out objects in front of the plane; the far plane
+  // clips out objects behind the plane. This way only what is drawn
+  // between the planes is actually rendered.
+  aRenderer->ResetCameraClippingRange();
+    
+
+  // Set a background color for the renderer and set the size of the
+  // render window (expressed in pixels).
+  // Initialize the event loop and then start it.
+  renWin->SetSize(1024, 1024);
+  renWin->SetWindowName("FieldDensity");
+  renWin->Render();
+
+  // Take a pointer to the FrameUpdater for keypress mgt.
+  //  KBfu = fu;
+  //  vtkNew<vtkCallbackCommand> keypressCallback;
+  //  keypressCallback->SetCallback(KeypressCallbackFunction);
+  //  iren->AddObserver(vtkCommand::KeyPressEvent,keypressCallback);
+
+  iren->Initialize();
+
+  if ( mpeg ) {
+#ifdef MPEG
+    vtkWindowToImageFilter *imageFilter = vtkWindowToImageFilter::New();
+    imageFilter->SetInput( renWin );
+    imageFilter->SetInputBufferTypeToRGB();
+    
+    vtkFFMPEGWriter *writer = vtkFFMPEGWriter::New();
+    writer->SetFileName("movie.avi");
+    writer->SetRate(framerate);
+    writer->SetInputConnection(imageFilter->GetOutputPort());
+    writer->Start();
+
+    for(int i=0;i<fu->maxCount;i++){
+	fu->Execute(iren,vtkCommand::TimerEvent,nullptr);
+	imageFilter->Modified();
+	writer->Write();
+    }
+    writer->End();
+    writer->Delete();
+#else
+    assert(-1 && "MPEG support not compiled");
+#endif
+  } else { 
+  
+    // Add control of contour threshold
+    // Create a slider widget
+    vtkSmartPointer<vtkSliderRepresentation2D> sliderRep = vtkSmartPointer<vtkSliderRepresentation2D>::New();
+    sliderRep->SetMinimumValue(0.0);
+    sliderRep->SetMaximumValue(10.0);
+    sliderRep->SetValue(1.0);
+    sliderRep->SetTitleText("Fraction RMS");
+    // Set color properties:
+
+    // Change the color of the knob that slides
+    //  sliderRep->GetSliderProperty()->SetColor(colors->GetColor3d("Green").GetData());
+    sliderRep->GetTitleProperty()->SetColor(colors->GetColor3d("AliceBlue").GetData());
+    sliderRep->GetLabelProperty()->SetColor(colors->GetColor3d("AliceBlue").GetData());
+    sliderRep->GetSelectedProperty()->SetColor(colors->GetColor3d("DeepPink").GetData());
+
+    // Change the color of the bar
+    sliderRep->GetTubeProperty()->SetColor(colors->GetColor3d("MistyRose").GetData());
+    sliderRep->GetCapProperty()->SetColor(colors->GetColor3d("Yellow").GetData());
+    sliderRep->SetSliderLength(0.05);
+    sliderRep->SetSliderWidth(0.025);
+    sliderRep->SetEndCapLength(0.02);
+
+    sliderRep->GetPoint1Coordinate()->SetCoordinateSystemToNormalizedDisplay();
+    sliderRep->GetPoint1Coordinate()->SetValue(0.1, 0.1);
+    sliderRep->GetPoint2Coordinate()->SetCoordinateSystemToNormalizedDisplay();
+    sliderRep->GetPoint2Coordinate()->SetValue(0.9, 0.1);
+  
+    vtkSmartPointer<vtkSliderWidget> sliderWidget = vtkSmartPointer<vtkSliderWidget>::New();
+    sliderWidget->SetInteractor(iren);
+    sliderWidget->SetRepresentation(sliderRep);
+    sliderWidget->SetAnimationModeToAnimate();
+    sliderWidget->EnabledOn();
+    
+    // Create the slider callback
+    vtkSmartPointer<SliderCallback> slidercallback = vtkSmartPointer<SliderCallback>::New();
+    slidercallback->fu = fu;
+    sliderWidget->AddObserver(vtkCommand::InteractionEvent, slidercallback);
+
+    if ( NoTime==0 ) {
+      int timerId = iren->CreateRepeatingTimer(10000/framerate);
+      std::cout << "timerId "<<timerId<<std::endl;
+    }
+
+    // Start the interaction and timer
+    iren->Start();
+  }
+
+  Grid_finalize();
+
+  return EXIT_SUCCESS;
+}