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 elliptical then it is major axis angle
background=p(7); 

% rotate squished beam appropriately
Xn = (X-cx)*cos(theta) - (Y-cy)*sin(theta);
Yn = (X-cx)*sin(theta) + (Y-cy)*cos(theta);

% finally we calculate Gaussian profile
z = amp*exp(-2*((Xn/wx).^2+(Yn/wy).^2)) + background;