 }
 
 
-
-
 // Returns a * sc, where a is a complex-valued float and sc is a scalar.
 // This device kernel can only runs on a CUDA device.
-__device__ cufftComplex complexScale32(cufftComplex a, float sc)
-{
+__device__ cufftComplex complexScale32(cufftComplex a, float sc) {
 	cufftComplex tempResult;
 	tempResult.x = a.x*sc;
 	tempResult.y = a.y*sc;

//
// This kernel is used for finding which z-plane maximimzes a focus metric.
//
static __global__ void compareRealKernel(float* max_metric,
                                         short* index_of_max,
                                         float* slice_metric,
                                         const int comparison_index,
                                         const int M,
                                         const int N) {
	int xidx = threadIdx.x + blockDim.x * blockIdx.x;

		int idx = yidx * ypix + xidx; //the linear index of the data
		
		//set the value
		gradient_filter_in_y[idx].x = -2.0f * v * __expf(-(u * u + v * v) /
			            (2.0f * sigma * sigma));
        // The imaginary component is zero.
        gradient_filter_in_y[idx].y = 0;
	}
}



void reconstruct(float* focus_metric, short* index_of_max_focus_metric, float* min_intensity_cpu, short* index_of_min_intensity, float* hologram_cpu, 
					 double *z, const int nz, const int M, const int N, 
					 const double wavelength, const double pixel_size,
					 const float power_filter_order, const float power_filter_offset,
					 const float steerable_filter_sigma, const float maximum_intensity,
					 float* intensity_at_best_focus_cpu, bool rintmode, bool record_intensity_at_best_focus) {
	float* hologram_gpu, *min_intensity_gpu, *simmaxdev, *gradient_magnitude_gpu, *image_magnitude, *sim, *intensity_at_best_focus_gpu;
	short* index_of_min_intensity_gpu, *index_of_max_focus_metric_gpu;
	float2* chologram_gpu, *reconstruction_gpu, *gradient_filter_in_x, *gradient_filter_in_y, *sx, *sy;
	float adudu, cutoff2, fftscale;
	float2 offsetPhaseExp;
	float sigma;

    // A block is the number of threads which are launched at each round. A size
    // of 16x16 is used to step through the image 256 threads at a time.
	dim3 myblock(16, 16);
    // TODO(nloomis): doc
	dim3 mygrid( roundUpDiv(N,16), roundUpDiv(M,16) );

	mexPrintf("Uploading hologram\n");

	//upload the hologram and convert it to a complex-valued matrix
	cudaMalloc( (void**)&hologram_gpu, sizeof(float)*M*N);
	cudaMalloc( (void**)&chologram_gpu, sizeof(float2)*M*N);
	cudaMemcpy( hologram_gpu, hologram_cpu, sizeof(float)*M*N, cudaMemcpyHostToDevice);
	float_to_complexKernel<<<mygrid, myblock>>>(chologram_gpu, hologram_gpu, M, N);
	cudaFree(hologram_gpu);

	mexPrintf("Creating memory on GPU\n");

	//create variables to track the max focus_metric, min int
	cudaMalloc( (void**)&simmaxdev, sizeof(float)*M*N);
	cudaMalloc( (void**)&index_of_max_focus_metric_gpu, sizeof(short)*M*N);
	cudaMemset(index_of_max_focus_metric_gpu, 0, sizeof(short)*M*N);
	//fillFloatKernel<<<mygrid, myblock>>>(min_intensity_gpu, FLT_MIN, M, N); //don't need this
	if (rintmode){
		cudaMalloc( (void**)&min_intensity_gpu, sizeof(float)*M*N);
		cudaMalloc( (void**)&index_of_min_intensity_gpu, sizeof(short)*M*N);
		cudaMemset(index_of_min_intensity_gpu, 0, sizeof(short)*M*N);
	}

	mexPrintf("Creating FFT plan\n");
