function [calim] = ppfakecalim(xpix,ypix,pixsize,spotdist,spotsize,imtilt,maxlev,noiselev,bglev,bgnoise)

% function [calim] = ppfakecalim(xpix,ypix,pixsize,spotdist,spotsize,imtilt,noiselev,bglev,bgnoise)
%
%    PPFAKECALIM - generate a fake pepperpot calibration image with 4
%    Gaussian light spots.
%
%    [calim] = ppfakecalim(xpix,ypix,pixsize)
%    [...] = ppfakecalim(xpix,ypix,pixsize,spotdist)
%    [...] = ppfakecalim(xpix,ypix,pixsize,spotdist,spotsize)
%    [...] = ppfakecalim(xpix,ypix,pixsize,spotdist,spotsize,imtilt)
%    [...] = ppfakecalim(xpix,ypix,pixsize,spotdist,spotsize,imtilt,maxlev)
%    [...] = ppfakecalim(xpix,ypix,pixsize,spotdist,spotsize,imtilt,maxlev,noiselev)
%    [...] = ppfakecalim(xpix,ypix,pixsize,spotdist,spotsize,imtilt,maxlev,noiselev,bglev)
%    [...] = ppfakecalim(xpix,ypix,pixsize,spotdist,spotsize,imtilt,maxlev,noiselev,bglev,bgnoise)
%
%    PPFAKECALIM generates fake pepperpot calibration images, based upon
%    the new 4-light calibration system.  This is designed to test the
%    spot-finding used in the calibration calculation code by producing a
%    single image with 4 Gaussian light spots.
%
%    The calibration image is always assumed to be centred, with each spot
%    at an equal distance from the centre, which assumes that the real
%    calibration spots are symmetrical about the central axis.  It is
%    possible to specify a spot spacing in millimetres (see below): this is
%    the distance from top-to-bottom or left-to-right, not the short
%    distance (top-to-left etc.), since the light spots are assumed to be
%    top, bottom, left and right.  The real spot-to-spot distance is 140mm.
%
%    The sub-function GAUSSPIXSPOT is used to generate the individual light
%    distribution data for each spot: this assumes that each light spot
%    produces a Gaussian light distribution.
%
%    [calim] = ppfakecalim(xpix,ypix,pixsize) - output a fake pepperpot
%    calibration image, CALIM, based on the following input parameters:
%
%        XPIX - number of horizontal pixels in the image
%        YPIX - number of vertical pixels in the image
%        PIXSIZE - pixel size, in NUMBER OF PIXELS PER MM
%
%    This assumes that the default distance between top-and-bottom or
%    left-and-right light spots is 140 mm, and that the size of each light
%    spot is 500 microns (corresponding to a radius of 1 sigma).  It also
%    assumes that there is zero noise on both the background and the spots
%    themselves and that the image is perfectly aligned to the camera.
%
%    [...] = ppfakecalim(xpix,ypix,pixsize,spotdist) - specify the
%    distance, in mm, between top-and-bottom or left-and-right light spots 
%    with SPOTDIST (default is 140 mm).
%
%    [...] = ppfakecalim(xpix,ypix,pixsize,spotdist,spotsize) - also
%    specify the 2 sigma width of the Gaussian distribution, in mm, of each
%    light spot with SPOTSIZE ie. half width is 1 sigma (default width is
%    500 microns).
%
%    [...] = ppfakecalim(xpix,ypix,pixsize,spotdist,spotsize,imtilt) - tilt
%    the image with respect to the camera by IMTILT radians: this is given
%    from the positive X-axis.
%
%    [...] = ppfakecalim(xpix,ypix,pixsize,spotdist,spotsize,imtilt,maxlev)
%    - set the maximum signal level, MAXLEV, in counts, of the grayscale
%    image.  For example, for a 16-bit greyscale image, the maximum number
%    of counts is 2^16 = 65,536: this is the default value.
%
%    [...] = 
%    ppfakecalim(xpix,ypix,pixsize,spotdist,spotsize,imtilt,maxlev,noiselev)
%    - specify a noise level NOISELEV, in percent, for each light spot:
%    this adds random noise, up to +/- NOISELEV percent of each pixel, to
%    the light from each spot.  The default value is 0: no noise added.
%
%    [...] =
%    ppfakecalim(xpix,ypix,pixsize,spotdist,spotsize,imtilt,maxlev,noiselev,bglev)  
%    - add a constant level of background noise to the image of BGLEV
%    counts: this is equivalent to reducing the image signal level if the
%    value is constant (see BGNOISE below).  Default value is 0: no
%    background noise.
%
%    [...] =
%    ppfakecalim(xpix,ypix,pixsize,spotdist,spotsize,imtilt,maxlev,noiselev,bglev,bgnoise) 
%    - add Gaussian fluctuations to the level of background noise by
%    BGNOISE counts: this corresponds to 1 sigma of a Gaussian
%    distribution.  For example, for the PCO 2000 camera with a noise floor
%    of 100 counts which varied with a standard deviation of 10 counts,
%    BGLEV = 100 and BGNOISE = 10.

% File is the property of Simon Jolly, Imperial College London.  File
% released under the GNU public license.
%
% File history
%
%   04-Oct-2010 S. Jolly
%       Creation of fake pepperpot calibration image for new 4 light spot
%       calibration system.
%
%======================================================================


    if nargin < 3
        error('Need at least 3 input variables: [calim] = ppfakecalim(xpix,ypix,pixsize)') ;
    elseif nargin > 10
        error('Can only specify 10 input variables: [calim] = ppfakecalim(xpix,ypix,pixsize,spotdist,spotsize,imtilt,maxlev,noiselev,bglev,bgnoise)') ;
    elseif nargout > 1
        error('Can only specify 1 output variable: [calim] = ppfakecalim(xpix,ypix,pixsize)') ;
    end

    if nargin < 10 || isempty(bgnoise)
        bgnoise = 0 ;
    end

    if nargin < 9 || isempty(bglev)
        bglev = 0 ;
    end

    if nargin < 8 || isempty(noiselev)
        noiselev = 0 ;
    end

    if nargin < 7 || isempty(maxlev)
        maxlev = 2^16 ;
    end

    if nargin < 6 || isempty(imtilt)
        imtilt = 0 ;
    end

    if nargin < 5 || isempty(spotsize)
        spotsize = 500e-3 ;
    end

    if nargin < 4 || isempty(spotdist)
        spotdist = 140 ;
    end

%    Calculate central spot positions

    calim = zeros(ypix,xpix) ;
    midpoint = [( (ypix/2) - 0.5 ) ( (xpix/2) - 0.5 )] ;

    hspotdist = spotdist./2 ;

    topmm = [-(hspotdist.*sin(imtilt)) -(hspotdist.*cos(imtilt))] ;
    botmm = [(hspotdist.*sin(imtilt)) (hspotdist.*cos(imtilt))] ;
    leftmm = [-(hspotdist.*cos(imtilt)) (hspotdist.*sin(imtilt))] ;
    rightmm = [(hspotdist.*cos(imtilt)) -(hspotdist.*sin(imtilt))] ;

    toppixoff = topmm.*pixsize ;
    botpixoff = botmm.*pixsize ;
    leftpixoff = leftmm.*pixsize ;
    rightpixoff = rightmm.*pixsize ;

    toppix = midpoint + toppixoff ;
    botpix = midpoint + botpixoff ;
    leftpix = midpoint + leftpixoff ;
    rightpix = midpoint + rightpixoff ;

    if toppix(1) < 1 || toppix(1) > ypix || toppix(2) < 1 || toppix(2) > xpix
        warning('Centre of top spot falls outside the image') ;
    end

    if botpix(1) < 1 || botpix(1) > ypix || botpix(2) < 1 || botpix(2) > xpix
        warning('Centre of bottom spot falls outside the image') ;
    end

    if leftpix(1) < 1 || leftpix(1) > ypix || leftpix(2) < 1 || leftpix(2) > xpix
        warning('Centre of left spot falls outside the image') ;
    end

    if rightpix(1) < 1 || rightpix(1) > ypix || rightpix(2) < 1 || rightpix(2) > xpix
        warning('Centre of right spot falls outside the image') ;
    end

%    Create spot images

    topim = calim ;
    botim = calim ;
    leftim = calim ;
    rightim = calim ;

    spotpix = spotsize.*pixsize ;

    [topim] = gausspixspot([xpix ypix],spotpix,[],[],toppixoff) ;
    [botim] = gausspixspot([xpix ypix],spotpix,[],[],botpixoff) ;
    [leftim] = gausspixspot([xpix ypix],spotpix,[],[],leftpixoff) ;
    [rightim] = gausspixspot([xpix ypix],spotpix,[],[],rightpixoff) ;

%    calim = calim + topim + botim + leftim + rightim ;
    calim = calim + topim ;
    calim = calim + botim ;
%    calim = calim + leftim ;
    calim = calim + rightim ;

%    Add noise

    if noiselev > 0

        calsize = size(calim) ;
        spotnoise = ( (rand(calsize).*2) - 1 ).*(noiselev./100) + 1 ;
        calim = calim.*spotnoise ;

    end
    
    maxim = max(calim(:)) ;
    calim = round(calim.*(maxlev./maxim)) ;
    gaussnoise = 0 ;

    if bglev > 0

        if bgnoise > 0
            gaussnoise = randn(calsize).*bgnoise ;
        end
        
        bgim = bglev + gaussnoise ;

        bginds = bgim > calim ;
        calim(bginds) = bgim(bginds) ;

        toobiglev = calim > maxlev ;
        calim(toobiglev) = maxlev ;

    end

    return






