I need help to make graphs please
显示 更早的评论
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
采纳的回答
Torsten
2023-1-3
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 个）
Constantino Carlos Reyes-Aldasoro
2023-1-3
0 个投票
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.
类别
在 帮助中心 和 File Exchange 中查找有关 2-D and 3-D Plots 的更多信息
Translated by 
