 	dev_scaleHalfRect<<<iDivUp(nrX*nrY,128), 128>>>(d_color_tmp, 50.0, nrX*nrY);
 	cutilSafeCall(cudaMemcpy(blue_yellow,d_color_tmp,sizeof(float)*nrX*nrY,GPUpointers?cudaMemcpyDeviceToDevice:cudaMemcpyDeviceToHost));
 
-	// shift d_stimBuf in time by 1 frame
+
+
+
+	/* ***** MOTION ENERGY MODEL ***** */
+
+	// ME responses do not depend on color responses, but on grayscale values found in d_stimBuf[]
+
+	// shift d_stimBuf in time by 1 frame, from frame i to frame i-1
 	for(int i=1;i<nrT;i++)
 		cutilSafeCall(cudaMemcpy(&d_stimBuf[nrX*nrY*(i-1)],&d_stimBuf[nrX*nrY*i],sizeof(float)*nrX*nrY,cudaMemcpyDeviceToDevice));
 
+	// allocate d_resp, which will contain the response to all 28 (nrDirs) space-time orientation at 3 (nrScales) scales
+	// for every pixel location (x,y)
 	cutilSafeCall(cudaMemset (d_resp, 0, sizeof(float)*nrX*nrY*nrDirs*nrScales));
 
-	//copy d_stimBuf to d_scalingStimBuf
+	// working copy of grayscale values: copy d_stimBuf to d_scalingStimBuf
 	cutilSafeCall(cudaMemcpy(d_scalingStimBuf,d_stimBuf,sizeof(float)*nrX*nrY*nrT,cudaMemcpyDeviceToDevice));
+
+	// compute the V1 simple cell responses at 3 different spatial scales
 	for (int scale=1; scale<=nrScales; scale++) {
+		// blur/scale the image... each time this is called stim is blurred more
+		// scale 1 == original image resolution (space/time)
 		if (scale > 1) {
 			float* tmp;
 			cutilSafeCall(cudaMalloc((void**)&tmp, sizeof(float)*nrX*nrY*nrT));

	}
}

__global__ void dev_scaleHalfRect(float *data, float scale, int len) {
	const int     tid = IMUL(blockDim.x, blockIdx.x) + threadIdx.x;
	const int threadN = IMUL(blockDim.x, gridDim.x);

	}
}

// parallel mulitplying with a scale factor
__global__ void dev_scale(float *data, float scale, int len) {
	const int     tid = IMUL(blockDim.x, blockIdx.x) + threadIdx.x;
	const int threadN = IMUL(blockDim.x, gridDim.x);
