classdef DensityProfileBEC2DModel < handle %% DensityProfileBEC2DModel % Author: Jianshun Gao % Date: 2025-09-12 % Version: 1.0 % % Description: % This class provides methods to model, fit, and analyze 2D Bose-Einstein % condensate (BEC) density profiles, including both condensed (Thomas-Fermi) % and thermal (polylog) components. % % Properties: % - params: Structure storing parameter values and bounds. % - fitResult: MATLAB fit object after fitting. % - gof: Goodness-of-fit structure from the fit. % % Methods: % - cal_cond_frac(X, Y): Compute condensate fraction from fitted profile. % - return_atom_number(X, Y): Compute total, BEC, and thermal atom numbers. % - return_temperature(tof, omega): Estimate temperature from thermal cloud width. % % Static Methods: % - ThomasFermi_2d: 2D Thomas-Fermi parabolic profile. % - polylog: Series approximation of polylogarithm function. % - polylog2_2d: 2D thermal cloud profile using polylog(2) function. % - density_profile_BEC_2d: Full 2D density (BEC + thermal), optionally rotated. % - density_1d: 1D slice of the density profile combining BEC and thermal parts. % % Usage: % obj = DensityProfileBEC2DModel(); % atom_n = obj.return_atom_number(X, Y); % Atom numbers % cond_frac = obj.cal_cond_frac(X, Y); % Condensate fraction % T = obj.return_temperature(tof, omega); % Temperature % % Notes: % - All static methods can be called independently without creating an object. properties % Conversion factors and default constants fwhm_factor = 2*sqrt(2*log(2)); % FWHM factor for Gaussian height_factor = 1/(2*pi); % Height normalization factor prefix = ''; % Optional parameter prefix tiny = 1e-15; % Small number to avoid division by zero params; % Struct for parameter definitions fitResult; % Result of 2D fit cond_frac; % Condensate fraction gof; % Goodness-of-fit structure atom_n_conv = 144; % Conversion factor for atom number pre_check = false; % Flag for pre-fit check post_check = false; % Flag for post-fit check is_debug = false; % Flag for debug plots end methods function obj = DensityProfileBEC2DModel(prefix) % Constructor if nargin > 0 obj.prefix = prefix; end obj.set_paramhints_prefix(); end function set_paramhints_prefix(obj) % Initialize parameter hints and expressions obj.params = struct(); % Minimum values obj.params.([obj.prefix 'amp_bec']).min = 0; obj.params.([obj.prefix 'amp_th']).min = 0; obj.params.([obj.prefix 'x0_bec']).min = 0; obj.params.([obj.prefix 'y0_bec']).min = 0; obj.params.([obj.prefix 'x0_th']).min = 0; obj.params.([obj.prefix 'y0_th']).min = 0; obj.params.([obj.prefix 'sigmax_bec']).min = 0; obj.params.([obj.prefix 'sigmay_bec']).min = 0; obj.params.([obj.prefix 'sigma_th']).min = 0; obj.params.([obj.prefix 'rot_angle']).min = -90; obj.params.([obj.prefix 'rot_angle']).max = 90; % Expressions for derived parameters obj.params.([obj.prefix 'atom_number_bec']).expr = ... [obj.prefix 'amp_bec / 5 * 2 * pi * ' obj.prefix 'sigmax_bec * ' obj.prefix 'sigmay_bec']; obj.params.([obj.prefix 'atom_number_th']).expr = ... [obj.prefix 'amp_th * 2 * pi * 1.20206 / 1.643 * ' obj.prefix 'sigma_th * ' obj.prefix 'sigma_th']; obj.params.([obj.prefix 'condensate_fraction']).expr = ... [obj.prefix 'atom_number_bec / (' obj.prefix 'atom_number_bec + ' obj.prefix 'atom_number_th)']; end function params = guess(obj, data, x, y, varargin) % Estimate initial parameters for BEC + thermal cloud fit. % Performs 1D slicing along the shorter axis of the image and initial amplitude guesses. p = inputParser; addParameter(p, 'negative', false); addParameter(p, 'pureBECThreshold', 0.5); addParameter(p, 'noBECThThreshold', 0.0); addParameter(p, 'rot_angle', 0); addParameter(p, 'vary_rot', false); addParameter(p, 'is_debug', false); addParameter(p, 'pre_check', false); addParameter(p, 'post_check', false); parse(p, varargin{:}); obj.is_debug = p.Results.is_debug; obj.pre_check = p.Results.pre_check; obj.post_check = p.Results.post_check; x_width = length(unique(x)); y_width = length(unique(y)); x_1d = linspace(x(1), x(end), x_width); y_1d = linspace(y(1), y(end), y_width); data_2d = reshape(data, y_width, x_width)'; % Rotate image if needed rot_angle = p.Results.rot_angle; if rot_angle ~= 0 data_2d = imrotate(data_2d, rot_angle, 'bilinear', 'crop'); end % Debug plot of input data if obj.is_debug figure; imagesc(x_1d, y_1d, data_2d'); axis equal tight; colorbar; colormap('jet'); title('Input Data'); xlabel('x-axis'); ylabel('y-axis'); end % Binarize and locate BEC thresh = obj.calc_thresh(data_2d, 0.5); center_pix = obj.calc_cen_pix(thresh); center = obj.center_pix_conv(center_pix, x_1d, y_1d); BEC_width_guess = obj.guess_BEC_width(thresh, center_pix); if obj.is_debug figure; imagesc(x_1d, y_1d, thresh'); hold on; plot(center(1), center(2), 'gx', 'MarkerSize', 25, 'LineWidth', 2); axis equal tight; colorbar; colormap('jet'); title(sprintf('Binarized Image (BEC Width: x=%.0f, y=%.0f pixels)', BEC_width_guess(1), BEC_width_guess(2))); xlabel('x-axis'); ylabel('y-axis'); end % 1D slicing along the shorter axis if BEC_width_guess(1) < BEC_width_guess(2) s_width_ind = 1; x_fit = x_1d; slice_range = round(center_pix(2) - BEC_width_guess(2)/2) : round(center_pix(2) + BEC_width_guess(2)/2); X_guess = sum(data_2d(:, slice_range), 2) / length(slice_range); else s_width_ind = 2; x_fit = y_1d; slice_range = round(center_pix(1) - BEC_width_guess(1)/2) : round(center_pix(1) + BEC_width_guess(1)/2); X_guess = sum(data_2d(slice_range, :), 1) / length(slice_range); end max_val = max(X_guess); % Construct initial parameter struct params_1d = struct(); params_1d.x0_bec.value = center(s_width_ind); params_1d.x0_bec.min = center(s_width_ind) - 10; params_1d.x0_bec.max = center(s_width_ind) + 10; params_1d.x0_th.value = center(s_width_ind); params_1d.x0_th.min = center(s_width_ind) - 10; params_1d.x0_th.max = center(s_width_ind) + 10; params_1d.amp_bec.value = 0.5 * max_val; params_1d.amp_bec.min = 0; params_1d.amp_bec.max = 1.3 * max_val; params_1d.amp_th.value = 0.5 * max_val; params_1d.amp_th.min = 0; params_1d.amp_th.max = 1.3 * max_val; params_1d.sigma_bec.value = BEC_width_guess(s_width_ind) / 1.22; params_1d.sigma_bec.min = 0; params_1d.sigma_bec.max = 2 * BEC_width_guess(s_width_ind); params_1d.sigma_th.value = 0.632 * params_1d.sigma_bec.value + 0.518 * 3 * BEC_width_guess(s_width_ind); params_1d.sigma_th.min = 0; params_1d.sigma_th.max = 3 * (0.632 * params_1d.sigma_bec.value + 0.518 * 3 * BEC_width_guess(s_width_ind)); % Perform 1D bimodal fit [fitResult_1d, gof_1d] = fit_1d_bimodal(x_fit, X_guess, params_1d); % Extract fit coefficients bval_1d = coeffvalues(fitResult_1d); paramNames_1d = coeffnames(fitResult_1d); for i = 1:length(paramNames_1d) bval_1d_struct.(paramNames_1d{i}) = bval_1d(i); end if obj.is_debug disp('Initial 1D fit parameters:'); disp(bval_1d); figure; plot(x_fit, X_guess, 'b-', 'LineWidth', 2); hold on; plot(x_fit, obj.density_1d(x_fit, bval_1d.x0_bec, bval_1d.x0_th, ... bval_1d.amp_bec, bval_1d.amp_th, bval_1d.sigma_bec, bval_1d.sigma_th), 'r-', 'LineWidth', 2); legend('Data', 'Fit'); if s_width_ind == 1 xlabel('x-axis'); title('1D Fit along x-axis'); else xlabel('y-axis'); title('1D Fit along y-axis'); end ylabel('Intensity'); end % Final parameter preparation (pre_check logic omitted for brevity) params = params_1d; obj.params = params; end function [fitResult, gof] = fit(obj, data, x, y, varargin) % Perform 2D nonlinear least squares fit for BEC + thermal cloud. % Optionally includes rotation of the cloud profile. data = double(data); figure imagesc(data); axis equal tight; colorbar; colormap('jet'); if isempty(obj.params) obj.guess(data, x, y, varargin{:}); end % Prepare grid and data [X, Y] = meshgrid(x, y); xyData = [X(:), Y(:)]; zData = double(data(:)); % Create fit options and fittype options = fitoptions('Method', 'NonlinearLeastSquares'); if obj.params.([obj.prefix 'rot_angle']).vary % Include rotation parameter paramOrder = {[obj.prefix 'amp_bec'], [obj.prefix 'amp_th'], ... [obj.prefix 'x0_bec'], [obj.prefix 'y0_bec'], ... [obj.prefix 'x0_th'], [obj.prefix 'y0_th'], ... [obj.prefix 'sigmax_bec'], [obj.prefix 'sigmay_bec'], ... [obj.prefix 'sigma_th'], [obj.prefix 'rot_angle']}; % Start point, lower and upper bounds startPoint = zeros(1, length(paramOrder)); lowerBounds = zeros(1, length(paramOrder)); upperBounds = zeros(1, length(paramOrder)); for i = 1:length(paramOrder) paramName = paramOrder{i}; startPoint(i) = obj.params.(paramName).value; lowerBounds(i) = obj.params.(paramName).min; upperBounds(i) = obj.params.(paramName).max; end options.StartPoint = startPoint; options.Lower = lowerBounds; options.Upper = upperBounds; ft = fittype(@(amp_bec, amp_th, x0_bec, y0_bec, x0_th, y0_th, ... sigmax_bec, sigmay_bec, sigma_th, rot_angle, x, y) ... obj.density_profile_BEC_2d(x, y, amp_bec, amp_th, x0_bec, y0_bec, ... x0_th, y0_th, sigmax_bec, sigmay_bec, sigma_th, rot_angle), ... 'independent', {'x', 'y'}, 'dependent', 'z'); else % No rotation paramOrder = {[obj.prefix 'amp_bec'], [obj.prefix 'amp_th'], ... [obj.prefix 'x0_bec'], [obj.prefix 'y0_bec'], ... [obj.prefix 'x0_th'], [obj.prefix 'y0_th'], ... [obj.prefix 'sigmax_bec'], [obj.prefix 'sigmay_bec'], ... [obj.prefix 'sigma_th']}; startPoint = zeros(1, length(paramOrder)); lowerBounds = zeros(1, length(paramOrder)); upperBounds = zeros(1, length(paramOrder)); for i = 1:length(paramOrder) paramName = paramOrder{i}; startPoint(i) = obj.params.(paramName).value; lowerBounds(i) = obj.params.(paramName).min; upperBounds(i) = obj.params.(paramName).max; end options.StartPoint = startPoint; options.Lower = lowerBounds; options.Upper = upperBounds; ft = fittype(@(amp_bec, amp_th, x0_bec, y0_bec, x0_th, y0_th, ... sigmax_bec, sigmay_bec, sigma_th, x, y) ... obj.density_profile_BEC_2d(x, y, amp_bec, amp_th, x0_bec, y0_bec, ... x0_th, y0_th, sigmax_bec, sigmay_bec, sigma_th, 0), ... 'independent', {'x', 'y'}, 'dependent', 'z'); end [obj.fitResult, obj.gof] = fit(xyData, zData, ft, options); fitResult = obj.fitResult; gof = obj.gof; % Compute condensate fraction obj.cond_frac = obj.cal_cond_frac(X, Y); end function thresh = calc_thresh(obj, data, thresh_val, sigma) % Binarize image for BEC detection if nargin < 3, thresh_val = 0.3; end if nargin < 4, sigma = 0.4; end blurred = imgaussfilt(data, sigma); thresh = blurred < max(blurred(:)) * thresh_val; thresh = double(~thresh); end function center_pix = calc_cen_pix(obj, thresh) % Compute center of binarized image [X, Y] = size(thresh); thresh = thresh / sum(thresh(:)); dx = sum(thresh, 2); dy = sum(thresh, 1); center_pix = [sum(dx .* (1:X)'), sum(dy .* (1:Y))]; end function center = center_pix_conv(obj, center_pix, x, y) % Convert pixel indices to coordinate values center = [x(round(center_pix(1))), y(round(center_pix(2)))]; end function BEC_width_guess = guess_BEC_width(obj, thresh, center) % Estimate BEC width along x and y [X, Y] = size(thresh); BEC_width_guess = [sum(thresh(:, round(center(2)))), sum(thresh(round(center(1)), :))]; BEC_width_guess(BEC_width_guess<=0) = 1; end function cond_frac = cal_cond_frac(obj, X, Y) % Compute condensate fraction from fitted BEC + thermal cloud profile bval = coeffvalues(obj.fitResult); paramNames = coeffnames(obj.fitResult); % Extract fit parameters for i = 1:length(paramNames) eval([paramNames{i} ' = bval(i);']); end if ~obj.params.([obj.prefix 'rot_angle']).vary rot_angle = 0; end % Thomas-Fermi fit for condensate tf_fit = obj.ThomasFermi_2d(X, Y, x0_bec, y0_bec, amp_bec, sigmax_bec, sigmay_bec); % Total density profile (BEC + thermal cloud) fit_total = obj.density_profile_BEC_2d(X, Y, amp_bec, amp_th, x0_bec, y0_bec, ... x0_th, y0_th, sigmax_bec, sigmay_bec, sigma_th, rot_angle); N_bec = sum(tf_fit(:)); N_ges = sum(fit_total(:)); cond_frac = N_bec / N_ges; end function atom_n = return_atom_number(obj, X, Y, is_print) % Compute total number of atoms, BEC atoms, thermal atoms, and condensate fraction if nargin < 4 is_print = true; end bval = coeffvalues(obj.fitResult); paramNames = coeffnames(obj.fitResult); % Extract fit parameters for i = 1:length(paramNames) eval([paramNames{i} ' = bval(i);']); end tf_fit = obj.ThomasFermi_2d(X, Y, x0_bec, y0_bec, amp_bec, sigmax_bec, sigmay_bec); th_fit = obj.polylog2_2d(X, Y, x0_th, y0_th, amp_th, sigma_th, sigma_th); N_bec = obj.atom_n_conv * sum(tf_fit(:)); N_th = obj.atom_n_conv * sum(th_fit(:)); N = N_bec + N_th; frac = N_bec / N; if is_print fprintf('\nAtom numbers:\n'); fprintf(' N_bec: %.0f\n', N_bec); fprintf(' N_th: %.0f\n', N_th); fprintf(' N_total: %.0f\n', N); fprintf(' Condensate fraction: %.2f %%\n', frac * 100); end atom_n = struct('N', N, 'N_bec', N_bec, 'N_th', N_th, 'cond_f', frac); end function T = return_temperature(obj, tof, omg, is_print, eff_pix) % Compute temperature from thermal cloud width and time-of-flight if nargin < 3 omg = []; end if nargin < 4 is_print = true; end if nargin < 5 eff_pix = 2.472e-6; end bval = coeffvalues(obj.fitResult); paramNames = coeffnames(obj.fitResult); % Extract fit parameters for i = 1:length(paramNames) eval([paramNames{i} ' = bval(i);']); end R_th = sigma_th * eff_pix * sqrt(2); % Thermal cloud radius % Physical constants u = 1.66053906660e-27; % Atomic mass unit [kg] k = 1.380649e-23; % Boltzmann constant [J/K] if isempty(omg) % Free expansion T = R_th^2 * 164 * u / k / tof^2; else % Trap expansion included T = R_th^2 * 164 * u / k / (1/omg^2 + tof^2); end if is_print fprintf('Temperature: %.2f nK\n', T * 1e9); end end end methods (Static) function res = ThomasFermi_2d(x, y, centerx, centery, amplitude, sigmax, sigmay) % 2D Thomas-Fermi distribution for a BEC (parabolic density profile) res = (1 - ((x - centerx)/sigmax).^2 - ((y - centery)/sigmay).^2); res(res < 0) = 0; res = amplitude * res.^(3/2); % Standard TF 3/2 exponent end function res = polylog(power, numerator) % Approximate polylogarithm function using truncated series order = 20; % Truncation order dataShape = size(numerator); numerator = repmat(numerator(:), 1, order); numerator = numerator .^ repmat(1:order, prod(dataShape), 1); denominator = repmat((1:order), prod(dataShape), 1); data = numerator ./ (denominator .^ power); res = sum(data, 2); res = reshape(res, dataShape); end function res = polylog2_2d(x, y, centerx, centery, amplitude, sigmax, sigmay) % 2D thermal cloud distribution using polylog(2) arg = exp(-((x - centerx).^2 / (2*sigmax^2)) - ((y - centery).^2 / (2*sigmay^2))); res = amplitude / 1.643 * FitModels.DensityProfileBEC2DModel.polylog(2, arg); end function res = density_profile_BEC_2d(x, y, amp_bec, amp_th, x0_bec, y0_bec, ... x0_th, y0_th, sigmax_bec, sigmay_bec, sigma_th, rot_angle) % Combined 2D density profile: BEC (Thomas-Fermi) + thermal cloud (polylog) if nargin < 12 rot_angle = 0; end % Rotate coordinates if needed if rot_angle ~= 0 rot_angle_rad = -rot_angle * pi/180; % Clockwise rotation % Rotate coordinates x_rot = x * cos(rot_angle_rad) + y * sin(rot_angle_rad); y_rot = -x * sin(rot_angle_rad) + y * cos(rot_angle_rad); % Rotate centers x0_bec_rot = x0_bec * cos(rot_angle_rad) + y0_bec * sin(rot_angle_rad); y0_bec_rot = -x0_bec * sin(rot_angle_rad) + y0_bec * cos(rot_angle_rad); x0_th_rot = x0_th * cos(rot_angle_rad) + y0_th * sin(rot_angle_rad); y0_th_rot = -x0_th * sin(rot_angle_rad) + y0_th * cos(rot_angle_rad); x = x_rot; y = y_rot; x0_bec = x0_bec_rot; y0_bec = y0_bec_rot; x0_th = x0_th_rot; y0_th = y0_th_rot; end % Thomas-Fermi part (BEC) TF_part = FitModels.DensityProfileBEC2DModel.ThomasFermi_2d(x, y, x0_bec, y0_bec, amp_bec, sigmax_bec, sigmay_bec); % Polylog thermal part poly_part = FitModels.DensityProfileBEC2DModel.polylog2_2d(x, y, x0_th, y0_th, amp_th, sigma_th, sigma_th); % Total density res = TF_part + poly_part; end function res = density_1d(x, x0_bec, x0_th, amp_bec, amp_th, sigma_bec, sigma_th) % 1D density profile (Thomas-Fermi + thermal polylog) thermal_part = amp_th / 1.643 * polylog_int(exp(-0.5 * (x - x0_th).^2 / sigma_th^2)); TF_part = amp_bec * (1 - ((x - x0_bec)/sigma_bec).^2); TF_part(TF_part < 0) = 0; TF_part = TF_part.^(3/2); res = thermal_part + TF_part; end end end function res = polylog_int(x) x_int = linspace(0, 1.00001, 1000); poly_tab = zeros(size(x_int)); for i = 1:length(x_int) poly_tab(i) = sum((x_int(i).^(1:20)) ./ (1:20).^2); end res = interp1(x_int, poly_tab, x, 'linear', 'extrap'); end function [fitResult, gof] = fit_1d_bimodal(x, y, initialParams) bimodal1d = @(amp_bec, amp_th, x0_bec, x0_th, sigma_bec, sigma_th, x) ... FitModels.DensityProfileBEC2DModel.density_1d(x, x0_bec, x0_th, amp_bec, amp_th, sigma_bec, sigma_th); paramNames = {'amp_bec', 'amp_th', 'x0_bec', 'x0_th', 'sigma_bec', 'sigma_th'}; ft = fittype(bimodal1d, 'independent', 'x', 'dependent', 'y'); options = fitoptions(ft); % paramNames = fieldnames(initialParams); startPoint = zeros(1, length(paramNames)); lowerBounds = zeros(1, length(paramNames)); upperBounds = zeros(1, length(paramNames)); for i = 1:length(paramNames) paramName = paramNames{i}; startPoint(i) = initialParams.(paramName).value; lowerBounds(i) = initialParams.(paramName).min; upperBounds(i) = initialParams.(paramName).max; end options.StartPoint = startPoint; options.Lower = lowerBounds; options.Upper = upperBounds; [fitResult, gof] = fit(x(:), y(:), ft, options); end