function [dVdz_fit, dVdz_std, cfit, gof] = fit_potential_slope(z_grid, dens_slope_si, omegaT_si, z_int)
% Obtain slope of linear potential by fitting density profile.
%
% Input arguments:
%   z_grid - position grid
%           should be in SI-units
%   dens_slope_si - equilibrium density profiles 
%           should be in SI-units
%           should have indices (position, shots)
%   g_si - 1D coupling strength
%           should be in SI-units
%   z_int - (optional) interval of spatial points to use for fit
%           should be in SI-units
%
% Output arguments:
%   dVdz_fit - slope of potential obtained from fit
%   dVdz_std - "standard deviation" (68% CI) of fit
%   cfit - curve-fitting object
%   gof - goodness-of-fit data


if nargin < 4
    z_int = [z_grid(1), z_grid(end)];
end

% slope_dens_si should have indices (position, shots) 
dens_slope_mean = mean(dens_slope_si, 2, "omitnan"); 
dens_slope_std  = std(dens_slope_si, 0, 2, "omitnan");

% setup fitoptions
% fo = fitoptions('Method', 'NonlinearLeastSquares',...
%            'Lower', [0,0],...
%            'Upper', [Inf,max(cdate)],...
%            'StartPoint', [1 1]);

fo = fitoptions('Method', 'NonlinearLeastSquares',...
                'Exclude', z_grid < z_int(1) | z_grid > z_int(end) );

% linear fit of the density profile 
[cfit, gof]  = fit(z_grid(:), dens_slope_mean(:), 'poly1', fo)

% get 95% prediction interval for plotting
pred_int = predint(cfit, z_grid, 0.95);

% get slope and confidence interval of slope
vals    = coeffvalues(cfit);
confs   = confint(cfit, 0.95);
stds    = (confs(2,:)-confs(1,:))/(2 * 1.96); % convert from CI to std

% plot the original data and the linear fit
figure
hold on
box on
plot(cfit, z_grid, dens_slope_mean)
plot(z_grid, pred_int, 'm--')
xline(z_int(1), 'k:')
xline(z_int(2), 'k:')
legend('Data','Linear Fit','95% Prediction Interval','')
xlabel('position')
ylabel('density')

% calculate coupling based on the mean density
m_si        = 87*1.6605402e-27;
hbar_si     = 1.054571726e-34;
as_si       = 5.2e-9;

n0_si       = vals(2);

trans_broad = 0.5*(2 + 3*as_si*n0_si)/(1 + 2*as_si*n0_si)^(3/2);

g_si        = 2*as_si*hbar_si*omegaT_si;
g_eff       = g_si * trans_broad;


dg_dn       = g_si*(-0.375*as_si*(1 + as_si*n0_si))/((0.5 + as_si*n0_si)^2 * sqrt(1 + 2*as_si*n0_si));

% use \mu = g*n to get potential slope from density slope;

dVdz_fit    = vals(1)*g_eff; % in SI units: [energy]/[length]
dVdz_std    = sqrt( stds(1)^2 * g_eff^2 + dg_dn^2 *stds(2)^2 ); 

 

end