%% a function to fit a thermal cloud 2-D
function fp = fitgaussian2D(m,cx,cy,tol);
%% m = image
%% tol = fitting tolerance

%Here we prepare a set of options to do our fits
options = optimset('Display','off','TolFun',tol,'LargeScale','off');
theta = 0;
background=0;

%We determine the size of an image
[sizey sizex] = size(m);
% guessing the widths
[tx,ty,sx,sy] = centerofmass(m);
sx = 2.*sx;
sy = 2.*sy;
%peak value
pOD = max(max(m));

%This is to avoid any possible access to some sells outside of the matrix
tmx = round(cy);
if isnan(tmx) > 0
    tmx = sizey/2;
end;
if round(cy) > sizey
    tmx = sizey;
end;
if round(cy) < 1
    tmx = 1;
end;
mx = m(tmx,:);


%Here we do the 1D horizontal fit
%This is to give better initial parameters to our 2D fit
x1D = 1:sizex;
%Initiall parameters
ip1D = [cx,sx,pOD,background];
fp1D = fminunc(@fitGaussian1D,ip1D,options,mx,x1D);

cx = fp1D(1);
sx = fp1D(2);
PeakOD = fp1D(3);
background = fp1D(4);

tmy = round(cx);
if isnan(tmy) > 0
    tmy = sizex/2;
end;
if round(cx) > sizex
    tmy = sizex;
end;
if round(cx) < 1
    tmy = 1;
end;
my = m(:,tmy)';

y1D = 1:sizey;
ip1D = [cy,sy,PeakOD,background];
%1D vertical fit
fp1D = fminunc(@fitGaussian1D,ip1D,options,my,y1D);

cy = fp1D(1);
sy = fp1D(2);
PeakOD = fp1D(3);
background = fp1D(4);
[X,Y] = meshgrid(1:sizex,1:sizey);

%Lower and upper limits for the 2D fit
LB = [1,1,1,1,-inf,-1.58,0];
UB = [sizex,sizey,sizex*10,sizey*10,inf,1.58,inf];

% Here we do the 2D fit
options = optimset('lsqnonlin');
options.Display = 'iter';
options.TolFun = tol;
initpar = [cx,cy,sx,sy,PeakOD,theta,background];
%this is to avoid errors when doing fits on the image full of zeros
if isnan(sum(sum(Gaussian2Dff(initpar,X,Y,m)))) < 1
    % Here is the actual fit
    fp = lsqnonlin( @(pt) Gaussian2Dff(pt,X,Y,m),initpar,LB,UB,options);
else
    fp = [0,0,0,0,0,0,0];
end;



% 2D Gauss function
function [z] = Gaussian2Dff(p,X,Y,m);

cx = p(1);
cy = p(2);
wx = p(3);
wy = p(4);
amp = p(5);
theta = p(6);
background=p(7);
Xn = (X-cx)*cos(theta) - (Y-cy)*sin(theta);
Yn = (X-cx)*sin(theta) + (Y-cy)*cos(theta);

z = amp*(exp(-2*((Xn).^2./(wx^2)+(Yn).^2./(wy^2)))) + background - m;