%%% Script for plotting Transversal mode Intensity distributions %%

clc, clear all, close all

%% Parameters %%

lam = 632*10^-9;       %Wavelength [m]
w0 = 400*10^-6;        %Beam waist [m]
k = 2*pi/lam;          %Wave number[m]
M2 = 1;                %Beam Quality factor
z = 0.001;             %The beam 1mm away from waist w0

phi = atan(lam*z/(pi*(w0^2)));  % Phase

X = [-0.001:0.00001:0.001];     % Specify x and y intervals.
Y = [-0.001:0.00001:0.001];

R = z.*(1+((pi*w0^2)./(M2*lam.*z)).^2);         % Equation 4.7.17 b
w = w0*sqrt(1+((M2*lam.*z)./(pi*(w0^2))).^2);   % Equation 4.7.17 a
[x,y] = meshgrid(X,Y);                          % Create a x-y grid

%%% Mode E field (U) given by eq 4.7.31

Hl0 = 1;
Hm0 = 1;                % Hermite Polynomials 
Hl1 = 2*sqrt(2)*(y./w);
Hm1 = 2*sqrt(2)*(x./w);
Hl2 = 8*((y.^2)/(w^2)) - 2;
Hm2 = 8*((x.^2)/(w^2)) - 2;
Hm3 = (8*(sqrt(2).*x/w).^3) - 12*(sqrt(2).*x/w);

U00 = (w0./w).*(Hl0)*(Hm0).*exp(-(x.^2 + y.^2)./(w.^2)).*exp(-i*k*(x.^2 + y.^2)./(2.*R) +i*(1)*phi);
U01 = (w0./w).*(Hl1)*(Hm0).*exp(-(x.^2 + y.^2)./(w.^2)).*exp(-i*k*(x.^2 + y.^2)./(2.*R) +i*(2)*phi);
U11 = (w0./w).*(Hl1)*(Hm1).*exp(-(x.^2 + y.^2)./(w.^2)).*exp(-i*k*(x.^2 + y.^2)./(2.*R) +i*(3)*phi);
U21 = (w0./w).*(Hl1)*(Hm2).*exp(-(x.^2 + y.^2)./(w.^2)).*exp(-i*k*(x.^2 + y.^2)./(2.*R) +i*(4)*phi);
U30 = (w0./w).*(Hl0)*(Hm3).*exp(-(x.^2 + y.^2)./(w.^2)).*exp(-i*k*(x.^2 + y.^2)./(2.*R) +i*(4)*phi);

I00 = abs(U00);             % Intensity is the magnitude of the E-field.
I01 = abs(U01);
I11 = abs(U11);
I21 = abs(U21);
I30 = abs(U30);

figure;
mesh(x,y,I00);
title('TEM00 Intensity Distribution, w0 = 0.4mm')
xlabel('X [m]')
ylabel('Y [m]')
view(0,90)

figure;
mesh(x,y,I01);
view(0,90)
title('TEM01 Intensity Distribution')
xlabel('X [m]')
ylabel('Y [m]')

figure;
mesh(x,y,I11);
view(0,90)
title('TEM11 Intensity Distribution')
xlabel('X [m]')
ylabel('Y [m]')

figure;
mesh(x,y,I21);
view(0,90)
title('TEM21 Intensity Distribution')
xlabel('X [m]')
ylabel('Y [m]')

figure;
mesh(x,y,I30);
view(0,90)
title('TEM30 Intensity Distribution')
xlabel('X [m]')
ylabel('Y [m]')





