hermite gaussian beam

%this is an example script showing a use for the hermiteh function
%this script creates the Hermite-Gaussian 4,4 mode and displays it
w0=1; %2x the standard deviation of the gaussian
x0=linspace(-4,4,501);
[X,Y]=meshgrid(x0,x0); %create a grid of width 4 units x 4 units with 501x501 mesh
E0=exp(-(X.^2+Y.^2)/w0^2); %create the gaussian
A=hermiteH(2,sqrt(2)*X/w0).*hermiteH(0,sqrt(2)*Y/w0).*E0; %create the 16,16 hermite gaussian
%display the absoute value of the result
a=abs(A);
a=a./max(max(a));
imagesc(abs(a).^2);
function a=hermiteH(n,x)
%generate the nth hermite polynomial of x
m=0:floor(n/2);
[q1,q2]=size(m);
s=ndims(x);
[g1,g2]=size(x);
X=repmat(x,[ones(1,s),q1, q2]);
m=permute(repmat(m,[g1,1,g2]),[1,3,2]);
a=factorial(n)*sum((-1).^m./(factorial(m).*factorial(n-2*m)).*(2*X).^(n-2*m),3);
end