run_text_p.m

计算在光纤长度一定时

function s=run_text_P
N=5.51*10^24;
bbap=5*10^(-25);bbas=0.1*10^(-24);bbes=2.5*10^(-25);
Taup=0.8964;
%Taup=0.01;
Taus=0.7520;
ap=0.00035;as=0.7;

L=1.5;
R4=0.9;R3=1;
a_mirror=-log(R3*R4)/(2*L);
%a_mirror=0;

Pin=0:0.05:1;
m=length(Pin);

for jj=1:m

z=0:0.01:0.4;
n=length(z);
P=zeros(n,1);
S1=zeros(n,1);
S2=zeros(n,1);




for i=1:n
    
    P(i)=Pin(jj);
    
    x=newton(P(i),S1(i),S2(i));
    N1(i)=x(1);
    N2(i)=x(2);
    N4(i)=N-x(1)-x(2);
    
    %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
    
    Ps01=[0.1,0.1/R3];
    tn=z(i)+0.001;
    tspan=[0,tn];

        P(1)=Pin(jj);
        %P(i)=Pin(jj)*exp(-(Taup*(bbap*N1(i))+ap)*z(i));
        [v,Pp]=ode45('eq2',tspan,P(i),[],N1(i));
        u=length(v);
        P(i)=Pp(u);
        %P(i)=sum(Pp)/u;
        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*N2(i)-bbas*N1(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;
        %S2(1)=0.001;
        S2(i)=S2(i)*exp(-(Taus*(bbes*N2(i)-bbas*N1(i))+as-a_mirror)*z(i));

        if i<n
            S2(i+1)=S2(i);
        else
            S2(i)=S2(i);
        end
        S(jj)=S1(n)+S2(n)-0.1;
        %%%%%%%%%%%%%
        if S(jj)<0
            S(jj)=0;
        else
            S(jj)=S(jj);
        end
end
    
end
hold on
plot(Pin,S,'-*')