/////////////////////////////////////////////////////////////////////
// (c) Gatan Inc.
/////////////////////////////////////////////////////////////////////
//	This script simulates a focal series of theoretical diffractograms.
//////////////////////////////////////////////////////////////////////
//	last modified 08-July-2014 BS (comments only)

// Define required parameters; spherical abberation co-efficient, 
// incident electron wavelength and the Nyquist limit 
// and input from the user reuiqred parameters such as image size, defocus, etc.
// Note if (!GetNumber) exit(0) denotes that if cancel is selected, procedure will stop.
number pi = 3.1415926, Cs = 0.001, lambda = 0.00164, nyquist = 10   
number defocus, focalStep, width, height 							
if (!GetNumber("Start defocus (nm): ", -200, defocus)) exit (0)		
if (!GetNumber("Focal Step (nm): ", 40, focalStep)) exit (0)	
if (!GetNumber("Image size in pixels (X): ", 256, width)) exit (0)	
if (!GetNumber("Image size in pixels (Y): ", 256, height)) exit (0)	

// Post a message to the user alerting them to hit stop to finish the example.
if (!ContinueCancelDialog("Press spacebar to stop ...")) exit (0)	

					
// Calculate the square of the spatial frequency (zero frequency at image centre)
// and scale to go out to 4 times the Nyquist frequency				 
image freqSquared := RealImage("", 4, width, height)				
freqSquared = (icol-(width/2))*(icol-(width/2)) +(irow-(height/2))*(irow-(height/2))				
freqSquared = 4*nyquist*freqSquared/(width*height)					
	
// Calculate the wave aberration function resulting from the objective lens
image waveAbb := RealImage("", 4, width, height)
waveAbb = Cs*lambda*lambda*lambda*freqSquared*freqSquared*exp10(9)  

// Create an image to hold the phase contrast transfer function
// and display it
image PCTF := RealImage("Demo", 4, width, height)
ShowImage(PCTF)		

// While the space bar isn't pressed, keep recalculating and updating the PCTF with increasing focus step	
// Alternatively, one could also use: " while( !SpaceDown() )										
while( GetKey() != 32 )												
{
	// Calculate sin(phase shift) as a function of frequency, including the defocus term 
	PCTF = abs(sin(2*pi*(waveAbb+defocus*lambda*freqSquared)))		
	
	// Update the PCTF image, and increase the defocus by the focal step increment for the next loop
	UpdateImage(PCTF)												
	defocus = defocus + focalStep									
}