diff --git a/include/openmc/constants.h b/include/openmc/constants.h
index d249e4cb5..c2b52b25a 100644
--- a/include/openmc/constants.h
+++ b/include/openmc/constants.h
@@ -82,7 +82,10 @@ constexpr double EXTSRC_REJECT_FRACTION {0.05};
 constexpr double PI {3.1415926535898};
 const double SQRT_PI {std::sqrt(PI)};
 constexpr double INFTY {std::numeric_limits<double>::max()};
-
+constexpr double TWO {2.0};
+constexpr double HALF {0.5};
+ 
+ 
 // Physical constants
 constexpr double MASS_NEUTRON     {1.00866491588}; // mass of a neutron in amu
 constexpr double MASS_NEUTRON_EV  {939.5654133e6}; // mass of a neutron in eV/c^2
diff --git a/src/plot.cpp b/src/plot.cpp
index 925fbea05..bb4f1f756 100644
--- a/src/plot.cpp
+++ b/src/plot.cpp
@@ -2,6 +2,8 @@
 #include "openmc/constants.h"
 #include "openmc/settings.h"
 #include "openmc/error.h"
+#include "openmc/particle.h"
+#include "openmc/geometry.h"
 
 namespace openmc {
 
@@ -9,17 +11,21 @@ const int RED   = 1;
 const int GREEN = 2;
 const int BLUE  = 3;
 
-  int openmc_plot_geometry() {
-    int err;
+//===============================================================================
+// RUN_PLOT controls the logic for making one or many plots
+//===============================================================================
+  
+int openmc_plot_geometry() {
+  int err;
 
-    for(auto i : n_plots) {
-      ObjectPlot* pl = plots[i];
-
-      std::stringstream ss;
-      ss << "Processing plot " << pl->id_ << ": "
-         << pl->path_plot_ << "...";
+  for(auto i : n_plots) {
+    ObjectPlot* pl = plots[i];
+    
+    std::stringstream ss;
+    ss << "Processing plot " << pl->id_ << ": "
+       << pl->path_plot_ << "...";
       write_message(ss.str(), 5);
-
+      
       if (pl->type_ == PLOT_TYPE::SLICE) {
         // create 2D image
         // create_ppm(pl);
@@ -30,11 +36,80 @@ const int BLUE  = 3;
         continue;
       }
       
-    }
-
-    return 0;
   }
-    
+  
+  return 0;
+}
+
+//===============================================================================
+// CREATE_PPM creates an image based on user input from a plots.xml <plot>
+// specification in the portable pixmap format (PPM)
+//===============================================================================
+
+void create_ppm(ObjectPlot* pl) {
+
+  int width = pl->pixels_[0];
+  int height = pl->pixels_[1];
+
+  double in_pixel = (pl->width_[0])/double(width);
+  double out_pixel = (pl->width_[1])/double(height);
+
+  std::vector< std::vector< std::vector<int>>> data;
+
+  data.resize(width);
+  for (auto i : data) {
+    i.resize(height);
+    for (auto j : i) { j.resize(3); }
+  }
+
+  int in_i, out_i;
+  double xyz[3];
+  if (pl->basis_ == PLOT_BASIS::XY) {
+    in_i = 0;
+    out_i = 1;
+    xyz[0] = pl->origin_[0] - pl->width_[0] / TWO;
+    xyz[1] = pl->origin_[1] - pl->width_[1] / TWO;
+    xyz[2] = pl->origin_[2];
+  } else if (pl->basis_ == PLOT_BASIS::XZ) {
+    in_i = 0;
+    out_i = 2;
+    xyz[0] = pl->origin_[0] - pl->width_[0] / TWO;
+    xyz[1] = pl->origin_[1];
+    xyz[2] = pl->origin_[2] - pl->width_[1] / TWO;
+  } else if (pl->basis_ == PLOT_BASIS::YZ) {
+    in_i = 1;
+    out_i = 2;
+    xyz[0] = pl->origin_[0];
+    xyz[1] = pl->origin_[1] - pl->width_[0] / TWO;
+    xyz[2] = pl->origin_[2] - pl->width_[1] / TWO;
+  }
+
+  double dir[3];
+  Particle *p = new Particle();
+  p->initialize();
+  std::copy(xyz, xyz+3, p->coord[0].xyz);
+  std::copy(dir, dir+3, p->coord[0].uvw);
+  p->coord[0].universe = openmc_root_universe;
+
+  // local variables
+  int rgb[3];
+  int id;
+  for(int y = 0; y < height; y++) {
+    p->coord[0].xyz[out_i] = xyz[out_i] - out_pixel*(y);
+    for(int x = 0; x < width; x++) {
+      p->coord[0].xyz[in_i] = xyz[in_i] + in_pixel*(x);
+      // position_rgb(p, pl, rgb, id);
+
+      std::copy(rgb, rgb+3, &(data[x][y][0]));
+    }
+  }
+
+  //output_ppm(pl, data);
+  
+}
+
+  
+  
 void
 voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace, hid_t* dset,
            hid_t* memspace)