function Gbeam(w0,Dz, z, y, d, Factor, Xmax, Fz, What, How)

Lambd=638e-9;
stpX=Xmax/100; x=-Xmax: stpX: Xmax; 

K=2*pi/Lambd; z0=K*(w0^2)/2; ZZmax=z0*Fz; stpZ=ZZmax/100; 
zz=-(ZZmax-stpZ/2) :stpZ : (ZZmax+stpZ/2);
RR=zz.*(1+(z0./zz).^2); ww=w0.*sqrt(1+(zz./z0).^2);

w1=w0*sqrt(1+(z/z0)^2); w2=w0*sqrt(1+((z+Dz)/z0)^2); 
R1=z*(1+(z0/z)^2);R2=z*(1+(z0/(z+Dz))^2);

[X,Y] = meshgrid(x);
E0=1; 
Efield1=E0*exp(-(X.^2+Y.^2)/(w1^2));
Efield2=Factor*E0*exp(-((X-d).^2+Y.^2)/(w2^2));
Power1=Efield1.^2; Power2=Efield2.^2;

deltFi=-K/2/R2*((X-d).^2+Y.^2)+K/2/R1*(X.^2+Y.^2);
Interf=Efield1.*Efield2.*cos(deltFi);
Power=Power1+Power2+2*Interf;


hold off;
if What==1 plot(zz,RR,'r-');grid on; %pause; 
end;
if What==2 plot(zz,ww,'r-');grid on; %pause; 
end;

if What==3 
	if How==2 meshc(X,Y,Power); view(-5,30); %pause; 
	else
plot(x,Power1(round((y+Xmax)/stpX),:),'m:'); hold on;
plot(x,Power2(round((y+Xmax)/stpX),:),'g:'); 
plot(x,Power(round((y+Xmax)/stpX),:),'r-');
	grid on;
	end;
end;
if What==4 
plot(x,Power(round((y+Xmax)/stpX),:),'c-',x,Interf(round((y+Xmax)/stpX),:),'r:');
	grid on;
end;