diff --git a/calc_grad_perc_dev.cc b/calc_grad_perc_dev.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f0d00ef711debdd1efc62f1b3cd97fca29eab833
--- /dev/null
+++ b/calc_grad_perc_dev.cc
@@ -0,0 +1,110 @@
+/*  calc_grad_perc_dev.cc
+
+    Mark Jenkinson, FMRIB Image Analysis Group
+
+    Copyright (C) 2012 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+// Calculates the percentage deviation in the gradient fields due to non-linearities
+// Relies on the output of gradient_unwarp.py
+
+#define _GNU_SOURCE 1
+#define POSIX_SOURCE 1
+
+#include "newimage/newimageall.h"
+#include "miscmaths/miscmaths.h"
+#include "utils/options.h"
+
+using namespace MISCMATHS;
+using namespace NEWIMAGE;
+using namespace Utilities;
+
+// The two strings below specify the title and example usage that is
+//  printed out as the help or usage message
+
+string title="calc_grad_perc_dev (Version 1.0)\nCopyright(c) 2003, University of Oxford (Mark Jenkinson)";
+string examples="calc_grad_perc_dev [options] --fullwarp=<warp image> --out=<basename>";
+
+// Each (global) object below specificies as option and can be accessed
+//  anywhere in this file (since they are global).  The order of the
+//  arguments needed is: name(s) of option, default value, help message,
+//       whether it is compulsory, whether it requires arguments
+// Note that they must also be included in the main() function or they
+//  will not be active.
+
+Option<bool> verbose(string("-v,--verbose"), false, 
+		     string("switch on diagnostic messages"), 
+		     false, no_argument);
+Option<bool> help(string("-h,--help"), false,
+		  string("display this message"),
+		  false, no_argument);
+Option<string> fullwarp(string("--fullwarp"), string(""),
+		  string("full warp image from gradient_unwarp.py"),
+		  true, requires_argument);
+Option<string> outname(string("-o,--out"), string(""),
+		  string("output basename"),
+		  true, requires_argument);
+int nonoptarg;
+
+////////////////////////////////////////////////////////////////////////////
+
+// Local functions
+
+// for example ... print difference of COGs between 2 images ...
+int do_work(int argc, char* argv[]) 
+{
+  volume4D<float> fw, xpd, ypd, zpd;
+  read_volume4D(fw,fullwarp.value());
+  gradient(fw[0],xpd);
+  xpd/=fw.xdim();
+  gradient(fw[1],ypd);
+  ypd/=fw.ydim();
+  gradient(fw[2],zpd);
+  zpd/=fw.zdim();
+  string bname=fslbasename(outname.value());
+  save_volume4D(xpd,bname+"_x");
+  save_volume4D(ypd,bname+"_y");
+  save_volume4D(zpd,bname+"_z");
+  return 0;
+}
+
+////////////////////////////////////////////////////////////////////////////
+
+int main(int argc,char *argv[])
+{
+
+  Tracer tr("main");
+  OptionParser options(title, examples);
+
+  try {
+    // must include all wanted options here (the order determines how
+    //  the help message is printed)
+    options.add(fullwarp);
+    options.add(outname);
+    options.add(verbose);
+    options.add(help);
+    
+    nonoptarg = options.parse_command_line(argc, argv);
+
+    // line below stops the program if the help was requested or 
+    //  a compulsory option was not set
+    if ( (help.value()) || (!options.check_compulsory_arguments(true)) )
+      {
+	options.usage();
+	exit(EXIT_FAILURE);
+      }
+    
+  }  catch(X_OptionError& e) {
+    options.usage();
+    cerr << endl << e.what() << endl;
+    exit(EXIT_FAILURE);
+  } catch(std::exception &e) {
+    cerr << e.what() << endl;
+  } 
+
+  // Call the local functions
+
+  return do_work(argc,argv);
+}
+