new file mode 100644
--- /dev/null
+++ b/src/vtkplot.cpp
@@ -0,0 +1,179 @@
+#include <vtkPoints.h>
+#include <vtkPointData.h>
+#include <vtkDelaunay2D.h>
+#include <vtkLookupTable.h>
+
+#include <vtkPolyData.h>
+#include <vtkPolyDataNormals.h>
+#include <vtkPolyDataMapper.h>
+
+#include <vtkWarpScalar.h>
+
+#include <vtkCamera.h>
+#include <vtkActor.h>
+#include <vtkRenderWindow.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkRenderer.h>
+
+#include <vtkInteractorStyleTerrain.h>
+
+#include <vtkSmartPointer.h>
+#include <vtkDoubleArray.h>
+
+#include <vtkGraphicsFactory.h>
+#include <vtkImagingFactory.h>
+#include <vtkWindowToImageFilter.h>
+#include <vtkPNGWriter.h>
+
+#include <string>
+using std::string;
+
+#include "include/linalg.h"
+
+typedef double (*Func3d)(double, double);
+
+void InitOffscreen(bool offscreen);
+vtkSmartPointer<vtkDelaunay2D> TriangulateTerrain(Func3d func);
+
+void PlotPoints(vtkSmartPointer<vtkDelaunay2D> delaunay, 
+                bool offscreen=true,
+                const string& filename="vtkscreenshot.png");
+
+double F(double x, double y)
+{
+  return atan(x*y*sin(y));
+}
+
+int main (int argc, char ** argv)
+{
+	auto delaunay = TriangulateTerrain(&F);	
+  PlotPoints(delaunay,true);
+	return 0;
+}
+
+void InitOffscreen(bool offscreen)
+{
+  // Graphics Factory
+  vtkSmartPointer<vtkGraphicsFactory> graphics_factory
+    = vtkGraphicsFactory::New();
+  graphics_factory->SetOffScreenOnlyMode( offscreen );
+  graphics_factory->SetUseMesaClasses( offscreen );
+
+  // Imaging Factory
+  vtkSmartPointer<vtkImagingFactory>  imaging_factory
+    = vtkImagingFactory::New();
+  imaging_factory->SetUseMesaClasses( offscreen );
+}
+
+vtkSmartPointer<vtkDelaunay2D> TriangulateTerrain(Func3d func)
+{
+
+	vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
+  vtkSmartPointer<vtkDoubleArray> scalars = vtkSmartPointer<vtkDoubleArray>::New();
+ 
+	unsigned int GridSize = 50;
+	for(double x = -3; x < 3; x += 6.0/GridSize)
+	{
+		for(double y = -3; y < 3; y += 6.0/GridSize)
+		{
+      double z = func(x,y);
+			vtkIdType idx = points->InsertNextPoint(x, y, z);
+      scalars -> InsertTuple(idx, &z);
+		}
+	}
+ 
+	//add the grid points to a polydata object
+	vtkSmartPointer<vtkPolyData> polydata = vtkPolyData::New();
+	polydata->SetPoints(points);
+  polydata->GetPointData() -> SetScalars(scalars);
+ 
+	//triangulate the grid points
+	vtkSmartPointer<vtkDelaunay2D> delaunay = vtkDelaunay2D::New();
+	delaunay->SetInput(polydata);
+	delaunay->Update();
+
+  return delaunay;
+}
+
+void PlotPoints(vtkSmartPointer<vtkDelaunay2D> delaunay, 
+                bool offscreen, 
+                const string& filename)
+{
+  InitOffscreen(offscreen);
+
+  //Normals for Gourad shading
+  vtkSmartPointer<vtkPolyDataNormals> normals = vtkPolyDataNormals::New();
+  normals -> SetInputConnection(delaunay ->GetOutputPort() );
+  normals -> SetFeatureAngle(60);
+
+  //Set the colours for the rendering
+  vtkSmartPointer<vtkLookupTable> lut = vtkLookupTable::New();
+  lut -> SetHueRange(0.66667, 0.0);
+  lut -> SetNumberOfColors(256);
+  lut -> SetRampToLinear();
+  lut -> Build();
+
+  
+	// map the contours to graphical primitives
+	vtkSmartPointer<vtkPolyDataMapper> contMapper = vtkPolyDataMapper::New();
+  contMapper->SetInput(normals -> GetOutput() );
+  contMapper->SetScalarRange(0,1);
+  contMapper->SetLookupTable(lut);
+  contMapper->ImmediateModeRenderingOn();
+
+ 
+  // create an actor for the contours
+	vtkSmartPointer<vtkActor> contActor = vtkActor::New();
+	contActor->SetMapper(contMapper);
+ 
+  // a renderer and render window
+	vtkSmartPointer<vtkRenderer> ren1 = vtkRenderer::New();
+	vtkSmartPointer<vtkRenderWindow> renWin = vtkRenderWindow::New();
+  if(offscreen)
+  {
+    renWin->SetOffScreenRendering(1);
+  }
+	renWin->AddRenderer(ren1);
+
+  // an interactor
+  vtkSmartPointer<vtkRenderWindowInteractor> iren 
+    = vtkRenderWindowInteractor::New();
+  iren->SetRenderWindow(renWin);  
+
+	// add the actors to the scene
+	ren1->SetBackground(0.0,0.0,0.0);
+  ren1->GetActiveCamera()->SetViewUp(0,0,1);
+  ren1->GetActiveCamera()->SetPosition(1, 1, 1);
+  ren1->GetActiveCamera()->SetParallelProjection(1);
+	ren1->AddActor(contActor);
+  ren1->ResetCamera();
+ 
+  // render an image (lights and cameras are created automatically)
+	renWin->Render();
+
+  if(offscreen)
+  {
+    // Write to file
+    vtkSmartPointer<vtkWindowToImageFilter> windowToImageFilter = 
+      vtkSmartPointer<vtkWindowToImageFilter>::New();
+    windowToImageFilter->SetInput(renWin);
+    windowToImageFilter->Update();
+ 
+    vtkSmartPointer<vtkPNGWriter> writer = vtkSmartPointer<vtkPNGWriter>::New();
+    writer->SetFileName(filename.c_str());
+    writer->SetInput(windowToImageFilter->GetOutput());
+    writer->Write();
+  }
+  else
+  {
+    //The style
+    vtkSmartPointer<vtkInteractorStyleTerrain> terrain_style 
+      = vtkInteractorStyleTerrain::New();
+ 
+    // Set the style
+    iren -> SetInteractorStyle(terrain_style);
+  
+    // begin mouse interaction
+    iren->Start(); 
+  }
+}