%% 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 % Generic orthogonal basis given a wavevector getPolarizationBasis = @(kvec) deal( ... simplify( cross(kvec, [1;0;0]) ), ... simplify( cross(kvec, cross(kvec, [1;0;0])) ) ... ); % Get k-vectors k1vec = k1(theta); k2vec = k2(theta); % Get orthonormal basis for each polarization [e1_a, e1_b] = getPolarizationBasis(k1vec); [e2_a, e2_b] = getPolarizationBasis(k2vec); % Normalize them (optional if you want unit vectors) e1_a = simplify(e1_a / norm(e1_a)); e1_b = simplify(e1_b / norm(e1_b)); e2_a = simplify(e2_a / norm(e2_a)); e2_b = simplify(e2_b / norm(e2_b)); % Construct rotated polarization vectors using gamma e_pol_1 = cos(gamma) * e1_a + sin(gamma) * e1_b; e_pol_2 = cos(gamma) * e2_a + sin(gamma) * e2_b; coords = [x; y; z]; rot1 = RotatedCoords1(theta); rot2 = RotatedCoords2(theta); E1 = sqrt((2 * P) / (pi * wx1 * wz1)) * ... e_pol_1 .* 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_2 .* exp(1i * (k2(theta).' * coords)) * ... exp(-(rot2(1)^2 / wx2^2) - (rot2(3)^2 / wz2^2)); Efield = simplify(E1 + E2); IntensityVector = simplify(1/2 * real(conj(Efield) .* Efield)); % 3-component IntensityScalar = simplify(1/2 * sum(abs(Efield).^2)); % Scalar total intensity (with polarization interference) %% ================ Plot lattice in Y-Z plane =================== %% % Define parameters theta_val = deg2rad(10); % half-angle between beams gamma_val = deg2rad(90); % tilt of linear polarization % Extract z-component of intensity at x = 0 Iplane_z = simplify(subs(IntensityVector(3), x, 0)); % Convert to function Iplane_z_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, 500)); % Evaluate intensity Ivals = Iplane_z_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 900 700]) contourf(ygrid, zgrid, Ivals, 200, 'LineColor', 'none'); colormap('turbo'); colorbar; 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 1-D intensity profiles 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 900 700]) plot(yvec, Iprop_cut, 'LineWidth', 2); title('Profile at x=0, z=0'); xlabel('y [\mum]'); ylabel('Intensity'); grid on; set(gca, 'FontSize', 12, 'Box', 'on'); % Plot -Ivert/2 along z figure(3); clf set(gcf,'Position',[50 50 900 700]) plot(zvec, Ivert_cut, 'LineWidth', 2); title('Profile at x=0, y=0'); xlabel('z [\mum]'); ylabel('Intensity'); grid on; set(gca, 'FontSize', 12, 'Box', 'on'); %% Plot scalar intensity to see effect from change in polarization theta_val = deg2rad(10); % half-angle between beams Iplane_scalar = simplify(subs(IntensityScalar, x, 0)); Iplane_scalar_func = matlabFunction(Iplane_scalar, 'Vars', {y, z, theta, wo, lambda, P, gamma}); % Prepare figure figure(4); clf set(gcf,'Position',[50 50 900 700]) hold on % Define gamma values to compare gammas = linspace(0, pi/2, 10); % from 0 (p) to pi/2 (s) % Find y = 0 index [~, idx_y0] = min(abs(ygrid(1,:))); % column index % Colors for curves colors = turbo(length(gammas)); % Preallocate contrast array contrasts = zeros(size(gammas)); for i = 1:length(gammas) gamma_val = gammas(i); % Evaluate intensity Ivals = Iplane_scalar_func(ygrid, zgrid, theta_val, wo_val, lambda_val, P_val, gamma_val); Ivals = Ivals / max(Ivals(:)); % normalize % Extract z-cut at y = 0 Icut = Ivals(:, idx_y0); % Compute contrast Imax = max(Icut); Imin = min(Icut); contrasts(i) = (Imax - Imin) / (Imax + Imin); % Plot plot(zvec, Icut, 'LineWidth', 2, 'DisplayName', ... ['\gamma = ' num2str(rad2deg(gamma_val), '%0.1f') '^\circ'], ... 'Color', colors(i, :)); end % Finalize figure xlabel('z [\mum]', 'FontSize', 12); ylabel('Normalized Intensity', 'FontSize', 12); title('Intensity cut along z at y=0 for varying \gamma', 'FontSize', 14); legend('Location', 'best'); grid on box on set(gca, 'FontSize', 12); % --- New Figure for Contrast vs Gamma --- figure(5); clf set(gcf,'Position',[50 50 900 700]) plot(rad2deg(gammas), contrasts, '-o', 'LineWidth', 2); xlabel('\gamma [degrees]', 'FontSize', 14); ylabel('Fringe Contrast', 'FontSize', 14); title(['Contrast vs. Polarization Angle \gamma, \theta = ' num2str(rad2deg(theta_val)) '^\circ'], 'FontSize', 16); grid on box on set(gca, 'FontSize', 14); %% % --- Contrast vs Theta for fixed gamma = 0 and gamma = pi/2 --- % Define range of theta values (in radians) theta_vals = linspace(deg2rad(1), deg2rad(45), 50); % e.g., from 1° to 45° % Fixed gammas for contrast comparison gamma_fixed = [0, pi/2]; contrast_vs_theta = zeros(length(gamma_fixed), length(theta_vals)); for g_idx = 1:length(gamma_fixed) gamma_val = gamma_fixed(g_idx); for t_idx = 1:length(theta_vals) theta_val_local = theta_vals(t_idx); % Evaluate intensity on grid Ivals = Iplane_scalar_func(ygrid, zgrid, theta_val_local, wo_val, lambda_val, P_val, gamma_val); Ivals = Ivals / max(Ivals(:)); % normalize % Extract z-cut at y=0 Icut = Ivals(:, idx_y0); % Compute contrast Imax = max(Icut); Imin = min(Icut); contrast_vs_theta(g_idx, t_idx) = (Imax - Imin) / (Imax + Imin); end end % Plot contrast vs theta for the two fixed gammas figure(6); clf set(gcf,'Position',[50 50 900 700]) plot(rad2deg(theta_vals), contrast_vs_theta(1,:), '-o', 'LineWidth', 2, 'DisplayName', '\gamma = 0^\circ (p-pol)'); hold on plot(rad2deg(theta_vals), contrast_vs_theta(2,:), '-s', 'LineWidth', 2, 'DisplayName', '\gamma = 90^\circ (s-pol)'); xlabel('\theta [degrees]', 'FontSize', 14); ylabel('Fringe Contrast', 'FontSize', 14); title('Contrast vs. Half-angle \theta for fixed polarizations \gamma', 'FontSize', 16); legend('Location', 'southwest'); grid on box on set(gca, 'FontSize', 14);