Calculation of angular spectrum of a Gaussian laser beam

GaussianBeamAngles.m

clear all;

lambda=650e-9;

%beam radius (given by full width half maximum, FWHM)
radius=0.5e-3;
FWHM=2*radius;


x=(-1:1/40:1)*5*FWHM;
y=x;

[X Y]=meshgrid(x,y);

% VX=X /((x(2)-x(1))*(max(x)-min(x)));

vx=x/(x(2)-x(1))/(max(x)-min(x));
VX=X /(x(2)-x(1))/(max(x)-min(x));

vy=y /(y(2)-y(1))/(max(y)-min(y));
VY=Y /(y(2)-y(1))/(max(y)-min(y));

%standard deviation
std_f=FWHM/(2*sqrt(2*log(2)));

f=exp(-(X.^2+Y.^2)/(2*std_f.^2));
%%
subplot(121)
imagesc(x*1e3,y*1e3,f.^2);

hold on;
th = 0:pi/50:2*pi;
r=radius;
xunit = r * cos(th) ;
yunit = r * sin(th) ;
plot(xunit*1e3, yunit*1e3,'g:','linewidth',2);
hold off;
axis image;

title(['Spatial intensity distribution' sprintf('\nWavelength %d nm', round(lambda*1e9))]);
xlabel('x [mm]');
ylabel('y [mm]');



subplot(122)
phix=asin(vx*lambda);
phiy=asin(vy*lambda);

F=fft2(f);

imagesc(phix*1e3,phiy*1e3,abs(fftshift(F).^2));

colormap hot;
axis image;

hold on;
th = 0:pi/50:2*pi;
std_F=(1/(2*pi*std_f));
FWHM_F=std_F*(2*sqrt(2*log(2)));
radius_F=FWHM_F/2;
r=asin(lambda*(radius_F));
xunit = r * cos(th) ;
yunit = r * sin(th) ;
plot(xunit*1e3, yunit*1e3,'g:','linewidth',2);
hold off;
xlabel('\theta_x [milliradians]');
ylabel('\theta_y [milliradians]');
axis image;

title(['Angular intensity distribution' sprintf('\n')]);