From f6981b92337f2a820c903bbbe37f3befded23eb2 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Fri, 18 Jul 2025 22:10:05 +0200 Subject: [PATCH] New routine to plot insitu images across runs, new script to calculate interference pattern, potential for an optical accordion lattice. --- .../SpectralAnalysisRoutines/plotImages.m | 16 +- .../plotPhaseDiagram.m | 767 ++++++++++++++++++ Estimations/OpticalAccordionLattice.m | 141 ++++ 3 files changed, 916 insertions(+), 8 deletions(-) create mode 100644 Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/plotPhaseDiagram.m create mode 100644 Estimations/OpticalAccordionLattice.m diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/plotImages.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/plotImages.m index 5a5ea3e..4866156 100644 --- a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/plotImages.m +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/plotImages.m @@ -4,16 +4,16 @@ groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axi "/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", ... "/images/Vertical_Axis_Camera/in_situ_absorption"]; -folderPath = "//DyLabNAS/Data/TwoDGas/2025/04/02/"; +folderPath = "//DyLabNAS/Data/TwoDGas/2025/07/16/"; -run = '0007'; +run = '0002'; folderPath = strcat(folderPath, run); cam = 5; angle = 0; -center = [1285, 2100]; +center = [1430, 2025]; span = [200, 200]; fraction = [0.1, 0.1]; @@ -25,13 +25,13 @@ ImagingMode = 'HighIntensity'; PulseDuration = 5e-6; % Plotting and saving -scan_parameter = 'rot_mag_fin_pol_angle'; +scan_parameter = 'evap_rot_mag_field'; scan_groups = 0:10:50; -savefileName = 'DropletsToStripes'; +savefileName = 'Droplets'; font = 'Bahnschrift'; % Flags -skipUnshuffling = false; +skipUnshuffling = true; %% ===== Load and compute OD image, rotate and extract ROI for analysis ===== % Get a list of all files in the folder with the desired file name pattern. @@ -175,14 +175,14 @@ for k = 1 : length(od_imgs) 'Interpreter', 'tex', 'Units', 'normalized', ... 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top'); end - + colorbarHandle = colorbar; ylabel(colorbarHandle, 'Optical Density', 'Rotation', -90, 'FontSize', 14, 'FontName', font); xlabel('x (\mum)', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font); ylabel('y (\mum)', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font); title('OD Image', 'FontSize', 16, 'FontWeight', 'bold', 'Interpreter', 'tex', 'FontName', font); - + drawnow; pause(0.5); end diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/plotPhaseDiagram.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/plotPhaseDiagram.m new file mode 100644 index 0000000..900bc5d --- /dev/null +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/plotPhaseDiagram.m @@ -0,0 +1,767 @@ +% === Parameters === +baseFolder = '//DyLabNAS/Data/TwoDGas/2025/04/'; + +dates = ["01", "02"]; +runs = { + ["0059", "0060", "0061"], + ["0007", "0008", "0009", "0010", "0011"] +}; + +scan_groups = 0:10:50; +scan_parameter = 'rot_mag_fin_pol_angle'; +cam = 5; + +% Image cropping and alignment +angle = 0; +center = [1285, 2100]; +span = [200, 200]; +fraction = [0.1, 0.1]; + +% Imaging and calibration parameters +pixel_size = 5.86e-6; % in meters +magnification = 23.94; +removeFringes = false; +ImagingMode = 'LowIntensity'; +PulseDuration = 5e-6; + +% Optional visualization / zooming +options.zoom_size = 50; + +% Optional flags or settings struct +skipUnshuffling = false; +skipPreprocessing = true; +skipMasking = true; +skipIntensityThresholding = true; +skipBinarization = true; + +groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis_Camera/in_situ_absorption", ... + "/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", ... + "/images/Vertical_Axis_Camera/in_situ_absorption"]; +%% + +allData = {}; % now a growing list of structs per B field +dataCounter = 1; + +for i = 1:length(dates) + dateStr = dates(i); + runList = runs{i}; + + for j = 1:length(runList) + folderPath = fullfile(baseFolder, dateStr, runList{j}); + filePattern = fullfile(folderPath, '*.h5'); + files = dir(filePattern); + refimages = zeros(span(1) + 1, span(2) + 1, length(files)); + absimages = zeros(span(1) + 1, span(2) + 1, length(files)); + + for k = 1 : length(files) + baseFileName = files(k).name; + fullFileName = fullfile(files(k).folder, baseFileName); + + fprintf(1, 'Now reading %s\n', fullFileName); + + atm_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/atoms")), angle)); + bkg_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle)); + dark_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle)); + + refimages(:,:,k) = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)'; + absimages(:,:,k) = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)'; + end + + % ===== Fringe removal ===== + + if removeFringes + optrefimages = removefringesInImage(absimages, refimages); + absimages_fringe_removed = absimages(:, :, :) - optrefimages(:, :, :); + + nimgs = size(absimages_fringe_removed,3); + od_imgs = cell(1, nimgs); + for k = 1:nimgs + od_imgs{k} = absimages_fringe_removed(:, :, k); + end + else + nimgs = size(absimages(:, :, :),3); + od_imgs = cell(1, nimgs); + for k = 1:nimgs + od_imgs{k} = absimages(:, :, k); + end + end + + %% ===== Get rotation angles ===== + scan_parameter_values = zeros(1, length(files)); + + % Get information about the '/globals' group + for k = 1 : length(files) + baseFileName = files(k).name; + fullFileName = fullfile(files(k).folder, baseFileName); + info = h5info(fullFileName, '/globals'); + for i = 1:length(info.Attributes) + if strcmp(info.Attributes(i).Name, scan_parameter) + if strcmp(scan_parameter, 'rot_mag_fin_pol_angle') + scan_parameter_values(k) = 180 - info.Attributes(i).Value; + else + scan_parameter_values(k) = info.Attributes(i).Value; + end + end + if strcmp(info.Attributes(i).Name, "rot_mag_field") + B = info.Attributes(i).Value; + end + end + end + + % ===== Unshuffle if necessary to do so ===== + + if ~skipUnshuffling + n_values = length(scan_groups); + n_total = length(scan_parameter_values); + + % Infer number of repetitions + n_reps = n_total / n_values; + + % Preallocate ordered arrays + ordered_scan_values = zeros(1, n_total); + ordered_od_imgs = cell(1, n_total); + + counter = 1; + + for rep = 1:n_reps + for val = scan_groups + % Find the next unused match for this val + idx = find(scan_parameter_values == val, 1, 'first'); + + % Assign and remove from list to avoid duplicates + ordered_scan_values(counter) = scan_parameter_values(idx); + ordered_od_imgs{counter} = od_imgs{idx}; + + % Mark as used by removing + scan_parameter_values(idx) = NaN; % NaN is safe since original values are 0:5:45 + od_imgs{idx} = []; % empty cell so it won't be matched again + + counter = counter + 1; + end + end + + % Now assign back + scan_parameter_values = ordered_scan_values; + od_imgs = ordered_od_imgs; + end + % === Reshape === + od_imgs_reshaped = reshape(od_imgs, [length(scan_groups), n_reps]); + + % === Store === + allData{dataCounter} = struct(... + 'B', B, ... + 'theta_vals', scan_groups, ... + 'od_imgs', od_imgs_reshaped ... + ); + dataCounter = dataCounter + 1; + end +end + +%% === % Plot PD - 1st rep of each θ per B-field === +[theta_vals, ~, idx] = unique(scan_parameter_values); +nB = numel(allData); +nTheta = numel(theta_vals); + +% Select every 2nd B-field index +idxToPlot = 1:2:nB; % indices 1, 3, 5, ... + +% Update number of B-fields to plot +nB_new = numel(idxToPlot); + +figure(101); clf; + +% Make the figure wider to fit the colorbar comfortably +set(gcf, 'Position', [100, 100, 1300, 800]); + +% Create tiled layout with some right padding to reserve space for colorbar +t = tiledlayout(nB_new, nTheta, 'TileSpacing', 'compact', 'Padding', 'compact'); + +font = 'Bahnschrift'; +allAxes = gobjects(nB_new, nTheta); + +for new_i = 1:nB_new + i = idxToPlot(new_i); % original index in allData + data = allData{i}; + for j = 1:nTheta + ax = nexttile((new_i-1)*nTheta + j); + allAxes(new_i,j) = ax; + + od = data(j).od_imgs; + imagesc(od, 'Parent', ax); + set(ax, 'YDir', 'normal'); + axis(ax, 'image'); + ax.XTick = []; + ax.YTick = []; + + colormap(ax, Colormaps.inferno()); + end +end + +% Use colorbar associated with the last image tile +cb = colorbar('Location', 'eastoutside'); +cb.Layout.Tile = 'east'; % Attach it to the layout edge +cb.FontName = font; +cb.FontSize = 18; +cb.Label.FontSize = 20; +cb.Label.Rotation = 90; +cb.Label.VerticalAlignment = 'bottom'; +cb.Label.HorizontalAlignment = 'center'; +cb.Direction = 'normal'; % Ensure ticks go bottom-to-top + + +% Add x and y tick labels along bottom and left +% Use bottom row for θ ticks +for j = 1:nTheta + ax = allAxes(end, j); + ax.XTick = size(od,2)/2; + ax.XTickLabel = sprintf('%d°', theta_vals(j)); + ax.XTickLabelRotation = 0; + ax.FontName = font; + ax.FontSize = 20; +end + +% Use first column for B ticks (only the plotted subset) +for new_i = 1:nB_new + i = idxToPlot(new_i); + ax = allAxes(new_i, 1); + ax.YTick = size(od,1)/2; + ax.YTickLabel = sprintf('%.2f G', allData{i}(1).B); + ax.FontName = font; + ax.FontSize = 20; +end + +%% Helper Functions +function [IMGFFT, IMGPR] = computeFourierTransform(I, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization) + % computeFourierSpectrum - Computes the 2D Fourier power spectrum + % of binarized and enhanced lattice image features, with optional central mask. + % + % Inputs: + % I - Grayscale or RGB image matrix + % + % Output: + % F_mag - 2D Fourier power spectrum (shifted) + + if ~skipPreprocessing + % Preprocessing: Denoise + filtered = imgaussfilt(I, 10); + IMGPR = I - filtered; % adjust sigma as needed + else + IMGPR = I; + end + + if ~skipMasking + [rows, cols] = size(IMGPR); + [X, Y] = meshgrid(1:cols, 1:rows); + % Elliptical mask parameters + cx = cols / 2; + cy = rows / 2; + + % Shifted coordinates + x = X - cx; + y = Y - cy; + + % Ellipse semi-axes + rx = 0.4 * cols; + ry = 0.2 * rows; + + % Rotation angle in degrees -> radians + theta_deg = 30; % Adjust as needed + theta = deg2rad(theta_deg); + + % Rotated ellipse equation + cos_t = cos(theta); + sin_t = sin(theta); + + x_rot = (x * cos_t + y * sin_t); + y_rot = (-x * sin_t + y * cos_t); + + ellipseMask = (x_rot.^2) / rx^2 + (y_rot.^2) / ry^2 <= 1; + + % Apply cutout mask + IMGPR = IMGPR .* ellipseMask; + end + + if ~skipIntensityThresholding + % Apply global intensity threshold mask + intensity_thresh = 0.20; + intensity_mask = IMGPR > intensity_thresh; + IMGPR = IMGPR .* intensity_mask; + end + + if ~skipBinarization + % Adaptive binarization and cleanup + IMGPR = imbinarize(IMGPR, 'adaptive', 'Sensitivity', 0.0); + IMGPR = imdilate(IMGPR, strel('disk', 2)); + IMGPR = imerode(IMGPR, strel('disk', 1)); + IMGPR = imfill(IMGPR, 'holes'); + F = fft2(double(IMGPR)); % Compute 2D Fourier Transform + IMGFFT = abs(fftshift(F))'; % Shift zero frequency to center + else + F = fft2(double(IMGPR)); % Compute 2D Fourier Transform + IMGFFT = abs(fftshift(F))'; % Shift zero frequency to center + end +end + +function [k_rho_vals, S_radial] = computeRadialSpectralDistribution(IMGFFT, kx, ky, thetamin, thetamax, num_bins) + % IMGFFT : 2D FFT image (fftshifted and cropped) + % kx, ky : 1D physical wavenumber axes [μm⁻¹] matching FFT size + % thetamin : Minimum angle (in radians) + % thetamax : Maximum angle (in radians) + % num_bins : Number of radial bins + + [KX, KY] = meshgrid(kx, ky); + K_rho = sqrt(KX.^2 + KY.^2); + Theta = atan2(KY, KX); + + if thetamin < thetamax + angle_mask = (Theta >= thetamin) & (Theta <= thetamax); + else + angle_mask = (Theta >= thetamin) | (Theta <= thetamax); + end + + power_spectrum = abs(IMGFFT).^2; + + r_min = min(K_rho(angle_mask)); + r_max = max(K_rho(angle_mask)); + r_edges = linspace(r_min, r_max, num_bins + 1); + k_rho_vals = 0.5 * (r_edges(1:end-1) + r_edges(2:end)); + S_radial = zeros(1, num_bins); + + for i = 1:num_bins + r_low = r_edges(i); + r_high = r_edges(i + 1); + radial_mask = (K_rho >= r_low) & (K_rho < r_high); + full_mask = radial_mask & angle_mask; + S_radial(i) = sum(power_spectrum(full_mask)); + end +end + +function [theta_vals, S_theta] = computeAngularSpectralDistribution(IMGFFT, r_min, r_max, num_bins, threshold, sigma, windowSize) + % Apply threshold to isolate strong peaks + IMGFFT(IMGFFT < threshold) = 0; + + % Prepare polar coordinates + [ny, nx] = size(IMGFFT); + [X, Y] = meshgrid(1:nx, 1:ny); + cx = ceil(nx/2); + cy = ceil(ny/2); + R = sqrt((X - cx).^2 + (Y - cy).^2); + Theta = atan2(Y - cy, X - cx); % range [-pi, pi] + + % Choose radial band + radial_mask = (R >= r_min) & (R <= r_max); + + % Initialize angular structure factor + S_theta = zeros(1, num_bins); + theta_vals = linspace(0, pi, num_bins); + + % Loop through angle bins + for i = 1:num_bins + angle_start = (i-1) * pi / num_bins; + angle_end = i * pi / num_bins; + angle_mask = (Theta >= angle_start & Theta < angle_end); + bin_mask = radial_mask & angle_mask; + fft_angle = IMGFFT .* bin_mask; + S_theta(i) = sum(sum(abs(fft_angle).^2)); + end + + % Smooth using either Gaussian or moving average + if exist('sigma', 'var') && ~isempty(sigma) + % Gaussian convolution + half_width = ceil(3 * sigma); + x = -half_width:half_width; + gauss_kernel = exp(-x.^2 / (2 * sigma^2)); + gauss_kernel = gauss_kernel / sum(gauss_kernel); + % Circular convolution + S_theta = conv([S_theta(end-half_width+1:end), S_theta, S_theta(1:half_width)], ... + gauss_kernel, 'same'); + S_theta = S_theta(half_width+1:end-half_width); + elseif exist('windowSize', 'var') && ~isempty(windowSize) + % Moving average via convolution (circular) + pad = floor(windowSize / 2); + kernel = ones(1, windowSize) / windowSize; + S_theta = conv([S_theta(end-pad+1:end), S_theta, S_theta(1:pad)], kernel, 'same'); + S_theta = S_theta(pad+1:end-pad); + end +end + +function contrast = computeRadialSpectralContrast(IMGFFT, r_min, r_max, threshold) + % Apply threshold to isolate strong peaks + IMGFFT(IMGFFT < threshold) = 0; + + % Prepare polar coordinates + [ny, nx] = size(IMGFFT); + [X, Y] = meshgrid(1:nx, 1:ny); + cx = ceil(nx/2); + cy = ceil(ny/2); + R = sqrt((X - cx).^2 + (Y - cy).^2); + + % Ring region (annulus) mask + ring_mask = (R >= r_min) & (R <= r_max); + + % Squared magnitude in the ring + ring_power = abs(IMGFFT).^2 .* ring_mask; + + % Maximum power in the ring + ring_max = max(ring_power(:)); + + % Power at the DC component + dc_power = abs(IMGFFT(cy, cx))^2; + + % Avoid division by zero + if dc_power == 0 + contrast = Inf; % or NaN or 0, depending on how you want to handle this + else + contrast = ring_max / dc_power; + end +end + +function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction) + % image must be a 2D numerical array + [dim1, dim2] = size(img); + + s1 = img(1:round(dim1 * y_fraction), 1:round(dim2 * x_fraction)); + s2 = img(1:round(dim1 * y_fraction), round(dim2 - dim2 * x_fraction):dim2); + s3 = img(round(dim1 - dim1 * y_fraction):dim1, 1:round(dim2 * x_fraction)); + s4 = img(round(dim1 - dim1 * y_fraction):dim1, round(dim2 - dim2 * x_fraction):dim2); + + ret = mean([mean(s1(:)), mean(s2(:)), mean(s3(:)), mean(s4(:))]); +end + +function ret = subtractBackgroundOffset(img, fraction) + % Remove the background from the image. + % :param dataArray: The image + % :type dataArray: xarray DataArray + % :param x_fraction: The fraction of the pixels used in x axis + % :type x_fraction: float + % :param y_fraction: The fraction of the pixels used in y axis + % :type y_fraction: float + % :return: The image after removing background + % :rtype: xarray DataArray + + x_fraction = fraction(1); + y_fraction = fraction(2); + offset = getBkgOffsetFromCorners(img, x_fraction, y_fraction); + ret = img - offset; +end + +function ret = cropODImage(img, center, span) + % Crop the image according to the region of interest (ROI). + % :param dataSet: The images + % :type dataSet: xarray DataArray or DataSet + % :param center: The center of region of interest (ROI) + % :type center: tuple + % :param span: The span of region of interest (ROI) + % :type span: tuple + % :return: The cropped images + % :rtype: xarray DataArray or DataSet + + x_start = floor(center(1) - span(1) / 2); + x_end = floor(center(1) + span(1) / 2); + y_start = floor(center(2) - span(2) / 2); + y_end = floor(center(2) + span(2) / 2); + + ret = img(y_start:y_end, x_start:x_end); +end + +function imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime) +%CALCULATEODIMAGE Calculates the optical density (OD) image for absorption imaging. +% +% imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime) +% +% Inputs: +% imageAtom - Image with atoms +% imageBackground - Image without atoms +% imageDark - Image without light +% mode - 'LowIntensity' (default) or 'HighIntensity' +% exposureTime - Required only for 'HighIntensity' [in seconds] +% +% Output: +% imageOD - Computed OD image +% + + arguments + imageAtom (:,:) {mustBeNumeric} + imageBackground (:,:) {mustBeNumeric} + imageDark (:,:) {mustBeNumeric} + mode char {mustBeMember(mode, {'LowIntensity', 'HighIntensity'})} = 'LowIntensity' + exposureTime double = NaN + end + + % Compute numerator and denominator + numerator = imageBackground - imageDark; + denominator = imageAtom - imageDark; + + % Avoid division by zero + numerator(numerator == 0) = 1; + denominator(denominator == 0) = 1; + + % Calculate OD based on mode + switch mode + case 'LowIntensity' + imageOD = -log(abs(denominator ./ numerator)); + + case 'HighIntensity' + if isnan(exposureTime) + error('Exposure time must be provided for HighIntensity mode.'); + end + imageOD = abs(denominator ./ numerator); + imageOD = -log(imageOD) + (numerator - denominator) ./ (7000 * (exposureTime / 5e-6)); + end + +end + +function drawODOverlays(x1, y1, x2, y2) + + % Parameters + tick_spacing = 10; % µm between ticks + tick_length = 2; % µm tick mark length + line_color = [0.5 0.5 0.5]; + tick_color = [0.5 0.5 0.5]; + font_size = 10; + + % Vector from start to end + dx = x2 - x1; + dy = y2 - y1; + L = sqrt(dx^2 + dy^2); + + % Unit direction vector along diagonal + ux = dx / L; + uy = dy / L; + + % Perpendicular unit vector for ticks + perp_ux = -uy; + perp_uy = ux; + + % Midpoint (center) + xc = (x1 + x2) / 2; + yc = (y1 + y2) / 2; + + % Number of positive and negative ticks + n_ticks = floor(L / (2 * tick_spacing)); + + % Draw main diagonal line + plot([x1 x2], [y1 y2], '--', 'Color', line_color, 'LineWidth', 1.2); + + for i = -n_ticks:n_ticks + d = i * tick_spacing; + xt = xc + d * ux; + yt = yc + d * uy; + + % Tick line endpoints + xt1 = xt - 0.5 * tick_length * perp_ux; + yt1 = yt - 0.5 * tick_length * perp_uy; + xt2 = xt + 0.5 * tick_length * perp_ux; + yt2 = yt + 0.5 * tick_length * perp_uy; + + % Draw tick + plot([xt1 xt2], [yt1 yt2], '--', 'Color', tick_color, 'LineWidth', 1); + + % Label: centered at tick, offset slightly along diagonal + if d ~= 0 + text(xt, yt, sprintf('%+d', d), ... + 'Color', tick_color, ... + 'FontSize', font_size, ... + 'HorizontalAlignment', 'center', ... + 'VerticalAlignment', 'bottom', ... + 'Rotation', atan2d(dy, dx)); + end + + end +end + +function drawPSOverlays(kx, ky, r_min, r_max) +% drawFFTOverlays - Draw overlays on existing FFT plot: +% - Radial lines every 30° +% - Annular highlight with white (upper half) and gray (lower half) circles between r_min and r_max +% - Horizontal white bands at ky=0 in annulus region +% - Scale ticks and labels every 1 μm⁻¹ along each radial line +% +% Inputs: +% kx, ky - reciprocal space vectors (μm⁻¹) +% r_min - inner annulus radius offset index (integer) +% r_max - outer annulus radius offset index (integer) +% +% Example: +% hold on; +% drawFFTOverlays(kx, ky, 10, 30); + + hold on + + % === Overlay Radial Lines + Scales === + [kx_grid, ky_grid] = meshgrid(kx, ky); + [~, kr_grid] = cart2pol(kx_grid, ky_grid); % kr_grid in μm⁻¹ + + max_kx = max(kx); + max_ky = max(ky); + + for angle = 0 : pi/6 : pi + x_line = [0, max_kx] * cos(angle); + y_line = [0, max_ky] * sin(angle); + + % Plot radial lines + plot(x_line, y_line, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2); + plot(x_line, -y_line, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2); + + % Draw scale ticks along positive radial line + drawTicksAlongLine(0, 0, x_line(2), y_line(2)); + + % Draw scale ticks along negative radial line (reflect y) + drawTicksAlongLine(0, 0, x_line(2), -y_line(2)); + end + + % === Overlay Annular Highlight: White (r_min to r_max), Gray elsewhere === + theta_full = linspace(0, 2*pi, 500); + + center_x = ceil(size(kr_grid, 2) / 2); + center_y = ceil(size(kr_grid, 1) / 2); + + k_min = kr_grid(center_y, center_x + r_min); + k_max = kr_grid(center_y, center_x + r_max); + + % Upper half: white dashed circles + x1_upper = k_min * cos(theta_full(theta_full <= pi)); + y1_upper = k_min * sin(theta_full(theta_full <= pi)); + x2_upper = k_max * cos(theta_full(theta_full <= pi)); + y2_upper = k_max * sin(theta_full(theta_full <= pi)); + plot(x1_upper, y1_upper, 'k--', 'LineWidth', 1.2); + plot(x2_upper, y2_upper, 'k--', 'LineWidth', 1.2); + + % Lower half: gray dashed circles + x1_lower = k_min * cos(theta_full(theta_full > pi)); + y1_lower = k_min * sin(theta_full(theta_full > pi)); + x2_lower = k_max * cos(theta_full(theta_full > pi)); + y2_lower = k_max * sin(theta_full(theta_full > pi)); + plot(x1_lower, y1_lower, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.0); + plot(x2_lower, y2_lower, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.0); + + % === Highlight horizontal band across k_y = 0 === + x_vals = kx; + xW1 = x_vals((x_vals >= -k_max) & (x_vals < -k_min)); + xW2 = x_vals((x_vals > k_min) & (x_vals <= k_max)); + + plot(xW1, zeros(size(xW1)), 'k--', 'LineWidth', 1.2); + plot(xW2, zeros(size(xW2)), 'k--', 'LineWidth', 1.2); + + hold off + + + % --- Nested helper function to draw ticks along a radial line --- + function drawTicksAlongLine(x_start, y_start, x_end, y_end) + % Tick parameters + tick_spacing = 1; % spacing between ticks in μm⁻¹ + tick_length = 0.05 * sqrt((x_end - x_start)^2 + (y_end - y_start)^2); % relative tick length + line_color = [0.5 0.5 0.5]; + tick_color = [0.5 0.5 0.5]; + font_size = 8; + + % Vector along the line + dx = x_end - x_start; + dy = y_end - y_start; + L = sqrt(dx^2 + dy^2); + ux = dx / L; + uy = dy / L; + + % Perpendicular vector for ticks + perp_ux = -uy; + perp_uy = ux; + + % Number of ticks (from 0 up to max length) + n_ticks = floor(L / tick_spacing); + + for i = 1:n_ticks + % Position of tick along the line + xt = x_start + i * tick_spacing * ux; + yt = y_start + i * tick_spacing * uy; + + % Tick endpoints + xt1 = xt - 0.5 * tick_length * perp_ux; + yt1 = yt - 0.5 * tick_length * perp_uy; + xt2 = xt + 0.5 * tick_length * perp_ux; + yt2 = yt + 0.5 * tick_length * perp_uy; + + % Draw tick + plot([xt1 xt2], [yt1 yt2], '-', 'Color', tick_color, 'LineWidth', 1); + + % Label with distance (integer) + text(xt, yt, sprintf('%d', i), ... + 'Color', tick_color, ... + 'FontSize', font_size, ... + 'HorizontalAlignment', 'center', ... + 'VerticalAlignment', 'bottom', ... + 'Rotation', atan2d(dy, dx)); + end + end + +end + +function [optrefimages] = removefringesInImage(absimages, refimages, bgmask) + % removefringesInImage - Fringe removal and noise reduction from absorption images. + % Creates an optimal reference image for each absorption image in a set as + % a linear combination of reference images, with coefficients chosen to + % minimize the least-squares residuals between each absorption image and + % the optimal reference image. The coefficients are obtained by solving a + % linear set of equations using matrix inverse by LU decomposition. + % + % Application of the algorithm is described in C. F. Ockeloen et al, Improved + % detection of small atom numbers through image processing, arXiv:1007.2136 (2010). + % + % Syntax: + % [optrefimages] = removefringesInImage(absimages,refimages,bgmask); + % + % Required inputs: + % absimages - Absorption image data, + % typically 16 bit grayscale images + % refimages - Raw reference image data + % absimages and refimages are both cell arrays containing + % 2D array data. The number of refimages can differ from the + % number of absimages. + % + % Optional inputs: + % bgmask - Array specifying background region used, + % 1=background, 0=data. Defaults to all ones. + % Outputs: + % optrefimages - Cell array of optimal reference images, + % equal in size to absimages. + % + + % Dependencies: none + % + % Authors: Shannon Whitlock, Caspar Ockeloen + % Reference: C. F. Ockeloen, A. F. Tauschinsky, R. J. C. Spreeuw, and + % S. Whitlock, Improved detection of small atom numbers through + % image processing, arXiv:1007.2136 + % Email: + % May 2009; Last revision: 11 August 2010 + + % Process inputs + + % Set variables, and flatten absorption and reference images + nimgs = size(absimages,3); + nimgsR = size(refimages,3); + xdim = size(absimages(:,:,1),2); + ydim = size(absimages(:,:,1),1); + + R = single(reshape(refimages,xdim*ydim,nimgsR)); + A = single(reshape(absimages,xdim*ydim,nimgs)); + optrefimages=zeros(size(absimages)); % preallocate + + if not(exist('bgmask','var')); bgmask=ones(ydim,xdim); end + k = find(bgmask(:)==1); % Index k specifying background region + + % Ensure there are no duplicate reference images + % R=unique(R','rows')'; % comment this line if you run out of memory + + % Decompose B = R*R' using singular value or LU decomposition + [L,U,p] = lu(R(k,:)'*R(k,:),'vector'); % LU decomposition + + for j=1:nimgs + b=R(k,:)'*A(k,j); + % Obtain coefficients c which minimise least-square residuals + lower.LT = true; upper.UT = true; + c = linsolve(U,linsolve(L,b(p,:),lower),upper); + + % Compute optimised reference image + optrefimages(:,:,j)=reshape(R*c,[ydim xdim]); + end +end \ No newline at end of file diff --git a/Estimations/OpticalAccordionLattice.m b/Estimations/OpticalAccordionLattice.m new file mode 100644 index 0000000..e432ba2 --- /dev/null +++ b/Estimations/OpticalAccordionLattice.m @@ -0,0 +1,141 @@ +%% 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');