% Performs the convolution of a pattern with a beam shape
% All dimensions are in pixelx which corresponds approximatly 25nm
% The base shape function is "square"
% Other functions (chess, imaxis, slitpattern) are builf from the previous

function ebeamprox();
    tic
    
    ims=2^8; %image size in pixel

    %define the beam shape
    p = spotshape(4,.5,90,ims); 
    figure(4);dispim(p);title('Beam shape - color levels');
%    figure(4);dispcont(p);title('Beam shape - contour levels');
    
    %d1 = 1;
    %d2=1.1; d3=1.1;
    %d2=d1; d3=d1*.8;
    %d2=d1; d3=d1*1.1;
    %dose=[[d3 d1 d3]; [d2 d1 d2]; [d1 d1 d1]; [d2 d1 d2]; [d3 d1 d3]];
    %f=chess(40,40,3,5,2,2,dose,ims);

    f=slitpattern(40,80,150,20,120,ims);
    
    %f=f+imaxis(100,ims); %draw an axis in the middle of the image

    %draw the initial pattern
    figure(1);dispim(f);title('Intial pattern');

    %draw the convolution result
    r = convpattern(f,p,ims);
    figure(2);dispim(r);title('Convolution - color levels');
    figure(3);dispcont(r);title('Convolution - contour levels');
    
    toc; beep;beep;beep;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set of functions used to define the patterns    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
function f = chess(w,h,npx,npy,nx,ny,dose,ims)
    ww = ceil(w/npx);
    hh = ceil(h/npy);
    ox=0; oy=0;
    f=zeros(ims,ims);
    for k=0:ny-1
        for l=0:nx-1
            ox=-nx*w+l*2*w;
            oy=-ny*h+k*2*h;
            for j=0:npy-1
                for i=0:npx-1
                    f=f + square(dose(j+1,i+1),ox+i*ww,oy+j*hh,ww,hh,ims,ims);
                    f=f + square(dose(j+1,i+1),ox+w+i*ww,oy+h+j*hh,ww,hh,ims,ims);
                end
            end
        end
    end

function f = slitpattern(w,sw,sl,fw,fl,ims)
    d2=0.55; d3=0.9; d4=0.45; d5=1.3; d6=1.05; d7=1.05; d8=0.45;
    oy=(sl-w)/2;
    %left side
    f=square(d7,-sw/2-w,oy,w/2,w/2,ims,ims);
    f=f + square(d2,-sw/2-w/2,oy,w/2,w/2,ims,ims);
    f=f + square(d7,-sw/2-w/2,w/2+oy,w/2,w/2,ims,ims);
    f=f + square(d3,-sw/2-w,w/2+oy,w/2,w/2,ims,ims);
    f=f + square(3,-sw/2-w,7*w/8+oy,w/8,w/8,ims,ims);
    f=f + square(1,-sw/2,oy,(sw-fw)/2,w,ims,ims);
    f=f + square(1,-sw/2-w,-sl+oy,w,sl,ims,ims);
    %right side
    f=f + square(d2,sw/2,oy,w/2,w/2,ims,ims);
    f=f + square(1,sw/2+w/2,oy,w/2,w/2,ims,ims);
    f=f + square(1,sw/2+w/2,w/2+oy,w/2,w/2,ims,ims);
    f=f + square(1,sw/2,w/2+oy,w/2,w/2,ims,ims);
    f=f + square(1,fw/2,oy,(sw-fw)/2,w,ims,ims);
    f=f + square(1,sw/2,-sl+oy,w,sl,ims,ims);
    %finger
    f=f + square(d6,-fw/2,w/2+oy,fw,w/2,ims,ims);
    f=f + square(d8,-fw/2,oy,fw,w/2,ims,ims);
    f=f + square(d5,-fw/2,-1*fl/5+oy,fw,fl/5,ims,ims);
    f=f + square(d5,-fw/2,-2*fl/5+oy,fw,fl/5,ims,ims);
    f=f + square(d5,-fw/2,-3*fl/5+oy,fw,fl/5,ims,ims);
    f=f + square(d5,-fw/2,-4*fl/5+oy,fw,fl/5,ims,ims);
    f=f + square(d5*1.1,-fw/2,-5*fl/5+oy,fw,fl/5,ims,ims);

% function that draws a thin cross l long at 0,0
function f = imaxis(l,ims)
    f=square(10,-l/2,0,l,1,ims,ims);
    f=f + square(10,0,-l/2,1,l,ims,ims);
    
% creates a square of size lx,ly and of amplitude z 
% at position x,y (lower left corner)
% in center coordinate int(mx/2),int(my/2) of a matrix of size mx,my
function s = square(z, x, y, lx, ly, mx, my)
    xx=ceil(mx/2)+x+1;
    yy=ceil(my/2)-y;
    s=zeros(mx,my);
    s(yy-abs(ly)+1:yy,xx:xx+abs(lx)-1)=z;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set of imaging functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display a matrix image
function m = dispim(im)
    imshow(im/max(max(im)),'notruesize'); colormap(jet);

% display a matrix contour
function m = dispcont(im)
    ims=size(im); ims=ims(1:1);
    [C,h]=contour(im/max(max(im)),[0.25 0.5 0.75]);
    clabel(C,[0.25 0.75]);
    set(gca,'YDir','reverse')
    axis('square')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set of imaging functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = spotshape(psize, pastratio, pastangle, ims)
    deltax = 1/ims; 
    [xx,yy] = meshgrid(-0.5+deltax:deltax:0.5,-0.5+deltax:deltax:0.5);
    %psize = 75;
    %pastratio = .25 /2; % astigmatism ratio y/x
    %pastangle = -45; % astigmatism angle in degrees
    p = exp(-(psize/ims)^-2*(pastratio/2*xx.^2+(1-pastratio/2)*yy.^2));
    p = p/max(max(p));
    p = imrotate(p, pastangle,'crop');

function r = convpattern(f,p,ims)
    r = conv2(f,p);
    r = imcrop(r,[ims/2 ims/2 ims-1 ims-1]);
