tic
% Preparation
format default
clear; clc
sub = readmatrix('substrate_trans.txt');
samp = readmatrix('sample_trans.txt');
%% Preamble
% Fundamental constants
h = 4.1356692*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
% Data
T_sample = samp(:, 2);
T_substrate = sub(:, 2);
lambda = sub(:, 1)*10^-9;
E = h*c./lambda;
Absorbance = log10(T_substrate./T_sample);
% Data for fitting
min_E = 1.5;
max_E = 2.0;
Range_E = E >= min_E & E <= max_E;
E_p = E(Range_E);
Abs0 = Absorbance(find(max(E_p) == E):find(min(E_p) == E));
Abs = Abs0 - min(Abs0);
% Fitting function: Elliott's Model for Steady State with non-parabolic factor
function F = EM_SS_wR(p, e_p)
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*...
(integral(@(E)sech(((E_p - E)./p(6)))*(1 + 10*p(5)*(E - p(3)) + ...
126*(p(5))^2*(E - p(3))^2)/(1 - exp(-2*pi*sqrt(p(4)/(E - p(3))))), p(3), Inf, 'ArrayValued', 1)) + ...
p(2)*((2*pi*p(4)^3/2)/E_p)*1/p(6)*(...
(1/1^3)*sech((E_p - p(3) + p(4)/1^2)./p(6)) + ...
(1/2^3)*sech((E_p - p(3) + p(4)/2^2)./p(6)) + ...
(1/3^3)*sech((E_p - p(3) + p(4)/3^2)./p(6)) + ...
(1/4^3)*sech((E_p - p(3) + p(4)/4^2)./p(6)) + ...
(1/5^3)*sech((E_p - p(3) + p(4)/5^2)./p(6)) + ...
(1/6^3)*sech((E_p - p(3) + p(4)/6^2)./p(6)) + ...
(1/7^3)*sech((E_p - p(3) + p(4)/7^2)./p(6)));
end
F = F(:);
end
% Carrier Contribution
function F = EM_SS_wR_CC(p, e_p)
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*...
(integral(@(E)sech(((E_p - E)./p(6)))*(1 + 10*p(5)*(E - p(3)) + ...
126*(p(5))^2*(E - p(3))^2)/(1 - exp(-2*pi*sqrt(p(4)/(E - p(3))))), p(3), Inf, 'ArrayValued', 1));
end
F = F(:);
end
% Excitoninc Contribution
function F = EM_SS_wR_EC(p, E_p)
F = p(2)*((2*pi*p(4)^3/2)/E_p)*1/p(6)*(...
(1/1^3)*sech((E_p - p(3) + p(4)/1^2)./p(6)) + ...
(1/2^3)*sech((E_p - p(3) + p(4)/2^2)./p(6)) + ...
(1/3^3)*sech((E_p - p(3) + p(4)/3^2)./p(6)) + ...
(1/4^3)*sech((E_p - p(3) + p(4)/4^2)./p(6)) + ...
(1/5^3)*sech((E_p - p(3) + p(4)/5^2)./p(6)) + ...
(1/6^3)*sech((E_p - p(3) + p(4)/6^2)./p(6)) + ...
(1/7^3)*sech((E_p - p(3) + p(4)/7^2)./p(6)));
end
%% Curve fitting options
% Initial parameter guess and bounds
lb = [-Inf, -Inf, 0, 0, -Inf, -Inf]; ub = [Inf, Inf, 3, 0.030, Inf, Inf];
p0 = [0.13, 0.11, 1.61, 0.025, 0.12, 0.035]; % refer to the next line for their order
% p0 = [A1, A2, Eg, Eb, R, g]
% lsqcurvefit and choose between different algorithm that lsqcurvefit employs (3C1, comment those lines that are not choosen and uncomment the line that is choosen, if not, matlab will take the last line of "optim_lsq" by default)
optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt', 'MaxFunctionEvaluations',1000, 'MaxIterations', 1000, 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', 'MaxFunctionEvaluations',1000, 'MaxIterations', 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'interior-point', 'MaxFunctionEvaluations',1000, 'MaxIterations', 1000, 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% Solver for lsqcurvefit
[p, residualnorm, residual, exitflag, output, lambda, jacobian] = lsqcurvefit(@(p, e_p) EM_SS_wR(p, e_p), p0, E_p, Abs, lb, ub);
%% Error bars calculation
ci = nlparci(p, residual, 'Jacobian', jacobian);
PCI = table(ci(:,1), p(:), ci(:,2),'VariableNames',{'CI 0.025','p','CI 0.975'});
Parameter_CI = table2array(PCI);
upper_bar = (Parameter_CI(:,3) - Parameter_CI(:,2))./2;
lower_bar = (Parameter_CI(:,2) - Parameter_CI(:,1))./2;
%% Parameter values (refer to command window)
p1 = p(1,1);
p2 = p(1,2);
p3 = p(1,3);
p4 = p(1,4);
p5 = p(1,5);
p6 = p(1,6);
e1 = lower_bar(1,1);
e2 = lower_bar(2,1);
e3 = lower_bar(3,1);
e4 = lower_bar(4,1);
e5 = lower_bar(5,1);
e6 = lower_bar(6,1);
X1 = [' A1 = ', num2str(p1), char(177), num2str(e1)];
X2 = [' A2 = ', num2str(p2), char(177), num2str(e2)];
X3 = [' Eg = ', num2str(p3), char(177), num2str(e3)];
X4 = [' Eb = ', num2str(p4), char(177), num2str(e4)];
X5 = [' R = ', num2str(p5), char(177), num2str(e5)];
X6 = [' g = ', num2str(p6), char(177), num2str(e6)];
disp(X1);
disp(X2);
disp(X3);
disp(X4);
disp(X5);
disp(X6);
% Table of parameter values
parameter_name = {'A1'; 'A2'; 'Eg'; 'Eb'; 'R'; 'g'};
parameter_values = [p1; p2; p3; p4; p5; p6];
error = [e1; e2; e3; e4; e5; e6];
T = table(parameter_values, error, 'RowNames',parameter_name);
disp(T)
%% Plot command
plot(E_p, Abs, 'o', 'Color','blue') % scatter plot
hold on
plot(E_p, EM_SS_wR(p, E_p), 'LineWidth', 1.0, 'Color','red') % curve fitting
plot(E_p, EM_SS_wR_CC(p, E_p), 'LineWidth', 1.0,'Color','black') % carrier contribution
plot(E_p, EM_SS_wR_EC(p, E_p), 'LineWidth', 1.0, 'Color','green') % excitonic contribution
% Table of parameter values
TString = evalc('disp(T)');
TString = strrep(TString,'<strong>','\bf');
TString = strrep(TString,'</strong>','\rm');
TString = strrep(TString,'_','\_');
FixedWidth = get(0,'FixedWidthFontName');
annotation(gcf,'Textbox','String',TString,'Interpreter','tex',...
'FontName',FixedWidth,'Units','Normalized','Position',[0.15, 0.8, 0.1, 0.1]);
xlabel('Probe energy (eV)')
ylabel('Absorbance.O.D')
legend('Experimental Data', 'Fitted Curve', 'Carrier Contribution','Excitonic Contribution', 'Location','southeast')
hold off
toc
Hi all, as seen I have a code for curve fitting, and everything is going well but at the last step of graph plotting, I am supposed to get 3 lines but the figure only shows 2. Any idea what is going wrong?
