diff --git a/Estimations/OpticalAccordionLattice.m b/Estimations/OpticalAccordionLattice.m
index e432ba2..8d1a811 100644
--- a/Estimations/OpticalAccordionLattice.m
+++ b/Estimations/OpticalAccordionLattice.m
@@ -30,6 +30,7 @@ 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)];
@@ -43,44 +44,63 @@ RotatedCoords2 = @(theta) R(-theta) * [x; y; z];
 
 %% Define E fields
 
-k1vec = k1(theta);
-coords = [x; y; z];
-rot1 = RotatedCoords1(theta);
-rot2 = RotatedCoords2(theta);
+% Generic orthogonal basis given a wavevector
+getPolarizationBasis = @(kvec) deal( ...
+    simplify( cross(kvec, [1;0;0]) ), ...
+    simplify( cross(kvec, cross(kvec, [1;0;0])) ) ...
+);
 
-% Polarization vector
-e_pol = cos(gamma)*[0; 0; 1] + sin(gamma)*[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 .* exp(1i * (k1(theta).' * coords)) * ...
+     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 .* exp(1i * (k2(theta).' * coords)) * ...
+     e_pol_2 .* exp(1i * (k2(theta).' * coords)) * ...
      exp(-(rot2(1)^2 / wx2^2) - (rot2(3)^2 / wz2^2));
 
-Efield = simplify(E1 + E2);
+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)
 
-%% Intensity expression
-Intensity = simplify(1/2 * real(conj(Efield) .* Efield));  % 3-component
-
-%% ================ Plot lattice =================== %%
+%% ================ Plot lattice in Y-Z plane =================== %%
 
 % Define parameters
-theta_val      = 10 * pi / 180;  % 10 degrees in radians
-gamma_val      = pi/2.0;              % tilt of linear polarization
+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(Intensity(3), x, 0));
+Iplane_z       = simplify(subs(IntensityVector(3), x, 0));
 
 % Convert to function
-Iplane_func    = matlabFunction(Iplane_z, 'Vars', {y, z, theta, wo, lambda, P, gamma});
+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, 300));
+[ygrid, zgrid] = meshgrid(linspace(-1000, 1000, 500), linspace(-100, 100, 500));
 
 % Evaluate intensity
-Ivals          = Iplane_func(ygrid, zgrid, theta_val, wo_val, lambda_val, P_val, gamma_val);
+Ivals          = Iplane_z_func(ygrid, zgrid, theta_val, wo_val, lambda_val, P_val, gamma_val);
 
 % Normalization
 Ivals          = Ivals / max(Ivals(:));  
@@ -88,19 +108,17 @@ Ivals          = Ivals / max(Ivals(:));
 % Plotting
 figure(1)
 clf
-set(gcf,'Position',[50 50 950 750])
+set(gcf,'Position',[50 50 900 700])
 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;
+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 =================== %%
+% ================ 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
@@ -121,21 +139,132 @@ 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);
+set(gcf,'Position',[50 50 900 700])
+plot(yvec, Iprop_cut, 'LineWidth', 2);
 title('Profile at x=0, z=0');
 xlabel('y [\mum]');
-ylabel('Depth');
+ylabel('Intensity');
 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);
+set(gcf,'Position',[50 50 900 700])
+plot(zvec, Ivert_cut, 'LineWidth', 2);
 title('Profile at x=0, y=0');
 xlabel('z [\mum]');
-ylabel('Depth');
+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);