 	//C is column-wise ordering (column-major)
 	
 	//retrieve the depths at which to do reconstructions
-	z = mxGetPr(prhs[1]);
-	nz = mxGetM(prhs[1]) * mxGetN(prhs[1]);
+	reconstruction_z_depths = mxGetPr(prhs[1]);
+	num_z_depths = mxGetM(prhs[1]) * mxGetN(prhs[1]);
 	
-	//mexPrintf("M: %i, N: %i, nz: %i\n", M, N, nz);
+	//mexPrintf("M: %i, N: %i, num_z_depths: %i\n", M, N, num_z_depths);
 
 	//get optional parameters; use defaults otherwise
 	if (nrhs>=3)

	if (record_intensity_at_best_focus)
		cudaMalloc( (void**)&intensity_at_best_focus_gpu, sizeof(float) * M * N);

	//fillFloatKernel<<<mygrid, myblock>>>(image_magnitude, 1.0f, M, N); //don't need this

	mexPrintf("Computing reconstructions\n");

	//reconstruct the hologram at each depth in z, compare it to the previous
	//minimum intensity
	for (int i=0; i<nz; i++){
		//copy the fft spectrum to a variable that's temporary for each plane
		cudaMemcpy(reconstruction_gpu, chologram_gpu, sizeof(float2)*M*N, cudaMemcpyDeviceToDevice);

		//set the propagation parameters
		setPropagationParams( (float)(z[i]), (float)wavelength, (float)pixel_size, M, N, 
			&adudu, &cutoff2, &offsetPhaseExp, &fftscale, false);

		//apply the propagation kernel (full... could speed this up by ignoring cutoff)

}


/******************************************************************************
FFT coordinates

CUDA's FFT follows the same pattern as the popular FFTW library and arranges its
output data in a computationally-efficient but confusing pattern. The first half
of the data corresponds to the positive frequency components, while the last
half of the data corresponds to the negative frequency components. CUDA also
makes use of the fact that the component at the positive Nyquist frequency is
the same as at the negative Nyquist frequency.

The output is arranged as:

  [0, 1, 2, 3, ..., N/2 - 1, N/2, -N/2+1, -N/2+2, ... -1]

  |<------------------------>fn<----------------------->|
       positive frequency          negative frequency
           components                 components

where the N/2 sample is the Nyquist frequency component and is the same for both
the positive and negative Nyquist frequency.

Functions which compute values directly in the Fourier domain first find the
frequency component for each array entry.
******************************************************************************/


// Applies a specialized low-pass filter to data which has been 2D Fourier
// transformed.
//
// A "power" filter is a low-pass filter with a controllable ramp down, and
// follows
//
//   P(u, v) = exp(-(u^p + v^p) / sigma)
//
// where p is a positive even integer, sigma is a scale parameter, and u and
// v are normalized frequences. Normalization is against the Nyquist frequency,
// so that u and v range from [-1, 1]. Setting p = 2 gives a Gaussian low-pass
// filter; the higher p, the higher the cut-off frequency becomes and the
// sharper the ramp down.
//
// The sigma parameter can be used to set the scale, and functions like an
// additive offset for high frequencies (so that a small, controllable amount of
// the higher spatial frequencies can be kept). At |u| = |v| = 1, the filter has
// a value of P = exp(-1 / sigma).
//
// The power filter is applied directly to the 2D FFT data in fsrc, multiplying
// the magnitude of the component. The power_order is p, and should be an even
// integer. M is the number of rows and N is the number of columns in fsrc.
// The data is modified in-place.
static __global__ void PowerKernel(float2* fsrc,
                                   const int N, const int M,
								   const float power_order,
								   const float sigma) {
	int xidx = threadIdx.x + blockDim.x * blockIdx.x;
	int yidx = threadIdx.y + blockDim.y * blockIdx.y;
	if (xidx < N && yidx < M){
		int idx = yidx * N + xidx;

		//get the FFTW-shifted u,v coordinates (in numbers of samples)
		//option 4) use normalized freqs (just the magnitude) - note that uidx->u here in the nomenclature.
        // Find the sample corresponding to the Nyquist frequency.
		int nyquist_sample_x = N / 2;
		int nyquist_sample_y = M / 2;

        // Compute the normalized frequency for each pixel.
        // Each pixel's corresponding Fourier frequency is found using the
        // shifted coordinates discussed above, then normalized so that the
        // Nyquist frequency is 1.
		float normalized_freq_u = (xidx >= nyquist_sample_x ?
            (float) (N - xidx) / (float) nyquist_sample_x :
		    (float) xidx / (float) nyquist_sample_x);
		float normalized_freq_v = (yidx >= nyquist_sample_y ?
            (float) (M - yidx) / (float) nyquist_sample_y :
		    (float)yidx / (float) nyquist_sample_y);

        // Pull out the coefficient for this pixel's frequency.
		float2 value = fsrc[idx];

        // Compute the power kernel's magnitude.
		float fkeep = __expf(-( __powf(normalized_freq_u, power_order) +
		                        __powf(nonrmalized_freq_v, power_order)) /
		                     sigma);

        // Scale the FFT coefficient's magnitude and store it back to its
        // original array.
		value.x = value.x * fkeep;
		value.y = value.y * fkeep; 
		fsrc[idx] = value;
	}
}
