Clear Filters
Clear Filters

Graph not showing as it should

2 views (last 30 days)
Jack
Jack on 20 Sep 2024
Commented: Jack on 20 Sep 2024

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?

Accepted Answer

Epsilon
Epsilon on 20 Sep 2024
Edited: Epsilon on 20 Sep 2024
Hi Jack,
It seems the Excitonic Contribution line is missing but it is being plotted. It is just not visible due to being plotted parallel to the horizontal axis. On replacing it with circle markers, it becomes visible.
Original vs Replaced with circle markers:
This is the due to the returned value of the ‘defined function EM_SS_wR_EC’ being a scalar, consider using the following definition instead:
% Excitonic Contribution
function F = EM_SS_wR_EC(p, E_p)
F = zeros(size(E_p));
% Compute excitonic contribution
for n = 1:7
F = F + (1/n^3) * sech((E_p - p(3) + p(4)/n^2) ./ p(6));
end
% Apply scaling factor
F = p(2) * ((2 * pi * p(4)^3 / 2) ./ E_p) * (1 / p(6)) .* F;
end
Plot with replaced function definition:
Glad to help!

More Answers (0)

Categories

Find more on Mathematics in Help Center and File Exchange

Products


Release

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!