 			if (writeIndex+1<10) ss << "0"; 
 			ss << writeIndex+1 << "psi.fits";
 			writeFitsImg(ss.str().c_str(), dimensions, 
-				pixelsPerImage, psiImages[writeIndex]);
+				pixelsPerImage, phiImages[writeIndex]);
 			ss.str("");
 			ss.clear();
 		}
 	}
 	std::cout << "Done.\n";
 
-	/* Write psiImagess file */
+	/* Write results to file */
 	// std::cout needs to be rerouted to output to console after this...
 	std::freopen("results.txt", "w", stdout);
 	std::cout << "# t0_x t0_y theta_par theta_perp v_x v_y likelihood est_flux\n";

	trajectory best = { .xVel = 0.0, .yVel = 0.0, .lh = 0.0, 
		.x = x, .y = y, .itCount = trajectoryCount };
	
	for (int t=0; t<trajectoryCount; ++t)
	{
		float xVel = tests[t].xVel;
		float yVel = tests[t].yVel;
		float currentLikelyhood = 0.0;
		for (int i=0; i<imageCount; ++i)
		{
			currentLikelyhood += images[ i*width*height + 
				(y+int( yVel*float(i)))*width +
				 x+int( xVel*float(i)) ] / mean; 	
		}
		
		if ( currentLikelyhood > best.lh )
		{
			best.lh = currentLikelyhood;
			best.xVel = xVel;
			best.yVel = yVel;
		}		
	}	
	
	psiImagess[ y*width + x ] = best;	
}



	CUDA_CHECK_RETURN(cudaMemcpy(deviceImgTimes, imageTimes,
			sizeof(float)*imageCount, cudaMemcpyHostToDevice));

	// Copy over interleaved buffer of psi and phi images
	CUDA_CHECK_RETURN(cudaMemcpy(devicePsiPhi, interleavedPsiPhi,
		2*sizeof(float)*pixelsPerImage*imageCount, cudaMemcpyHostToDevice));
	}

	dim3 blocks(dimensions[0]/32+1,dimensions[1]/32+1);
	dim3 threads(32,32);

	// Launch Search
	searchImages<<<blocks, threads>>> (trajCount, dimensions[0], 
		dimensions[1], imageCount, padding, devicePsiPhi,
		deviceTests, deviceSearchResults, deviceImageTimes);

	// Read back results
	CUDA_CHECK_RETURN(cudaMemcpy(bestTrajects, deviceSearchResults,
				sizeof(trajectory)*pixelsPerImage, cudaMemcpyDeviceToHost));

	CUDA_CHECK_RETURN(cudaFree(deviceTests));
	CUDA_CHECK_RETURN(cudaFree(deviceImgTimes));
	CUDA_CHECK_RETURN(cudaFree(deviceSearchResults));
	CUDA_CHECK_RETURN(cudaFree(devicePsiPhi));

	
	// Sort results by likelihood
	qsort(bestTrajects, pixelsPerImage, sizeof(trajectory), compareTrajectory);
	if (debug)
	{
