diff --git a/stim/cuda/cuda_texture.cuh b/stim/cuda/cuda_texture.cuh
index b967cd5..0ae6bf3 100644
--- a/stim/cuda/cuda_texture.cuh
+++ b/stim/cuda/cuda_texture.cuh
@@ -7,7 +7,7 @@
 #include <cuda_runtime.h>
 #include <cublas_v2.h>
 #include <stdio.h>
-#include <GL/glew.h>
+//#include <GL/glew.h>
 #include <GL/glut.h>
 #include <sstream>
 #include <stim/visualization/colormap.h>
diff --git a/stim/cuda/cudatools/glbind.h b/stim/cuda/cudatools/glbind.h
index 85a8f31..4ef1e2c 100644
--- a/stim/cuda/cudatools/glbind.h
+++ b/stim/cuda/cudatools/glbind.h
@@ -1,8 +1,9 @@
 #ifndef RTS_GL_BIND_H
 #define RTS_GL_BIND_H
 
-#include <GL/glew.h>
-#include <GL/gl.h>
+//#include <GL/glew.h>
+#include <GL/GL.h>
+#include <GL/GLU.h>
 
 #include <stdio.h>
 #include <cstring>
@@ -18,12 +19,12 @@ static void InitGLEW()
 {
 	//Initialize the GLEW toolkit
 
-	GLenum err = glewInit();
+	/*GLenum err = glewInit();
 	if(GLEW_OK != err)
 	{
 		printf("Error starting GLEW.");
 	}
-	fprintf(stdout, "Status: Using GLEW %s\n", glewGetString(GLEW_VERSION));
+	fprintf(stdout, "Status: Using GLEW %s\n", glewGetString(GLEW_VERSION));*/
 }
 
 static void cudaSetDevice(int major = 1, int minor = 3)
diff --git a/stim/gl/gl_spider.h b/stim/gl/gl_spider.h
index b4f7afd..6b893fa 100644
--- a/stim/gl/gl_spider.h
+++ b/stim/gl/gl_spider.h
@@ -1,7 +1,7 @@
 #ifndef STIM_GL_SPIDER_H
 #define STIM_GL_SPIDER_H
 
-#include <GL/glew.h>
+//#include <GL/glew.h>
 #include <GL/glut.h>
 #include <cuda.h>
 #include <cuda_gl_interop.h>
@@ -55,13 +55,13 @@ class gl_spider : public virtual gl_texture<T>
 	private:
 
 		#ifdef TIMING
-			double branch_time = 0;
-			double direction_time = 0;
-			double position_time = 0;
-			double size_time = 0;
-			double cost_time = 0;
-			double network_time = 0;
-			double hit_time = 0;
+			double branch_time;// = 0;
+			double direction_time;// = 0;
+			double position_time;// = 0;
+			double size_time;// = 0;
+			double cost_time;// = 0;
+			double network_time;// = 0;
+			double hit_time;// = 0;
 		#endif
 
 		//
@@ -99,7 +99,7 @@ class gl_spider : public virtual gl_texture<T>
 		int numSamplesPos;
 		int numSamplesMag;
 
-		float stepsize = 5.0;			//Step size.
+		float stepsize;// = 5.0;			//Step size.
 //		float stepsize = 3.0;			//Step size.
 		int current_cost;			//variable to store the cost of the current step.
 		
@@ -121,8 +121,8 @@ class gl_spider : public virtual gl_texture<T>
 		stim::vec3<float> ups;
 		stim::vec3<float> ds;
 
-		static const float t_length = 16.0;
-
+		//static const float t_length = 16.0;
+		float t_length;
 
 		//cuda texture variables that keep track of the binding.
 		stim::cuda::cuda_texture t_dir;
@@ -929,6 +929,18 @@ class gl_spider : public virtual gl_texture<T>
 		void
 		attachSpider(GLuint id)
 		{
+#ifdef TIMING
+			branch_time = 0;
+			direction_time = 0;
+			position_time = 0;
+			size_time = 0;
+			cost_time = 0;
+			network_time = 0;
+			hit_time = 0;
+#endif
+			stepsize = 3.0;
+			t_length = 16.0;
+
 			srand(100);	
 			texID = id;
 			//GenerateFBO(16, numSamples*8);
diff --git a/stim/math/circle.h b/stim/math/circle.h
index a395f91..92f9335 100644
--- a/stim/math/circle.h
+++ b/stim/math/circle.h
@@ -5,6 +5,7 @@
 #include <stim/math/plane.h>
 #include <stim/math/vector.h>
 #include <stim/math/triangle.h>
+#include <stim/math/constants.h>
 #include <assert.h>
 #include <algorithm>
 #include <iostream>
@@ -153,8 +154,8 @@ public:
 		float step = 360.0/(float) n;
 		for(float j = 0; j <= 360.0; j += step)
 		{
-			y = 0.5*cos(j*2.0*M_PI/360.0)+0.5;
-			x = 0.5*sin(j*2.0*M_PI/360.0)+0.5;
+			y = 0.5*cos(j*stim::TAU/360.0)+0.5;
+			x = 0.5*sin(j*stim::TAU/360.0)+0.5;
 			result.push_back(p(x,y));
 		}
 		return result;
@@ -167,8 +168,8 @@ public:
 	p(T theta)
 	{
 		T x,y;
-		y = 0.5*cos(theta*2.0*M_PI/360.0)+0.5;
-		x = 0.5*sin(theta*2.0*M_PI/360.0)+0.5;
+		y = 0.5*cos(theta*stim::TAU/360.0)+0.5;
+		x = 0.5*sin(theta*stim::TAU/360.0)+0.5;
 		return p(x,y);
 	}
 
diff --git a/stim/math/complex.h b/stim/math/complex.h
index 237b3ca..3ccc7d0 100644
--- a/stim/math/complex.h
+++ b/stim/math/complex.h
@@ -11,7 +11,7 @@
 
 namespace stim
 {
-    enum complexComponentType {complexReal, complexImaginary, complexMag};
+    enum complexComponentType {complexReal, complexImaginary, complexMag, complexIntensity};
 
 template <class T>
 struct complex
@@ -346,6 +346,13 @@ static void abs(T* m, complex<T>* c, size_t n){
 		m[i] = c[i].abs();
 }
 
+/// Calculate the intensity of an array of complex values
+template<typename T>
+static void intensity(T* m, complex<T>* c, size_t n){
+	for(size_t i = 0; i < n; i++)
+		m[i] = pow(c[i].abs(), 2);
+}
+
 }	//end RTS namespace
 
 //addition
diff --git a/stim/optics/scalarfield.h b/stim/optics/scalarfield.h
index 397f49e..31b74f2 100644
--- a/stim/optics/scalarfield.h
+++ b/stim/optics/scalarfield.h
@@ -291,6 +291,8 @@ public:
 
 	/// Propagate the field along its orthogonal direction by a distance d
 	void propagate(T d, T k){
+		//X[2] += d;																		//move the plane position along the propagation direction
+		//Y[2] += d;
 		cpu_scalar_propagate(E, E, X.len(), Y.len(), d, k, R[0], R[1]);
 	}
 
@@ -299,6 +301,35 @@ public:
 		cpu_scalar_lowpass(E, E, X.len(), Y.len(), highest, R[0], R[1]);
 	}
 
+	/// Crop an image based on a given padding parameter (crop out the center)
+	void crop(size_t padding, stim::scalarfield<T>& cropped){
+		size_t Cx = R[0] / (2 * padding + 1);										//calculate the size of the cropped image based on the padding value
+		size_t Cy = R[1] / (2 * padding + 1);
+
+		if(cropped.R[0] != Cx || cropped.R[1] != Cy){
+			std::cout<<"Error: cropped field resolution ("<<cropped.R[0]<<" x "<<cropped.R[1]<<") does not match the required resolution ("<<Cx<<" x "<<Cy<<")."<<std::endl;
+			exit(1);
+		}
+
+		if(loc == CPUmem){
+			size_t x, y;
+			size_t sx, sy, si, di;
+			for(y = 0; y < Cy; y++){
+				sy = y + Cy * padding;								//calculate the y-index into the source image
+				for(x = 0; x < Cx; x++){
+					sx = x + Cx * padding;							//calculate the x-index into the source image
+					si = sy * R[0] + sx;							//calculate the 1D index into the source image
+					di = y * Cx + x;
+					cropped.E[di] = E[si];
+				}
+			}
+		}
+		else{
+			std::cout<<"Error: cropping not supported for GPU memory yet."<<std::endl;
+			exit(1);
+		}
+	}
+
 	std::string str(){
 		std::stringstream ss;
 		ss<<rect<T>::str()<<std::endl;
@@ -367,11 +398,49 @@ public:
 			break;
 		case complexImaginary:
 			stim::imag(image, E, size());
+			break;
+		case complexIntensity:
+			stim::intensity(image, E, size());
+			break;
 		}
 		stim::cpu2image(image, filename, R[0], R[1], cmap);			//save the resulting image
 		free(image);												//free the real image
 	}
 
+	void image(T* img, stim::complexComponentType type = complexMag){
+		if(loc == GPUmem) to_cpu();									//if the field is in the GPU, move it to the CPU
+
+		switch(type){												//get the specified component from the complex value
+		case complexMag:
+			stim::abs(img, E, size());
+			break;
+		case complexReal:
+			stim::real(img, E, size());
+			break;
+		case complexImaginary:
+			stim::imag(img, E, size());
+			break;
+		case complexIntensity:
+			stim::intensity(img, E, size());
+			break;
+		}
+		//stim::cpu2image(image, filename, R[0], R[1], cmap);			//save the resulting image
+		//free(image);												//free the real image
+	}
+
+	//adds the field intensity to the output array (useful for calculating detector response to incoherent fields)
+	void intensity(T* out){
+		if(loc == GPUmem) to_cpu();									//if the field is in the GPU, move it to the CPU
+		T* image = (T*) malloc( sizeof(T) * size() );				//allocate space for the real image
+		stim::intensity(image, E, size());							//calculate the intensity
+		
+		size_t N = size();											//calculate the number of elements in the field
+		for(size_t n = 0; n < N; n++)								//for each point in the field
+			out[n] += image[n];										//add the field intensity to the output image
+
+		free(image);												//free the temporary intensity image
+	}
+
 };				//end class scalarfield
 }
 
diff --git a/stim/visualization/cylinder.h b/stim/visualization/cylinder.h
index 6dec73b..d6dc68c 100644
--- a/stim/visualization/cylinder.h
+++ b/stim/visualization/cylinder.h
@@ -203,8 +203,10 @@ class cylinder
 		T
 		r(T pvalue)
 		{
-			if(pvalue < 0.0 || pvalue > 1.0)
-				return;
+			if(pvalue < 0.0 || pvalue > 1.0){
+				std::cerr<<"Error, value "<<pvalue<<" is outside of [0 1]."<<std::endl;
+				exit(1);
+			}
 			T l = pvalue*L[L.size()-1];
 			int idx = findIdx(l);
 			return (e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*((l-L[idx])/(L[idx+1]- L[idx])));
diff --git a/stim/visualization/glObj.h b/stim/visualization/glObj.h
index 1ea4e90..9f66d6a 100644
--- a/stim/visualization/glObj.h
+++ b/stim/visualization/glObj.h
@@ -1,7 +1,7 @@
 #ifndef STIM_GLOBJ_H
 #define STIM_GLOBJ_H
 
-#include <GL/glew.h>
+//#include <GL/glew.h>
 #include <GL/glut.h>
 #include <stim/math/vector.h>
 #include <stim/visualization/obj.h>
--
libgit2 0.21.4