function output = generate_simulated_spots(npixels, digitizeFLAG, ...
                          IntensityNoiseFLAG, positionNoiseFLAG)
%
% a function to generate an image of spots to centroid and to provide the original 
% locations of the spots. 
% 
% input
% -----
% npixels - pixels in output image
% digitizeFLAG:          0 - floating point array is output
%                        1 - intensity array is rounded down to nearest
%                       integer (12-bit)
% IntensityNoiseFLAG:     0 - no photo-electron shot noise
%                         1 - photo-electron shot noise is added
% positionNoiseFLAG:      0 - position noise off
%                  :      1 - position noise ON - deviation = positionDeviationP
% output
% ------
% output.image      - the 2D intensity array [npixels, npixels]
% output.positions  - the nominal positions of the spots: [N, 2] array 
%
% EXAMPLE
% --------
% 
% output = generate_simulated_spots(1024, 1, 1, 0, 0)
%   - generates an image with photo-electron intensity noise, but with
%   the theoretical spot positions undeviated from a regular array. Fast
%   Gaussian intensity calculated.
%
% Gaussian Integration 
% ---------------------
% 
% - the intensity, Ip, in a pixel is given by the following
% Ip = integrate( exp{ -2 [(x-xc)^2 + (y-yc)^2]/w^2 }, {x, xa, xb}, {y, ya, yb} )
%    = 1/8 * pi * w^2
%       * [ erf( sqrt(2)*(x0-xc)/w ) - erf( sqrt(2)*(x1-xc)/w ) ]
%       * [ erf( sqrt(2)*(y0-yc)/w ) - erf( sqrt(2)*(y1-yc)/w ) ]
%
%
%  Coordinate representation of pixels
%  -----------------------------------
%
%   |          |          |
%   |          |          |
% 2 +----------+----------+---
%   |          |          |
%   |          |          |
%   |          |      <==========  this pixel has indicies [2, 2] but is      
%   |          |          |        said to have coordinates [1.5, 1.5], 
%   |          |          |        corresponding to the center.
% 1 +----------+----------+---
%   |          |          |
%   |          |          |
%   |          |          |
%   |          |          |
%   |          |          |
% 0 +----------+----------+---
%   0          1          2
%  
%
% Written by Aidan Brooks. 19th Mar 2010

% define the physical parameters of the HWS
pixelSize = 12.0E-6;        % pixel size in m
maxElectronsPerPixel = 1.6E5;    % around 160,000 electrons before saturation
maxElectronsInExposure = 1.5E5;
holePitchM = 430E-6;         % separation of the holes in the Hartmann plate
holePitchP = holePitchM/pixelSize;         % separation of the holes in the Hartmann plate in Pixels
spotRadiusM = 75.0E-6;       % approximate spot radius on CCD = 1/2 hole size
spotRadiusP = spotRadiusM/pixelSize;
CCDbackground = 50.0;       % background level in 12-bit representation
positionDeviationP = 0.1;   % the standard deviation of random spot position deviation (in pixels)
intensityUnitsPerElectron = 4096/maxElectronsPerPixel;

% buffer (in pixels) around the edges 
buffer = ceil(holePitchP);

% generate the array of positions and the image if the array is big enough
if (npixels - 2*buffer - holePitchP > 0)
   % generate the blank image 
   image_arr = zeros(npixels, npixels);
   
   % generate the coordinate arrays - the pixel value is the center of the
   % pixel, therefore the first pixel is 0.5, the next 1.5, etc.
   x = 1:npixels;
   x2D = (x'*ones(1,npixels))' - 0.5;
   y2D = x2D';
    
   % get the number of spots across the top of a recti-linear array
   nspot1D = floor((npixels - 2*buffer)/holePitchP) + 1;
   % create an empty array of the nominal spot positions
   nominal_positions = [];
   
   % loop through the spots
   for ii = 1:nspot1D          % ii-th column
       for jj = 1:nspot1D      % jj-th row
           spot_position = ([(ii - 1), (jj-1)]*holePitchP) + [buffer buffer];
           % if its flag is present then add position noise
           if (positionNoiseFLAG ~= 0)
               spot_position = spot_position + ...
                               random('Normal', 0, positionDeviationP, 1, 2);
           end
           
           % add the spot_position to the nominal position list
           nominal_positions = [nominal_positions; spot_position];
           
           %-----------------------------------------------
           % IMAGE
           % in a square around the hole, add the intensity
           xc = spot_position(1);
           yc = spot_position(2);
           imin = floor(xc - holePitchP/2);
           imax =  ceil(xc + holePitchP/2);
           jmin = floor(yc - holePitchP/2);
           jmax =  ceil(yc + holePitchP/2); 
           

           % get the limits of integration
           % NOTE: [row number (= Y), column number (= X)]
           x0 = x2D(jmin:jmax, imin:imax) - 0.5;
           x1 = x2D(jmin:jmax, imin:imax) + 0.5;
           y0 = y2D(jmin:jmax, imin:imax) - 0.5;
           y1 = y2D(jmin:jmax, imin:imax) + 0.5;

           % get the integrated gaussian using erf
           intGauss = (1/8)*pi*spotRadiusP^2 * ...
               ( erf( sqrt(2) * (x0 - xc)/spotRadiusP ) - ...
                 erf( sqrt(2) * (x1 - xc)/spotRadiusP )).* ...
               ( erf( sqrt(2) * (y0 - yc)/spotRadiusP )  - ...
                 erf( sqrt(2) * (y1 - yc)/spotRadiusP ));

           %  add the integrated sub-array to the total image
           image_arr(jmin:jmax, imin:imax) = image_arr(jmin:jmax, imin:imax) + intGauss;
           %-----------------------------------------------          

       
       
       end
   end
   
   % rescale the image to be in photoelectrons
   image_arr = maxElectronsInExposure*image_arr/max(max(image_arr)); % in photo-electrons
   
   % add the intensity noise - a normal distribution where the deviation in photo-electrons is
   % proportional to the sqrt of the number of the photo-electrons
   if (IntensityNoiseFLAG ~= 0)
      image_arr = image_arr + ...
          random('Normal', 0, sqrt(image_arr), npixels, npixels);       
   end
   
   % rescale the image to be in 12-bit digital units
   image_arr = image_arr*intensityUnitsPerElectron;
   
   % add the background level
   image_arr = image_arr + CCDbackground;
   
   % digitize the image if this is requested.
   if (digitizeFLAG ~= 0)
       image_arr = floor(image_arr);
   end
   
   
   output.positions = nominal_positions;
   output.image = image_arr;
else
    disp(['The number of pixels must be greater than ', num2str(ceil(2*buffer + holePitchP))])
    output = 0;
end