You'll get my untested argument: Your free-widths parameters model include i and j and they might both be interpreted as (-1)^(1/2). I would argue that you should use fminsearch (or lsqnonlin) instead of fit. That would require you to explicitly write a function that calculates the sum-of-squared-residuals (or the residuals). That way you'll have explicit control over what goes on. I would also include a rotation-angle for your Gaussians. Something like this:
function err = towGaussianResiduals(pars,x,y,I)
Gaussians = pars(1);
for iP = 2:6:numel(pars)
Imax = pars(iP);
x0 = pars(iP+1);
sigmaX = pars(iP+2);
y0 = pars(iP+3);
sigmaY = pars(iP+4);
phi = pars(iP+5);
Gaussians = Gaussians + Imax*exp(-((+cos(phi)*(x-x0)+sin(phi)*(y-y0)).^2/sigmaX^2 + ...
(-sin(phi)*(x-x0)+cos(phi)*(y-y0)).^2/sigmaY^2));
end
err = sum(( I(:) - Gaussians(:)).^2);
end
That function you should be able to fit your parameters for with something like:
par0 = [0 I1 x1 sX1 y1 sY1 phi1 I2 x2 sX2 y2 sY2 phi2];
[x2D,y2D] = meshgrid(x,y);
pars = fminsearch(@(pars) twoGaussianResiduals(pars,x2D,y2D,I),par0);
You might have to check that the fitting has converged, check the help for fminsearch.
HTH
Best Answer