%% Create 2D Gaussian with noise
[X Y] = meshgrid(-5:0.1:5,-5:0.1:5);
a = 2.1; bx = 2.3; by = -1.2; cx = 1.3; cy = 0.8;
sig = 0.1;
Z = a*exp(-(((X-bx)/cx).^2 + ((Y-by)/cy).^2)) + sig*randn(size(X));
%% Do "brute force" fit
% Initial values around which to search
a0 = max(max(Z));
[ix iy] = find(Z == a0);
bx0 = X(ix,iy); by0 = Y(ix,iy);
cx0 = 1; cy0 = 1;
astep = 0.05;
bstep = 0.05;
cstep = 0.05;
chisqmin = 1e20;
for a = a0-5*astep:astep:a0+5*astep
    for bx = bx0-5*bstep:bstep:bx0+5*bstep
        for by = by0-5*bstep:bstep:by0+5*bstep
            for cx = cx0-5*cstep:cstep:cx0+5*cstep
                for cy = cy0-5*cstep:cstep:cy0+5*cstep
                   ztest = a*exp(-(((X-bx)/cx).^2 + ((Y-by)/cy).^2));
                   chisq = sum(sum(((ztest-Z)/sig).^2));
                   if (chisq < chisqmin)
                      chisqmin = chisq;
                      params = [a,bx,by,cx,cy];
                   end
                end
            end
        end
    end
end
   