# I need help to make graphs please

1 vue (au cours des 30 derniers jours)
Yordani le 3 Jan 2023
Commenté : Yordani le 3 Jan 2023
Hello! I have the following program 1 that depends on subroutines 2,3 and 4 I need to graph [flux(index) vs w] and also [Q_t(j) vs d] can someone help me? The graphs should be on a logarithmic scale on both axes.
program 1)
clear all
clc
tic
h = 1.054571596e-34;
global c0;
c0 = 2.99792458e+8;
kb = 1.3806503e-23;
qe = 1.602176462e-19;
e0 = 8.854187817e-12;
T1 = 300;
T2 = 299;
dc = 0.5e-9;
beta_max = pi/dc;
darray = [1e-8, 1e-7, 1e-6];
w1 = 2e12;
w2 = 2e15;
dw = 2e12;
nw = floor((w2-w1)/dw+1);
for j = 1:length(darray)
d=darray(j);
index=0;
for w=w1:dw:w2;
ep1=Lorentz_SiC(w);
ep2=Lorentz_SiC(w);
index = index+1;
a0 = 0;
a1 = (w-10)/c0;
nkp = 1000;
errp = 0.01;
ae1 = (w+10)/c0;
ae2 = 8*w/c0;
ae3 = 100*w/c0;
ae4 = beta_max;
nke0 = 1000;
nke1 = 1000;
nke2 = 1000;
erre0 = 1.0;
erre1 = 0.01;
erre2 = 0.1;
Q1=h*w/(exp(h*w/(kb*T1))-1);
Q2=h*w/(exp(h*w/(kb*T2))-1);
val_p = simpson_p(w, nkp, d, a0, a1, ep1, ep2, errp);
fluxp=val_p*(Q1-Q2);
val_e0= simpson_e(w, nke0, d, ae1, ae2, ep1, ep2, erre0);
fluxe0=val_e0*(Q1-Q2);
val_e1= simpson_e(w, nke1, d, ae2, ae3, ep1, ep2, erre1);
fluxe1=val_e1*(Q1-Q2);
val_e2= simpson_e(w, nke2, d, ae3, ae4, ep1, ep2, erre2);
fluxe2=val_e2*(Q1-Q2);
flux(j,index)=fluxp+fluxe0+fluxe1+fluxe2;
end
Q_t(j)=(flux(1)+4*sum(flux(2:2:(nw-1)))+2*sum(flux(3:2:(nw-2)))+flux(nw))*dw/3;
end
figure(1)
plot(darray,Q_t)
figure(2)
plot(w1:dw:w2,flux)
-----------------------------------------------------------------------------------------------------------------
subroutine 2
function [e_1] = Lorentz_SiC(w)
wL = 1.82652e+14;
wT = 1.49477e+14;
gamma = 8.9724e+11;
e_inf = 6.7;
e_1 = e_inf*(1+(wL*wL-wT*wT)./(wT*wT-w.*w-i*gamma*w));%1-wp*wp/(w*w+i*w*gal);
end
-------------------------------------------------------------------------------------------------------------------
subroutine 3
function [Int_val] = simpson_p(w, n, d, min, max, ep_1, ep_2, err)
n1 = n;
[s_p,dkx] = func_p(w, n1, d, min, max, ep_1, ep_2);
temp1=(s_p(1)+4*sum(s_p(2:2:n1))+2*sum(s_p(3:2:n1-1))+s_p(n1+1))*dkx/3;
if(temp1<=1e-50)
temp2 = 0;
else
n1 = n1*2;
[s_p,dkx] = func_p(w, n1, d, min, max, ep_1, ep_2);
temp2=(s_p(1)+4*sum(s_p(2:2:n1))+2*sum(s_p(3:2:n1-1))+s_p(n1+1))*dkx/3;
while ((abs(temp2-temp1)/temp2) >= err)
temp1 = temp2;
n1 = n1*2;
[s_p,dkx] = func_p(w, n1, d, min, max, ep_1, ep_2);
temp2=(s_p(1)+4*sum(s_p(2:2:n1))+2*sum(s_p(3:2:n1-1))+s_p(n1+1))*dkx/3;
end
end
Int_val = temp2;
end
function [s_p,dkx] = func_p(w, n, d, min, max, ep_1, ep_2)
global c0
dkx = (max-min)/n;
kx = zeros(n+1,1);
s_p = zeros(n+1,1);
for ind=1:n+1
kx(ind)= min+(ind-1)*dkx;
kz1 = sqrt(ep_1*w*w/(c0*c0)-kx(ind)^2);
kz2 = sqrt(ep_2*w*w/(c0*c0)-kx(ind)^2);
kz3 = sqrt(w*w/(c0*c0)-kx(ind)^2);
rs31 = ((kz3-kz1)/(kz3+kz1));
rp31 = ((ep_1*kz3-kz1)/(ep_1*kz3+kz1));
rs32 = ((kz3-kz2)/(kz3+kz2));
rp32 = ((ep_2*kz3-kz2)/(ep_2*kz3+kz2));
tps_temp = abs(1-rs31*rs32*exp(2*i*kz3*d));
tpp_temp = abs(1-rp31*rp32*exp(2*i*kz3*d));
s_p(ind) = kx(ind)*((1-abs(rs31)^2)*(1-abs(rs32)^2)/tps_temp^2+(1-abs(rp31)^2)*(1-abs(rp32)^2)/(tpp_temp^2))/(4*pi*pi);
end
end
----------------------------------------------------------------------------------------------------------
subroutine 4
function [Int_val] = simpson_e(w, n, d, min, max, ep_1, ep_2, err)
n1 = n;
[s_e,dkx] = func_e(w, n1, d, min, max, ep_1, ep_2);
temp1=(s_e(1)+4*sum(s_e(2:2:n1))+2*sum(s_e(3:2:n1-1))+s_e(n1+1))*dkx/3;
if(temp1<=1e-50)
temp2 = 0;
else
n1 = n1*2;
[s_e,dkx] = func_e(w, n1, d, min, max, ep_1, ep_2);
temp2=(s_e(1)+4*sum(s_e(2:2:n1))+2*sum(s_e(3:2:n1-1))+s_e(n1+1))*dkx/3;
while ((abs(temp2-temp1)/temp2) >= err)
temp1 = temp2;
n1 = n1*2;
[s_e,dkx] = func_e(w, n1, d, min, max, ep_1, ep_2);
temp2=(s_e(1)+4*sum(s_e(2:2:n1))+2*sum(s_e(3:2:n1-1))+s_e(n1+1))*dkx/3;
end
end
Int_val = temp2;
end
function [s_e,dkx] = func_e(w, n, d, min, max, ep_1, ep_2)
global c0
dkx = (max-min)/n;
kx = zeros(n+1,1);
s_e = zeros(n+1,1);
for ind=1:n+1
kx(ind)= min+(ind-1)*dkx;
kz1 = sqrt(ep_1*w*w/(c0*c0)-kx(ind)^2);
kz2 = sqrt(ep_2*w*w/(c0*c0)-kx(ind)^2);
kz3 = sqrt(w*w/(c0*c0)-kx(ind)^2);
rs31 = ((kz3-kz1)/(kz3+kz1));
rp31 = ((ep_1*kz3-kz1)/(ep_1*kz3+kz1));
rs32 = ((kz3-kz2)/(kz3+kz2));
rp32 = ((ep_2*kz3-kz2)/(ep_2*kz3+kz2));
e_temp = exp(-2*imag(kz3)*d);
s_e(ind) = kx(ind)*e_temp*(imag(rs31)*imag(rs32)/(abs(1-rs31*rs32*e_temp))^2+imag(rp31)*imag(rp32)/(abs(1-rp31*rp32*e_temp))^2)/(pi*pi);
end
end
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Réponse acceptée

Torsten le 3 Jan 2023
See above. I added four lines to your code to plot the results and changed the line
flux(index)=fluxp+fluxe0+fluxe1+fluxe2;
to
flux(j,index)=fluxp+fluxe0+fluxe1+fluxe2;
##### 1 commentaireAfficher -1 commentaires plus anciensMasquer -1 commentaires plus anciens
Yordani le 3 Jan 2023
Thnaks!

Connectez-vous pour commenter.

### Plus de réponses (1)

It is easier to help for specific questions, i.e., you tried something and it did not work. The way you have phrased this is more like trying to solve your homework and you will not get much help.
##### 1 commentaireAfficher -1 commentaires plus anciensMasquer -1 commentaires plus anciens
Yordani le 3 Jan 2023
OK thank you!

Connectez-vous pour commenter.

### Catégories

En savoir plus sur 2-D and 3-D Plots dans Help Center et File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by