I need help to make graphs please

Question

Yordani 2023-1-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1887897-i-need-help-to-make-graphs-please

评论： Yordani 2023-1-3

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2023-1-3

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1887897-i-need-help-to-make-graphs-please#answer_1140517

在 MATLAB Online 中打开

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Yordani 2023-1-3

Thnaks!

请先登录，再进行评论。

Answer 2

Constantino Carlos Reyes-Aldasoro 2023-1-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1887897-i-need-help-to-make-graphs-please#answer_1140457

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.