142 lines
4.0 KiB
Matlab
142 lines
4.0 KiB
Matlab
%% Physical Constants
|
|
|
|
PlanckConstant = 6.62606957e-34;
|
|
PlanckConstantReduced = PlanckConstant / (2 * pi);
|
|
FineStructureConstant = 7.2973525698e-3;
|
|
ElectronMass = 9.10938291e-31;
|
|
GravitationalConstant = 6.67384e-11;
|
|
ProtonMass = 1.672621777e-27;
|
|
AtomicMassUnit = 1.66053878283e-27;
|
|
BohrRadius = 0.52917721092e-10;
|
|
BohrMagneton = 927.400968e-26;
|
|
BoltzmannConstant = 1.3806488e-23;
|
|
StandardGravityAcceleration = 9.80665;
|
|
SpeedOfLight = 299792458;
|
|
StefanBoltzmannConstant = 5.670373e-8;
|
|
ElectronCharge = 1.602176565e-19;
|
|
VacuumPermeability = 4 * pi * 1e-7;
|
|
DielectricConstant = 1 / (SpeedOfLight^2 * VacuumPermeability);
|
|
ElectronGyromagneticFactor = -2.00231930436153;
|
|
AvogadroConstant = 6.02214129e23;
|
|
|
|
%% Parameters
|
|
syms x y z theta lambda P wo wx1 wz1 wx2 wz2 I gamma real
|
|
|
|
% Define constants
|
|
lambda_val = 0.532; % µm
|
|
P_val = 1;
|
|
wo_val = 100;
|
|
|
|
% Set beam waists equal for simplicity
|
|
wx1 = wo; wz1 = wo;
|
|
wx2 = wo; wz2 = wo;
|
|
%% Rotation matrix and k-vectors
|
|
% Rotation matrix
|
|
R = @(theta) [1 0 0; 0 cos(theta) -sin(theta); 0 sin(theta) cos(theta)];
|
|
|
|
% Define rotated coordinates and k-vectors
|
|
k1 = @(theta) R(theta) * [0; 1; 0];
|
|
k2 = @(theta) R(-theta) * [0; 1; 0];
|
|
|
|
RotatedCoords1 = @(theta) R(theta) * [x; y; z];
|
|
RotatedCoords2 = @(theta) R(-theta) * [x; y; z];
|
|
|
|
%% Define E fields
|
|
|
|
k1vec = k1(theta);
|
|
coords = [x; y; z];
|
|
rot1 = RotatedCoords1(theta);
|
|
rot2 = RotatedCoords2(theta);
|
|
|
|
% Polarization vector
|
|
e_pol = cos(gamma)*[0; 0; 1] + sin(gamma)*[1; 0; 0];
|
|
|
|
E1 = sqrt((2 * P) / (pi * wx1 * wz1)) * ...
|
|
e_pol .* exp(1i * (k1(theta).' * coords)) * ...
|
|
exp(-(rot1(1)^2 / wx1^2) - (rot1(3)^2 / wz1^2));
|
|
|
|
E2 = sqrt((2 * P) / (pi * wx2 * wz2)) * ...
|
|
e_pol .* exp(1i * (k2(theta).' * coords)) * ...
|
|
exp(-(rot2(1)^2 / wx2^2) - (rot2(3)^2 / wz2^2));
|
|
|
|
Efield = simplify(E1 + E2);
|
|
|
|
%% Intensity expression
|
|
Intensity = simplify(1/2 * real(conj(Efield) .* Efield)); % 3-component
|
|
|
|
%% ================ Plot lattice =================== %%
|
|
|
|
% Define parameters
|
|
theta_val = 10 * pi / 180; % 10 degrees in radians
|
|
gamma_val = pi/2.0; % tilt of linear polarization
|
|
|
|
% Extract z-component of intensity at x = 0
|
|
Iplane_z = simplify(subs(Intensity(3), x, 0));
|
|
|
|
% Convert to function
|
|
Iplane_func = matlabFunction(Iplane_z, 'Vars', {y, z, theta, wo, lambda, P, gamma});
|
|
|
|
% Grid for y and z
|
|
[ygrid, zgrid] = meshgrid(linspace(-1000, 1000, 500), linspace(-100, 100, 300));
|
|
|
|
% Evaluate intensity
|
|
Ivals = Iplane_func(ygrid, zgrid, theta_val, wo_val, lambda_val, P_val, gamma_val);
|
|
|
|
% Normalization
|
|
Ivals = Ivals / max(Ivals(:));
|
|
|
|
% Plotting
|
|
figure(1)
|
|
clf
|
|
set(gcf,'Position',[50 50 950 750])
|
|
contourf(ygrid, zgrid, Ivals, 200, 'LineColor', 'none');
|
|
colormap('turbo');
|
|
colorbar;
|
|
% Preserve physical aspect ratio
|
|
pbaspect([1 1 1]); % Set plot box aspect ratio to 1:1:1
|
|
axis tight;
|
|
xlabel('y [µm]', 'FontSize', 12);
|
|
ylabel('z [µm]', 'FontSize', 12);
|
|
title(['I_{plane}(y, z) at x = 0, \theta = ' num2str(rad2deg(theta_val)) '^\circ'], 'FontSize', 14);
|
|
set(gca, 'FontSize', 12, 'Box', 'on');
|
|
|
|
%% ================ Plot Potentials of lattice =================== %%
|
|
|
|
% Find indices closest to zero in y and z grids:
|
|
[~, idx_y0] = min(abs(ygrid(1,:))); % y=0 along columns
|
|
[~, idx_z0] = min(abs(zgrid(:,1))); % z=0 along rows
|
|
|
|
% Cut along y at z=0:
|
|
% z=0 corresponds to row idx_z0, extract entire column idx_z0 in y direction
|
|
Iprop_cut = Ivals(idx_z0, :); % 1D array vs y
|
|
|
|
% Cut along z at y=0:
|
|
% y=0 corresponds to column idx_y0, extract entire row idx_y0 in z direction
|
|
Ivert_cut = Ivals(:, idx_y0); % 1D array vs z
|
|
|
|
% Extract corresponding y and z vectors
|
|
yvec = ygrid(1, :);
|
|
zvec = zgrid(:, 1);
|
|
|
|
% Plot -Iprop/2 along y
|
|
figure(2);
|
|
clf
|
|
set(gcf,'Position',[50 50 950 750])
|
|
plot(yvec, -Iprop_cut/2, 'LineWidth', 2);
|
|
title('Profile at x=0, z=0');
|
|
xlabel('y [\mum]');
|
|
ylabel('Depth');
|
|
grid on;
|
|
set(gca, 'FontSize', 12, 'Box', 'on');
|
|
|
|
% Plot -Ivert/2 along z
|
|
figure(3);
|
|
clf
|
|
set(gcf,'Position',[50 50 950 750])
|
|
plot(zvec, -Ivert_cut/2, 'LineWidth', 2);
|
|
title('Profile at x=0, y=0');
|
|
xlabel('z [\mum]');
|
|
ylabel('Depth');
|
|
grid on;
|
|
set(gca, 'FontSize', 12, 'Box', 'on');
|