How to plot a graph correctly?

5 views (last 30 days)
Dmitry
Dmitry on 26 Jan 2023
Commented: Star Strider on 30 Jan 2023
I have a function of two variables F(t,d) and am building a plane, and I want to plot a function where F(t,d) - C_2=0, I use function 'contour', but I get a few curves and if I take contour(t,d,F) or contour(d,t,F) only the location of the coordinate axes changes, but the graph remains the same. In conclusion, you need to get a graph t(d) for F() - C_2=0.
My code:
%% initial conditions
global k0 h_bar ksi m E C_2
Ef = 2.77*10^3;
Kb = physconst('boltzmann'); % 1.38*10^(-23)
m = 9.1093837*10^(-31);
Tc = 1.2;
ksi = 10^(-9);
h_bar = (1.0545726*10^(-34));
k0 = (ksi/h_bar)*sqrt(2.*m.*pi.*Kb.*Tc);
C_2 = 6.81;
t = linspace(0.1, 2);
d = linspace(10, 1000);
%[TT,DD] = meshgrid(t,d);
%contour(TT, DD, @f_calc);
for i=1:numel(d)
for j = 1:numel(t)
F(i,j) = f_calc(t(j),d(i));
end
end
%plot(t,F)
contour(d,t,F)
%xlabel('d')
%ylabel('t')
%surf(t,d,F)
function F = f_calc(t, d)
global k0 h_bar ksi m C_2
nD = floor(375/(2*pi.*t*1.2) - 0.5);
F = 0;
for k = 0:nD
F = F + 1/(2*k+1).*imag(f_lg(k,t,d)+1i.*d.*k0.*((f_p1(k,t)-f_p2(k,t))./2)+1i*f_arg_1(k,t,d)-1i*f_arg_2(k,t,d));
end
F = -(1/d).*F - 7.826922968141167;
end
function p1 = f_p1(n,t)
p1 = ((1+1i)./sqrt(2)).*sqrt(t.*(2.*n+1));
end
function p2 = f_p2(n,t)
p2 = sqrt(3601+1i.*t.*(2.*n+1));
end
function arg_1 = f_arg_1(n,t,d)
global k0
arg_1 = angle(1+exp(-1i.*d*k0.*f_p1(n,t)));
end
function arg_2 = f_arg_2(n,t,d)
global k0
arg_2 = angle(1+exp(-1i.*d*k0.*f_p2(n,t)));
end
function n_lg = f_lg(n,t,d)
global k0;
arg_of_lg = (1+exp(-1i.*d.*k0.*f_p1(n,t)))/(1+exp(-1i.*d.*k0.*f_p2(n,t)));
n_lg = log(abs(arg_of_lg));
end

Answers (2)

Star Strider
Star Strider on 26 Jan 2023
I am not certain that I understand what you want to do.
If you only want the contour at ‘F=0’, supply that as the fourth, ‘Levels’ argument to contour:
contour(d,t,F,[0 0])
Try this —
%% initial conditions
global k0 h_bar ksi m E C_2
Ef = 2.77*10^3;
Kb = physconst('boltzmann'); % 1.38*10^(-23)
m = 9.1093837*10^(-31);
Tc = 1.2;
ksi = 10^(-9);
h_bar = (1.0545726*10^(-34));
k0 = (ksi/h_bar)*sqrt(2.*m.*pi.*Kb.*Tc);
C_2 = 6.81;
t = linspace(0.1, 2);
d = linspace(10, 1000);
%[TT,DD] = meshgrid(t,d);
%contour(TT, DD, @f_calc);
for i=1:numel(d)
for j = 1:numel(t)
F(i,j) = f_calc(t(j),d(i));
end
end
%plot(t,F)
contour(d,t,F,[0 0])
xlabel('d')
ylabel('t')
%surf(t,d,F)
function F = f_calc(t, d)
global k0 h_bar ksi m C_2
nD = floor(375/(2*pi.*t*1.2) - 0.5);
F = 0;
for k = 0:nD
F = F + 1/(2*k+1).*imag(f_lg(k,t,d)+1i.*d.*k0.*((f_p1(k,t)-f_p2(k,t))./2)+1i*f_arg_1(k,t,d)-1i*f_arg_2(k,t,d));
end
F = -(1/d).*F - 7.826922968141167;
end
function p1 = f_p1(n,t)
p1 = ((1+1i)./sqrt(2)).*sqrt(t.*(2.*n+1));
end
function p2 = f_p2(n,t)
p2 = sqrt(3601+1i.*t.*(2.*n+1));
end
function arg_1 = f_arg_1(n,t,d)
global k0
arg_1 = angle(1+exp(-1i.*d*k0.*f_p1(n,t)));
end
function arg_2 = f_arg_2(n,t,d)
global k0
arg_2 = angle(1+exp(-1i.*d*k0.*f_p2(n,t)));
end
function n_lg = f_lg(n,t,d)
global k0;
arg_of_lg = (1+exp(-1i.*d.*k0.*f_p1(n,t)))/(1+exp(-1i.*d.*k0.*f_p2(n,t)));
n_lg = log(abs(arg_of_lg));
end
.
  16 Comments
Dmitry
Dmitry on 30 Jan 2023
Star Strider, are these all contours for the zero level?
Star Strider
Star Strider on 30 Jan 2023
As far as I can tell, yes.
For whatever reason, contour occasionally breaks up a contour at one level into several different contours. Sometimes this is obvious, for example taking the contours at 0 of the peaks function, although here it is less obvious, at least to me, since I do not understand what you are doing. There are apparently non-zero regions between the zero segments here, and that breaks the 0 contour into disjointed segments.
You would likely have to look at a small section of the contour plot near the 0 region and plot all the available contours in that region to see what it is doing. (It may be necessary to specify those contours as described in the documentation section on levels. Use the xlim function to restrict the region of the contour plot, for example [400 600].)
I leave this to you because I do not understand what your code calculates, so I cannot interpret that result.

Sign in to comment.


Torsten
Torsten on 26 Jan 2023
Try "fimplicit".
  3 Comments
Torsten
Torsten on 26 Jan 2023
Edited: Torsten on 26 Jan 2023
Sure. At least for the ranges of t and d where we searched, there was no solution pair for the equation F(t,d) - C_2 = 0. Did you find one in the meantime ?
Dmitry
Dmitry on 26 Jan 2023
Edited: Walter Roberson on 29 Jan 2023
If we look at this solution, we can see that there is a line but if you plot contour(d,t,F,[0 0]) orcontour(d,t,F,[0 0]), do they not change and that's weird. I changed the constant that we subtract from F and now there are solutions there.
%% initial conditions
global k0 h_bar ksi m E C_2
Ef = 2.77*10^3;
Kb = physconst('boltzmann'); % 1.38*10^(-23)
m = 9.1093837*10^(-31);
Tc = 1.2;
ksi = 10^(-9);
h_bar = (1.0545726*10^(-34));
k0 = (ksi/h_bar)*sqrt(2.*m.*pi.*Kb.*Tc);
C_2 = 6.81;
t = linspace(0.1, 2);
d = linspace(10, 1000);
%[TT,DD] = meshgrid(t,d);
%contour(TT, DD, @f_calc);
for i=1:numel(d)
for j = 1:numel(t)
F(i,j) = f_calc(t(j),d(i));
end
end
%plot(t,F)
contour(d,t,F,[0 0])
xlabel('d')
ylabel('t')
%surf(t,d,F)
function F = f_calc(t, d)
global k0 h_bar ksi m C_2
nD = floor(375/(2*pi.*t*1.2) - 0.5);
F = 0;
for k = 0:nD
F = F + 1/(2*k+1).*imag(f_lg(k,t,d)+1i.*d.*k0.*((f_p1(k,t)-f_p2(k,t))./2)+1i*f_arg_1(k,t,d)-1i*f_arg_2(k,t,d));
end
F = -(1/d).*F - 7.826922968141167;
end
function p1 = f_p1(n,t)
p1 = ((1+1i)./sqrt(2)).*sqrt(t.*(2.*n+1));
end
function p2 = f_p2(n,t)
p2 = sqrt(3601+1i.*t.*(2.*n+1));
end
function arg_1 = f_arg_1(n,t,d)
global k0
arg_1 = angle(1+exp(-1i.*d*k0.*f_p1(n,t)));
end
function arg_2 = f_arg_2(n,t,d)
global k0
arg_2 = angle(1+exp(-1i.*d*k0.*f_p2(n,t)));
end
function n_lg = f_lg(n,t,d)
global k0;
arg_of_lg = (1+exp(-1i.*d.*k0.*f_p1(n,t)))/(1+exp(-1i.*d.*k0.*f_p2(n,t)));
n_lg = log(abs(arg_of_lg));
end

Sign in to comment.

Categories

Find more on Contour Plots in Help Center and File Exchange

Products

Community Treasure Hunt

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

Start Hunting!