HGEXAMPLE Hermite Gaussian Examples

by Chuck DiMarzio Northeastern University December 2008

!! This file may be copied, used, or modified for educational and
!! research purposes provided that this header information is not
!! removed or altered, and provided that the book is cited in
!! publications, as DiMarzio, Charles A., Optics for Engineers,
!! CRC Press, Boca Raton, FL, 2011.
!! http://www.crcpress.com
!! Other distribution is prohibited without permission.

Contents

First, let's do a circular aperture

gamma=1;
%
%
radius=2;
xmax=5;
dx=0.02;
ymax=xmax;
dy=dx;
mmax=20;nmax=20;
xaxis=(-xmax:dx:xmax);
yaxis=(-ymax:dy:ymax);
[x,y]=meshgrid(xaxis,yaxis);
rsq=x.^2+y.^2;
u=(rsq<radius^2)*1/pi/radius^2;
% Show source pattern
fig26=figure;imagesc(xaxis,yaxis,abs(u).^2);
axis image;colormap(flipud(gray));colorbar;
%print('-deps',[epspath,'9-26-hgex.eps']);
%
%  First, find the best w
%
wlist=(0.7:0.02:1.2)*radius;
c00=zeros(1,length(wlist));
for nw=1:length(wlist);
  h=hgm(0,0,x,y,wlist(nw),wlist(nw));
  c00(nw)=sum(sum(u.*conj(h)))*dx*dy;
end;
% Show optimization of w
fig27=figure;plot(wlist,c00,'k');grid on;
xlabel('w, Gaussian radius');ylabel('C_{00}, Coefficient');

w=wlist(find(c00==max(c00)));
c=zeros(mmax,nmax);
for m=0:2:mmax;
  for n=0:2:nmax;
    h=hgm(m,n,x,y,w,w);
    c(m+1,n+1)=sum(sum(u.*conj(h)))*dx*dy; % m+1 because subscripts
                                                % start at 1. c(1,1) is c00
  end;
end;

Show the coefficients

fig28=figure;imagesc(20*log10(abs(c)));axis image;colormap(flipud(gray));
caxis([-50,0]);colorbar;

Approximation

yslice=floor(length(xaxis)/2);
ua=zeros(size(x));
usave=zeros(size(x,1),mmax+1);
for m=0:2:mmax;
  for n=0:2:nmax;
    h=hgm(m,n,x,y,w,w);
    ua=ua+c(m+1,n+1)*h; % remember the index starts at 1
  end;
    usave(:,m+1)=ua(:,yslice);
end;

Image the original source

fig29=figure;imagesc(xaxis,yaxis,abs(ua).^2);
axis image;colormap(flipud(gray));colorbar;

Plot alternate even orders (1:4:end)

fig30=figure;plot(xaxis,20*log10(abs(usave(:,1:4:end))),'k');grid on;
xlabel('x, Radial Distance');ylabel('Irradiance, dB');
moose=axis;axis([moose(1:2),moose(4)+[-60,0]]);

Far field;

newxmax=10;
newdx=0.1;
newymax=newxmax;
newdy=newdx;
newxaxis=(-newxmax:newdx:newxmax);
newyaxis=(-newymax:newdy:newymax);
[newx,newy]=meshgrid(newxaxis,newyaxis);

yslice=floor(length(newxaxis)/2);
uaf=zeros(size(newx));
usavef=zeros(size(newx,1),mmax+1);
for m=0:2:mmax;
  for n=0:2:nmax;
    h=hgm(m,n,newx,newy,w,w);
    uaf=uaf+c(m+1,n+1)*exp(1i*pi/2*(1+m+n))*h; % remember the index starts at 1
  end;
    usavef(:,m+1)=uaf(:,yslice);
end;

plot diffraction pattern in db.

fig31=figure;imagesc(newxaxis,newyaxis,20*log10(abs(uaf)));
caxis([-50,-20]);
axis image;colormap(flipud(gray));colorbar;

Plot alternate even orders (1:4:end)

fig32=figure;plot(newxaxis*pi/w^2,20*log10(abs(usavef(:,1:4:end))),'k');
grid on;
moose=axis;axis([moose(1:2),moose(4)+[-60,0]]);
xlabel('x, Radial Distance');ylabel('Irradiance, dB');