clear all; close all; clc;

iFluid_path = ['..' filesep '..' filesep '..'];

addpath([iFluid_path filesep 'models' filesep])
addpath([iFluid_path filesep 'solvers' filesep])
addpath([iFluid_path filesep 'iFluid' filesep])
addpath([iFluid_path filesep 'utils' filesep])


set(0, 'DefaultTextInterpreter', 'latex');
set(0, 'DefaultLegendInterpreter', 'latex');
set(0, 'defaultAxesTickLabelInterpreter','latex');
set(0, 'defaultAxesFontSize',12);


%% Analyse experimental data files

data_dir    = 'Z:\Philipp\Drude_weight\step\exp_data';

filelist    = dir( [ data_dir filesep '*.mat'] );
Nfiles      = length(filelist); 


m_si        = 87*1.6605402e-27;
hbar_si     = 1.054571726e-34;
kB_si       = 1.38065e-23;
as_si       = 5.2e-9;


deltaN      = cell(1, Nfiles);
t_ax        = cell(1, Nfiles);

dn_fit      = zeros(1, Nfiles);
dmu_fit     = zeros(1, Nfiles);
dmu_err     = zeros(1, Nfiles);

for i = 1:Nfiles

    % load experimental data (stored in a struct named "clean_data")
    load([filelist(i).folder filesep filelist(i).name]) 

    % rename experimental data fields for cleaner code
    omegaT_si   = 2*pi*clean_data.transv_trap_freq_ver_Hz;
  
    z_axis      = clean_data.z_axis_si;
    t_axis      = clean_data.t_axis_si;
    n_shots     = clean_data.density_profiles_TOF_si_z_shots_t;
    n_mean      = squeeze(mean(n_shots,2,'omitnan'));

    if all(clean_data.scan_name(5:9) == '12054')
        disp('hej')
        n_mean = flip(n_mean, 1);

        t_axis = t_axis';
    end


    % fit position of box walls
    [z_wall, w_wall] = fit_box_walls(z_axis', mean(n_mean, 2));
    z_left      = z_wall(1) + 1*w_wall(1);
    z_right     = z_wall(2) - 1*w_wall(2); 

    % fit step to get mean density and midpoint
    z_fit_int   = 0.65*z_wall;
    [~, fitres] = fit_step_box(z_axis, n_mean(:,1), z_fit_int);
    fitcoeffs   = coeffvalues(fitres);
    ci          = confint(fitres, 0.95);

    n0          = fitcoeffs(4)*1e-6;
    dn          = fitcoeffs(1)/2 *1e-6;
    z0          = fitcoeffs(3)*1e6;
    
    % crop walls out of data and smoothen
    crop        = z_axis > z_left & z_axis < z_right;
    z           = z_axis(crop)*1e6 - z0;
    t           = t_axis*1e3;
    n           = n_mean(crop, :)*1e-6;
    % n           = smoothdata(n, 1, "movmean", 4);

    % calculate coupling
    lperp_si    = sqrt(hbar_si/m_si/omegaT_si);
    trans_broad = 0.5*(2 + 3*as_si*n0*1e6)/(1 + 2*as_si*n0*1e6)^(3/2);
    c_si        = 2*as_si/lperp_si^2 * trans_broad;
    g_si        = hbar_si^2*c_si/m_si;


    N_L         = trapz(z(z < 0), n(z < 0, :), 1);
    N_R         = trapz(z(z > 0), n(z > 0, :), 1);

    deltaN{i}   = N_L - N_R;
    t_ax{i}     = t;
    dn_fit(i)   = dn;
    dmu_fit(i)  = -dn*g_si*1e6 / (kB_si*1e-9);
    dmu_err(i)  = (ci(2,1) - fitcoeffs(1))/2*g_si/ (kB_si*1e-9);
end


%%



fig = figure;
fig.Units = 'centimeters';
fig.Position = [ 20 10 8.6 6 ];
set(fig,'renderer','Painters');

hold on
box on
for i = 1:4 %Nfiles

    test = deltaN{i}(1);
    if i == Nfiles
        test = 0.35*test;
    end

    plot(t_ax{i}, -(deltaN{i}-deltaN{i}(1))/test, ...
        'LineStyle', ':', ...
        'LineWidth', 1.25, ...
        'Marker', 'o', ...
        'MarkerSize', 3, ...
        'MarkerFaceColor',[1,1,1])
end

xlim([0, 20])
grid on

plot(t, 0.045*t, 'k--')


set(gca, 'XScale', 'log')
set(gca, 'YScale', 'log')
% xlim([1e0, 1e2])
ylim([1e-2, 1e0])

xlabel('$t \; [\mathrm{ms}]$')
ylabel('$\frac{\Delta N(t) - \Delta N(0)}{\Delta N(0)}$')


legendEntries = {};
for i = 1:4 %length(dmu_fit)
    legendEntries{i} = ['$\delta\mu = ' num2str(dmu_fit(i)) '\: [ \mathrm{nK} ]$' ];
end



legend(legendEntries, 'location', 'southeast', 'Interpreter','latex')