run_fone_p.m

来自「计算Tm」· M 代码 · 共 96 行

function x=run_fone_P
Ntm=6.18*10^(25);
%Nho=6.18*10^(24);
Nho=Ntm/10;
bbap=1.2*10^(-24);bbep=5*10^(-25);bbas=2*10^(-25);bbes=3*10^(-25); 
%Taup=0.8964;0.0043
Taup=0.01;
Taus=0.7520;
ap=0.00035;as=0;
%ap=0;as=0;
L=3;
R4=0.9;R3=1;
a_mirror=-log(R3*R4)/(2*L);
h=6.63*10^(-34);c=3*10^8;A=1.39*10^(-11);Lp=790*10^(-9);Ls=2*10^(-6);
n0=1.5;
P0=2*h*(c/n0)^1.33/(Ls^3);


Pin=0:0.01:1;
m=length(Pin)
for j=1:m


z=0:0.1:3;
n=length(z);
P=zeros(n,1);
S1=zeros(n,1);
S2=zeros(n,1);

P(1)=Pin(j);




for i=1:n
    
    x=fone(P(i),S1(i),S2(i));
    N2(i)=x(1);
    N4(i)=x(2);
    N1(i)=Ntm-x(1)-x(2);
    N6(i)=x(3);
    N8(i)=x(4);
    N9(i)=x(5);
    N5(i)=Nho-x(3)-x(4)-x(5);
    
        
        
        P(i)=Pin(j)*exp(-(Taup*(bbap*N1(i))+ap)*z(i));
        if i<n
            P(i+1)=P(i);
        else
            P(i)=P(i);
        end
        
        S1(1)=0.001;
        S1(i)=S1(i)*exp((Taus*(bbes*N6(i)-bbas*N5(i))-as-a_mirror)*z(i));
        if i<n
            S1(i+1)=S1(i);
        else
            S1(i)=S1(i);
        end
        
        S2(i)=S1(i)*R4;             %%%%% Attation
        S2(1)=0.001;
        S2(i)=S2(i)*exp(-(Taus*(bbes*N6(i)-bbas*N5(i))-as-a_mirror)*z(i));

        if i<n
            S2(i+1)=S2(i);
        else
            S2(i)=S2(i);
        end
        

        
    
    if (Taus*(bbes*N6(i)-bbas*N5(i))-as-a_mirror)<0
        Ps(i)=0;
        g(i)=0;
    else
        
        g(i)=Taus*(bbes*N6(i)-bbas*N5(i))-as-a_mirror;
    end
end
S=S1+S2;

B(j)=S(n);

end
Ps=B(1:m);
for ii=1:m
    SL(ii)=Ps(ii)/(Pin(ii)-0.15);
end
hold on
%plot(Pin,Ps,'-*')
plot(Pin,SL,'-*')