function [z] = Gaussian2D(p,X,Y);

% beam position (center of Gaussian)
cx = p(1);
cy = p(2); 
% waiste size 1/e in amplitude, 1/e^2 in intensity
wx = p(3);
wy = p(4);
amp = p(5); % beam amplitude or peak intensity
theta = p(6); % if beam is ellipticle then it is major axis angle
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;