diff --git a/EstimatesForAccordionLattice.m b/EstimatesForAccordionLattice.m new file mode 100644 index 0000000..41bb0ad --- /dev/null +++ b/EstimatesForAccordionLattice.m @@ -0,0 +1,114 @@ +%% Physical constants +PlanckConstant = 6.62607015E-34; +PlanckConstantReduced = 6.62607015E-34/(2*pi); +FineStructureConstant = 7.2973525698E-3; +ElectronMass = 9.10938291E-31; +GravitationalConstant = 6.67384E-11; +ProtonMass = 1.672621777E-27; +AtomicMassUnit = 1.660539066E-27; +BohrRadius = 5.2917721067E-11; +BohrMagneton = 9.274009994E-24; +BoltzmannConstant = 1.38064852E-23; +StandardGravityAcceleration = 9.80665; +SpeedOfLight = 299792458; +StefanBoltzmannConstant = 5.670373E-8; +ElectronCharge = 1.602176634E-19; +VacuumPermeability = 1.25663706212E-6; +DielectricConstant = 8.8541878128E-12; +ElectronGyromagneticFactor = -2.00231930436153; +AvogadroConstant = 6.02214076E23; +ZeroKelvin = 273.15; +GravitationalAcceleration = 9.80553; +VacuumPermittivity = 1 / (SpeedOfLight^2 * VacuumPermeability); +HartreeEnergy = ElectronCharge^2 / (4 * pi * VacuumPermittivity * BohrRadius); +AtomicUnitOfPolarizability = (ElectronCharge^2 * BohrRadius^2) / HartreeEnergy; % Or simply 4*pi*VacuumPermittivity*BohrRadius^3 + +% Dy specific constants +Dy164Mass = 163.929174751*1.660539066E-27; +Dy164IsotopicAbundance = 0.2826; +DyMagneticMoment = 9.93*9.274009994E-24; + +%% Lattice spacing + +Wavelength = 532e-9; +theta = linspace(1.5, 18.0, 100); +LatticeSpacing = Wavelength ./ (2.*sin((theta*pi/180)/2)); + +figure(1); +set(gcf,'Position',[100 100 950 750]) +plot(theta, LatticeSpacing * 1E6, LineWidth=2.0) +xlim([0 19]); +ylim([0.5 21]); +xlabel('Angle (deg)', FontSize=16) +ylabel('Lattice spacing (µm)', FontSize=16) +title(['\bf Upper bound = ' num2str(round(max(LatticeSpacing * 1E6),1)) ' µm ; \bf Lower bound = ' num2str(round(min(LatticeSpacing * 1E6),1)) ' µm'], FontSize=16) +grid on +%% Scaling of vertical trap frequency with lattice spacing +Wavelength = 532e-9; +a = 180 * (AtomicUnitOfPolarizability / (2 * SpeedOfLight * VacuumPermittivity)); +Power = 5; +waist_y = 250E-6; +waist_z = 50E-6; +TrapDepth = ((8 * a * Power) / (pi * waist_y * waist_z)) / (BoltzmannConstant * 1E-6); % in µK +thetas = linspace(1.5, 18.0, 100); +LatticeSpacings = zeros(1, length(thetas)); +Omega_z = zeros(1, length(thetas)); + +for idx = 1:length(thetas) + theta = 0.5 * thetas(idx) .* pi/180; + LatticeSpacings(idx) = Wavelength ./ (2.*sin(theta)); + Omega_z(idx) = sqrt(((16 * a * Power) / (pi * Dy164Mass * waist_y * waist_z)) * ... + ((2 * (cos(theta)/waist_z)^2) + ((Wavelength * sin(theta)/pi)^2 * ... + ((1/waist_y^4) + (1/waist_z^4))) + (pi / LatticeSpacings(idx))^2)); + +end + +nu_z = Omega_z ./ (2*pi); + +figure(2); +set(gcf,'Position',[100 100 950 750]) +plot(LatticeSpacings * 1E6, nu_z * 1E-3, LineWidth=2.0) +xlim([0.5 21]); +xlabel('Lattice spacing (µm)', FontSize=16) +ylabel('Trap frequency (kHz)', FontSize=16) +title(['\bf Upper bound = ' num2str(round(max(nu_z * 1E-3),2)) ' kHz ; \bf Lower bound = ' num2str(round(min(nu_z * 1E-3),2)) ' kHz'], FontSize=16) +grid on + +%% Scaling of Recoil Energy - All energy scales in an optical lattice are naturally parametrized by the lattice recoil energy + +LatticeSpacing = linspace(2E-6, 20E-6, 100); +RecoilEnergy = PlanckConstant^2 ./ (8 .* Dy164Mass .* LatticeSpacing.^2); + +figure(3); +set(gcf,'Position',[100 100 950 750]) +semilogy(LatticeSpacing * 1E6, RecoilEnergy/PlanckConstant, LineWidth=2.0, DisplayName=['\bf Max = ' num2str(round(max(RecoilEnergy / PlanckConstant),1)) ' Hz; Min = ' num2str(round(min(RecoilEnergy / PlanckConstant),1)) ' Hz']) +xlim([0.5 21]); +xlabel('Lattice spacing (µm)', FontSize=16) +ylabel('Recoil Energy (Hz)', FontSize=16) +title('\bf Scaling of Recoil Energy - All energy scales in an optical lattice are naturally parametrized by the lattice recoil energy', FontSize=12) +grid on +legend(FontSize=12) + +%% Interference pattern spacing in ToF - de Broglie wavelength associated with the relative motion of atoms + +ExpansionTime = linspace(1E-3, 20.0E-3, 100); + +figure(4); +set(gcf,'Position',[100 100 950 750]) +labels = []; + +for ls = [2E-6:2E-6:5E-6 6E-6:6E-6:20E-6] + InteferencePatternSpacing = (PlanckConstant .* ExpansionTime) ./ (Dy164Mass * ls); + plot(ExpansionTime*1E3, InteferencePatternSpacing* 1E6, LineWidth=2.0, DisplayName=['\bf Lattice spacing = ' num2str(round(max(ls * 1E6),1)) ' µm']) + hold on +end +xlim([0 22]); +xlabel('Free expansion time (milliseconds)', FontSize=16) +ylabel('Interference pattern period (µm)', FontSize=16) +title('\bf Interference of condensates - Fringe period is the de Broglie wavelength associated with the relative motion of atoms', FontSize=12) +legend(labels, 'Location','NorthWest', FontSize=12); +grid on +legend show + + + diff --git a/plotImages.m b/plotImages.m new file mode 100644 index 0000000..f835ff9 --- /dev/null +++ b/plotImages.m @@ -0,0 +1,223 @@ +%% Parameters + +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"]; + +folderPath = "C:/Users/Karthik/Documents/GitRepositories/Calculations/24/"; + +run = '0086'; + +folderPath = strcat(folderPath, run); + +cam = 5; + +angle = 90; +center = [2100, 1150]; +span = [500, 500]; +fraction = [0.1, 0.1]; + +pixel_size = 4.6e-6; + +%% Compute OD image, rotate and extract ROI for analysis +% Get a list of all files in the folder with the desired file name pattern. +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 = im2double(imrotate(h5read(fullFileName, append(groupList(cam), "/atoms")), angle)); + bkg_img = im2double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle)); + dark_img = im2double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle)); + + refimages(:,:,k) = subtract_offset(crop_image(bkg_img, center, span), fraction); + absimages(:,:,k) = subtract_offset(crop_image(calculate_OD(atm_img, bkg_img, dark_img), center, span), fraction); + +end +%% Fringe removal + +optrefimages = removefringesInImage(absimages, refimages); +absimages_fringe_removed = absimages(:, :, :) - optrefimages(:, :, :); + +nimgs = size(absimages_fringe_removed,3); +od_imgs = cell(1, nimgs); +for i = 1:nimgs + od_imgs{i} = absimages_fringe_removed(:, :, i); +end + +%% +figure(1) +clf +r = 120; +x = 250; +y = 250; +for k = 1 : length(od_imgs) + imagesc(xvals, yvals, od_imgs{k}) + hold on + th = 0:pi/50:2*pi; + xunit = r * cos(th) + x; + yunit = r * sin(th) + y; + h = plot(xunit, yunit, Color='yellow'); + xlabel('µm', 'FontSize', 16) + ylabel('µm', 'FontSize', 16) + axis equal tight; + hcb = colorbar; + hL = ylabel(hcb, 'Optical Density', 'FontSize', 16); + set(hL,'Rotation',-90); + colormap jet; + set(gca,'CLim',[0 1.0]); + set(gca,'YDir','normal') + title('DMD projection: Circle of radius 200 pixels', 'FontSize', 16); + + drawnow; +end + +%% Helper Functions + +function ret = get_offset_from_corner(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 = subtract_offset(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 = get_offset_from_corner(img, x_fraction, y_fraction); + ret = img - offset; +end + +function ret = crop_image(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 ret = calculate_OD(imageAtom, imageBackground, imageDark) + % Calculate the OD image for absorption imaging. + % :param imageAtom: The image with atoms + % :type imageAtom: numpy array + % :param imageBackground: The image without atoms + % :type imageBackground: numpy array + % :param imageDark: The image without light + % :type imageDark: numpy array + % :return: The OD images + % :rtype: numpy array + + numerator = imageBackground - imageDark; + denominator = imageAtom - imageDark; + + numerator(numerator == 0) = 1; + denominator(denominator == 0) = 1; + + ret = -log(double(abs(denominator ./ numerator))); + + if numel(ret) == 1 + ret = ret(1); + 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