Changing the Limits of the Graph When Trying to Model the Shockley Equation in MatLab Shifts the Position of the Curve [closed]

Question

I have been trying to create a piecewise linear model of a diode through curve fitting. While trying to implement the code below, I ran into a strange bug that is causing my output curve to shift to be roughly -0.2, the max value of the graph, as you can see in the picture below. The same problem occurs when the range of vd is shifted or an xlim is added to the graph. If you have any ideas, please let me know.

Thanks!

value :

%% variables
Is=1e-12;   % Saturation current (in Amperes)
n = 2;         % Ideality factor (typically between 1 and 2)
T = 300;       % Temperature in Kelvin
k = 1.380649e-23;  % Boltzmann constant (J/K)
q = 1.602176634e-19; % Electron charge (C)

% Calculate thermal voltage
vt=(k*T)/q;  % Thermal voltage (in Volts)
vd=0:0.01:10;  % Voltage range from -0.1V to 1V
%figure, plot(time,v1);

%% diode model calculations

%calculate the current voltage relation ship for the diode
Id_s=Is.*(exp((vd/(n*vt)))-1); %calcuate diode current using shockley's equation
%plot schokley eq
figure, 
plot(vd, Id_s);  % Use a logarithmic scale for current
xlabel('V_d (V)');  % Label for the x-axis (Voltage)
ylabel('I_d (A)');  % Label for the y-axis (Current)
title('Shockley equation diode model');
grid on;

%% curve fitting
range_1=0.74:0.01:0.88; %define range for the first curve piece
range_2=0.89:0.01:max(vd); %define range for the second curve piece
Id_fit1=polyfit(range_1,Id_s(find(vd==0.74):find(vd==0.88)),1); %fit a line to the first part of the curve
Id_fit2=polyfit(range_2,Id_s(find(vd==0.89):find(vd==max(vd))),1);%fit a line to the second part of the curve

line_1=polyval(Id_fit1,range_1);
line_2=polyval(Id_fit2,range_2);

Id_fit=[zeros(1,length(vd)-length(range_1)-length(range_2)),line_1, line_2]; %create a combine vector of each piece wise component

figure, 
plot(vd, Id_s);  % Use a logarithmic scale for current
xlabel('V_d (V)');  % Label for the x-axis (Voltage)
ylabel('I_d (A)');  % Label for the y-axis (Current)
title('Shockley equation diode model');
grid on;

hold on
plot(vd,Id_fit);

legend('Shockley Curve, Curve Fit')

Shockley curve with vd ranging from 0-2.5

Answer 1

summary of an approach

If you have any ideas, please let me know.

There are just three parameters to deal with:

• saturation current
• emission coefficient (non-ideality factor)
• bulk resistance

You can use part of the curve where the currents are highest to get an idea of the bulk resistance (the slope becomes far more linear there.)

So you really only need three well-placed data points to get some reasonable results. One at a very low current, one at the knee, and one at a high current. Four would be better still: one pair for the low current leg that defines the y-intercept (saturation current) and slope (emission coefficient) and one pair on the other side of the knee (higher current) to lock down the bulk resistance.

quantitative clarification example

Let me put together a simple diode model:

.model MyDiode ako:1N4148 D IS=1n N=2 RS=3

Look at the low current end:

I've added a simple line showing the y-axis intercept, which just happens to also be the saturation current (IS in the model.) Note that this matches the model.

The slope of the curve indicates the emission coefficient. \$V_T\approx 25.865 \:\text{mV}\$ at \$27^\circ\text{C}\$ (which is where I ran the curve.)

To work out the emission coefficient, take the voltage difference of \$718.36735\:\text{mV}-359.18367\:\text{mV}=0.35918368\:\text{mV}\$ and divide that by the logarithm of the current ratios times \$V_T\$:

$$\begin{align*} N &= \frac{718.36735\:\text{mV}-359.18367\:\text{mV}}{25.865\:\text{mV}\cdot \ln\left(\frac{1.012724\:\text{mA}}{1.0352707\:\mu\text{A}}\right)}\approx 2 \end{align*}$$

And that also matches up pretty close.

The only remaining step is the bulk resistance. For that, a lin-lin chart is better:

$$\begin{align*} R_\text{S} &= \frac{1.9953488\:\text{V}-1.5962791\:\text{V}}{327.05057\:\text{mA}-202.30955\:\text{mA}}=3.2\:\Omega \end{align*}$$

That estimate will be off by a tiny error, which is the emission coefficient times \$V_T\$ and divided by the average current for the span used, which is \$264.68\:\text{mA}\$ (halfway between.) So this is about \$\frac{2\,\cdot\, 25.865\:\text{mV}}{264.68\:\text{mA}}\approx 0.195\:\Omega\$.

So subtract that to correct the \$R_\text{S}=3.2\:\Omega-.195\:\Omega\approx 3\:\Omega\$ value.

It doesn't need to be more complicated than that.

summary

Not sure if any of that helps you write your MATLAB code. I'm not familiar with it, or its functions and their behaviors, and so I'm in no position to write your code. But perhaps the above brings up a detail to help you see what to do about the coding.