run_fone_n.m

计算Tm

function x=run_fone_N

%Nho=1.05*10^(25);
bbap=1.2*10^(-24);bbep=5*10^(-25);bbas=2*10^(-25);bbes=2.915*10^(-25); 
%Taup=0.8964;0.0043
Taup=0.01;
Taus=0.7520;
ap=0.00035;as=0.00035;
%ap=0;as=0;
L=3;
R4=0.04;R3=0.9;
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);
n=1.5;
P0=2*h*(c/n)^1.33/(Ls^3);

Ntm=1*10^(25):0.1*10^(25):10*10^(25);
m=length(Ntm);

for j=1:m
    Nho(j)=Ntm(j)/10;

z=0:0.1:3.1;
n=length(z);
P=zeros(n,1);
S1=zeros(n,1);
S2=zeros(n,1);
S=S1+S2;
Pin=1;
P(1)=1;




for i=1:n
    
    x=fone_N(Ntm(j),P(i),S(i));
    N2(i)=x(1);
    N4(i)=x(2);
    N1(i)=Ntm(j)-x(1)-x(2);
    N6(i)=x(3);
    N8(i)=x(4);
    N9(i)=x(5);
    N5(i)=Nho(j)-x(3)-x(4)-x(5);
    
        
        
        P(i)=Pin*exp(-(Taup*(bbap*N1(i))+ap)*z(i));
        if i<n
            P(i+1)=P(i);
        else
            P(i)=P(i);
        end
        
        S1(1)=10^(-10);
        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)=10^(-10);
        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
        S(i)=S1(i)+S2(i);

        
    
    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;
Ps(j)=S(n);
Pp(j)=P(n)
end
hold on
%plot(Ntm,Ps,'-x')
plot(Ntm,Pp,'-x')