diff --git a/xraysim/xraysim.cxx b/xraysim/xraysim.cxx
new file mode 100644
index 0000000..2e9a265
--- /dev/null
+++ b/xraysim/xraysim.cxx
@@ -0,0 +1,472 @@
+#include <cstdlib>
+#include <cstdio>
+#include <cstring>
+#include <unistd.h>
+
+#include <string>
+#include <iostream>
+#include <exception>
+
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+
+#include <math.h>
+
+#include "SimpleGVXR.h"
+
+#include <sys/time.h>
+#include <thread>
+
+#include <tiffio.h>
+
+float version = 3.14;
+
+// for synthetic noise 
+#include <random>
+std::mt19937 generator;
+double mean = 0.0;
+double stddev  = 1.0;
+std::normal_distribution<double> normal(mean, stddev);
+
+void usage () {
+  printf("xraysim v%f\n",version);
+  printf("usage: xraysim [-s(ourcePos) x y z] [-d(etPos) x y z] [-D(detUpVec) x y z]\n");
+  printf("               [-e(nergy) MeV] [-S(caleIntensity) min max range] [-n addgaussnoiseperc]\n");
+  printf("               [-r(otate-z) degree] [-R(otate) x y z]\n");
+  printf("               [-ct(-rotate-z-step) delta_degree] [-cr(-rotate-z-range) degree]\n");
+  printf("               [-w pixels]  [-h pixels] [-p pixelsize] [-g openglwinsize]\n");
+  printf("               [-show(_scene)] [-soft(ware_cpu_renderer)] [-o out[%%04d].tif/pgm]\n");
+  printf("               [-u density g/cm3] file.stl [-u density file2.stl]\n");
+}
+
+long int millitime() {
+  struct timeval tp;
+  gettimeofday(&tp, NULL);
+  long int ms = tp.tv_sec * 1000 + tp.tv_usec / 1000;
+  return ms;
+}
+
+int savePGM(std::string output,std::vector<std::vector<float>> image,
+            float scaleMin, float scaleMax, float scaleRange){
+  int iwidth=image.at(0).size(),
+      iheight=image.size();
+  // Save the image into a image file
+  printf("Saving Xray Image to file %s (Scale min=%f max=%f range=%f)\n",
+         output.c_str(),scaleMin,scaleMax,scaleRange);
+  // normalize to X bits / scaleRange
+  FILE *fd = fopen(output.c_str(),"w+");
+  fprintf(fd,"P2\n");
+  fprintf(fd,"%d %d\n",image.at(0).size(),image.size());
+  fprintf(fd,"%d\n",(int)scaleRange);
+  int k=0,length=image.size()*image.at(0).size();
+  for(int i=0;i<image.size();i++) {
+    for(int j=0;j<image.at(0).size();j++) {
+      float v=(image.at(i).at(j)-scaleMin)/(scaleMax-scaleMin);
+      fprintf(fd,"%d",int(v*scaleRange));
+      k++;
+      if ((k%16)==0) fprintf(fd,"\n");
+      else if (k<length) fprintf(fd," ");
+    }
+  }
+  fclose(fd);
+  return 0;
+}
+
+int saveTIFF(std::string output,std::vector<std::vector<float>> image,
+             float scaleMin, float scaleMax){
+  float scaleRange=65535;
+  int iwidth=image.at(0).size(),
+      iheight=image.size();
+  printf("Saving Xray Image to file %s (Scale min=%f max=%f range=%f)\n",
+         output.c_str(),scaleMin,scaleMax,scaleRange);
+  // Write 16 bit tiff
+  // Open the TIFF file
+  TIFF *output_image;
+  if((output_image = TIFFOpen(output.c_str(), "w")) == NULL){
+    std::cerr << "Unable to write tif file: " << output << std::endl;
+  }
+  TIFFSetField(output_image, TIFFTAG_IMAGEWIDTH, iwidth);
+  TIFFSetField(output_image, TIFFTAG_IMAGELENGTH, iheight);
+  TIFFSetField(output_image, TIFFTAG_SAMPLESPERPIXEL, 1);
+  TIFFSetField(output_image, TIFFTAG_BITSPERSAMPLE, 16);
+  TIFFSetField(output_image, TIFFTAG_ROWSPERSTRIP, iheight);
+  TIFFSetField(output_image, TIFFTAG_ORIENTATION, (int)ORIENTATION_TOPLEFT);
+  TIFFSetField(output_image, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
+  TIFFSetField(output_image, TIFFTAG_COMPRESSION, COMPRESSION_NONE);
+  TIFFSetField(output_image, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK);
+
+  tsize_t image_s;
+  short int *data=(short int*)malloc(iwidth*iheight*2);
+  int k=0;
+  for(int i=0;i<image.size();i++) {
+    for(int j=0;j<image.at(0).size();j++) {
+      float v=(image.at(i).at(j)-scaleMin)/(scaleMax-scaleMin);
+      data[k]=int(v*scaleRange);
+      k++;
+    }
+  }
+
+  if((image_s = TIFFWriteEncodedStrip(output_image, 0, (tdata_t) data, iwidth*iheight*2)) == -1)
+  {
+        std::cerr << "Unable to write tif file: " << output << std::endl;
+  }
+   
+  TIFFWriteDirectory(output_image);
+  TIFFClose(output_image);
+  return 0;
+}
+
+// usage: xraysim [-s x y z] [-e MeV] [-d x y z] [-w pixels] [-h pixels] [-p pixelsize] file.stl
+int main(int argc, char* argv[])
+{
+    int i;
+    
+
+    bool showScene=0;
+    bool softRenderer=0;
+    float sourcePosX=-40.0;
+    float sourcePosY=0.0;
+    float sourcePosZ=0.0;
+    float detPosX=10.0;
+    float detPosY=0.0;
+    float detPosZ=0.0;
+    float detUpX=0.0;
+    float detUpY=0.0;
+    float detUpZ=-1.0;
+    float hu[10] = {1000.0,1000.0,1000.0,1000.0,1000.0,
+                    1000.0,1000.0,1000.0,1000.0,1000.0};
+    float density[10] = {1.0,1.0,1.0,1.0,1.0,
+                         1.0,1.0,1.0,1.0,1.0};
+    
+    float energy  = 0.08;
+    int width     = 640;
+    int height    = 320;
+    int glwinsize = 500;
+    
+    float pixelSize = 0.5;
+    
+    float angleX   = 0.0;
+    float angleY   = 0.0;
+    float angleZ   = 0.0;
+    
+    float angleCT  = 0.0;
+    float angleCTRange = 360.0;
+    
+    float scaleK    = 0.0;
+    float scaleRange = 65535;
+    float scaleMin  = 0.0;
+    float scaleMax  = 0.0;
+    
+    float noise = 0.0;
+    
+    int parts   = 0;
+    int frameIndex  = 0;
+    int frames  = 1;
+    
+    std::string file[10];
+    std::string partLabel[10]={"A","B","C","D","E","F","G","H","I","J"};
+    std::string output="xray.tif";
+    
+    setenv("DISPLAY",":0.0",1);
+    
+    if (argc<2) { usage(); return EXIT_FAILURE; }
+    for (int i = 1; i < argc; ++i) {
+      if (std::string(argv[i]) == "-o") {
+        if (i + 1 < argc) {
+          i++;
+          output = argv[i];
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-u") {
+        if (i + 1 < argc) { 
+          i++;
+          density[parts] = std::stof(argv[i]);
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-e") {
+        if (i + 1 < argc) { 
+          i++;
+          energy = std::stof(argv[i]);
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-ct") {
+        if (i + 1 < argc) { 
+          i++;
+          angleCT = std::stof(argv[i]);
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-cr") {
+        if (i + 1 < argc) { 
+          i++;
+          angleCTRange = std::stof(argv[i]);
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-S") {
+        if (i + 3 < argc) { 
+          i++;
+          scaleMin = std::stof(argv[i++]);
+          scaleMax = std::stof(argv[i++]);
+          scaleRange = std::stof(argv[i]);
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-r") {
+        if (i + 1 < argc) { 
+          i++;
+          angleZ = std::stof(argv[i]);
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-p") {
+        if (i + 1 < argc) { 
+          i++;
+          pixelSize = std::stof(argv[i]);
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-w") {
+        if (i + 1 < argc) { 
+          i++;
+          width = std::stoi(argv[i]);
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-h") {
+        if (i + 1 < argc) { 
+          i++;
+          height = std::stoi(argv[i]);
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-g") {
+        if (i + 1 < argc) { 
+          i++;
+          glwinsize = std::stoi(argv[i]);
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-n") {
+        if (i + 1 < argc) { 
+          i++;
+          noise = std::stof(argv[i]);
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-s") {
+        if (i + 3 < argc) { 
+          i++;
+          sourcePosX = std::stof(argv[i++]);
+          sourcePosY = std::stof(argv[i++]);
+          sourcePosZ = std::stof(argv[i]);
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-d") {
+        if (i + 3 < argc) { 
+          i++;
+          detPosX = std::stof(argv[i++]);
+          detPosY = std::stof(argv[i++]);
+          detPosZ = std::stof(argv[i]);
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-D") {
+        if (i + 3 < argc) { 
+          i++;
+          detUpX = std::stof(argv[i++]);
+          detUpY = std::stof(argv[i++]);
+          detUpZ = std::stof(argv[i]);
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-R") {
+        if (i + 3 < argc) { 
+          i++;
+          angleX = std::stof(argv[i++]);
+          angleY = std::stof(argv[i++]);
+          angleZ = std::stof(argv[i]);
+          continue;
+        } else { usage(); return EXIT_FAILURE; }
+      } 
+      if (std::string(argv[i]) == "-show") {
+        showScene=1;
+        continue;
+      } 
+      if (std::string(argv[i]) == "-soft") {
+        softRenderer=1;
+        continue;
+      } 
+      if (argv[i][0]!='-') {
+        file[parts++]=argv[i];
+      }
+    }
+    try
+    {
+        if (softRenderer) setenv("LIBGL_ALWAYS_SOFTWARE","1",1);
+        // Print the libraries' version
+        std::cout << getVersionOfSimpleGVXR() << std::endl;
+        std::cout << getVersionOfCoreGVXR()   << std::endl;
+
+        printf("Creating OpenGL Context and Window\n");
+        // Create an OpenGL context
+        // createOpenGLContext(-1,3,2);
+        // createWindow(-1,1,"VULKAN");
+        createWindow();
+        
+        printf("Setting OpenGL Context Window size %d x %d \n",glwinsize,glwinsize);
+        setWindowSize(glwinsize, glwinsize);
+
+        // Set the position of the X-ray source
+        printf("Setting source position %f x %f x %f [cm]\n",sourcePosX,  sourcePosY, sourcePosZ);
+        setSourcePosition(sourcePosX,  sourcePosY, sourcePosZ, "cm");
+
+        // Set the shape of the X-ray source (here an infinitely small point)
+        usePointSource();
+
+        // The spectrum of the X-ray beam
+        // (here a monochromatic spectrum)
+        printf("Setting monochromatic energy %f MeV\n",energy);
+        setMonoChromatic(energy, "MeV", 1000);
+
+        // Set the position of the X-ray detector (film/cassette)
+        printf("Setting detector position %f x %f x %f [cm]\n",detPosX, detPosY, detPosZ);
+        setDetectorPosition(detPosX, detPosY, detPosZ, "cm");
+
+        // Set the orientation of the X-ray detector (film/cassette)
+        printf("Setting detector up vector %f x %f x %f\n",detUpX, detUpY, detUpZ);
+        setDetectorUpVector(detUpX, detUpY, detUpZ);
+
+        // Set the number of pixel of the X-ray detector
+        printf("Setting detector pixels %d x %d [%f mm]\n",width,height,pixelSize);
+        setDetectorNumberOfPixels(width, height);
+
+        // Set the distance between two successive pixels
+        setDetectorPixelSize(pixelSize, pixelSize, "mm");
+
+        printf("Reading polygon mesh files...\n");
+        // Load a polygon mesh from a STL file as a scenegraph
+ //       loadSceneGraph(file, "mm");
+        for(int i=0;i<parts;i++) {
+          printf("[%d] %s (%s) <= %f (g/cm3)\n",i,file[i].c_str(),partLabel[i].c_str(),density[i]);
+          loadMeshFile(partLabel[i], file[i], "mm");
+          // moveToCentre(partLabel[i]);
+          setHU(partLabel[i], 1000);
+          setDensity(partLabel[i], density[i], "g/cm3");
+        }
+        // Set the material properties of all the nodes within the scenegraph
+//        for (unsigned int i = 0; i < getNumberOfChildren("root"); ++i)
+//        {
+//            std::string label = getChildLabel("root", i);
+//            printf("Setting material property for %s <= %f (hu)\n",label.c_str(),hu);
+//            moveToCentre(label);
+//            setHU(label, hu);
+//        }
+
+
+/*
+        loadSceneGraph("second.stl", "mm");
+
+        for (unsigned int i = 0; i < getNumberOfChildren("root"); ++i)
+        {
+            std::string label = getChildLabel("root", i);
+            printf("Setting material property for %s <= %f (hu)\n",label.c_str(),hu);
+            moveToCentre(label);
+            setHU(label, hu);
+        }
+*/
+        if (angleX!=0.0 || angleY!=0 || angleZ!=0) {
+          printf("Rotating scene by %f %f %f\n",angleX,angleY,angleZ);
+          if (angleX!=0.0) rotateScene(angleX, 1.0, 0.0, 0.0);
+          if (angleY!=0.0) rotateScene(angleY, 0.0, 1.0, 0.0);
+          if (angleZ!=0.0) rotateScene(angleZ, 0.0, 0.0, 1.0);
+        }
+
+        if (angleCT!=0.0) frames=(int)((float)angleCTRange/angleCT);
+        int autoScale=0;
+        if (scaleMin==0.0 && scaleMax==0.0) autoScale=1;
+        
+        for(frameIndex=0;frameIndex<frames;frameIndex++) {
+          char path[1000];
+          snprintf(path, sizeof(path), output.c_str(), frameIndex);
+          std::string pathStr = path;
+          
+          printf("[%d] Computing Xray Image ...\n",frameIndex);
+          // Compute an X-ray image(planar projection/radiograph)
+          
+          long int start=millitime();
+          computeXRayImage();
+          long int stop=millitime();
+
+          printf("[%d] Computing Xray Image took %d ms...\n",frameIndex,stop-start);
+          
+          std::vector<std::vector<float>> image=getLastXRayImage();
+          
+          if (noise!=0.0) {
+            // add relative gaussian noise (rough approximation of xray noise)
+            for(int i=0;i<image.size();i++) {
+              for(int j=0;j<image.at(0).size();j++) {
+                float v = image.at(i).at(j),
+                      n = normal(generator)*noise,
+                      vn = v+abs(v)*n;
+                image.at(i).at(j)=vn;
+              }                            
+            }
+          }
+          int iwidth  = image.at(0).size(),
+              iheight = image.size();
+          if (autoScale) {
+            // auto-scale
+            scaleMin=scaleMax=image.at(0).at(0);
+            for(int i=0;i<image.size();i++) {
+              for(int j=0;j<image.at(0).size();j++) {
+                float v=image.at(i).at(j);
+                scaleMin=std::min(scaleMin,v);
+                scaleMax=std::max(scaleMax,v);
+              }
+            }
+          }
+              
+          float imin=1E30,imax=-1E30;
+          for(int i=0;i<image.size();i++) {
+            for(int j=0;j<image.at(0).size();j++) {
+              imin=std::min(imin,image.at(i).at(j));
+              imax=std::max(imax,image.at(i).at(j));
+            }
+          }
+          printf("[%d] Xray Image: %d x %d pixels [%f,%f]\n",frameIndex,image.at(0).size(),image.size(),imin,imax);
+          
+          if (path[strlen(path)-1]=='f') {
+            saveTIFF(pathStr,image,scaleMin,scaleMax);
+          }
+          if (path[strlen(path)-1]=='m') {
+            savePGM(pathStr,image,scaleMin,scaleMax,scaleRange);
+          }
+          if (angleCT!=0.0) {
+            printf("[%d] Rotating scene by 0 0 %f\n",frameIndex,angleCT);
+            rotateScene(angleCT, 0.0, 0.0, 1.0);
+          }
+        }
+        // saveLastSinogram();
+        // saveLastProjectionSet();
+        // saveLastLBuffer();
+    }
+    catch (const std::exception& err)
+    {
+        std::cerr << "FATAL ERROR:\t" << err.what() << std::endl;
+        return EXIT_FAILURE;
+    }
+
+    if (showScene) {
+      displayScene();
+      std::cout << "press any key to exit...";
+      std::cin.get();
+    }
+        
+    return EXIT_SUCCESS;
+}