run_text.m

计算在光纤长度一定时

function x=run_text
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%解沿掺Tm光纤长度方向各能级粒子数分布
%%%%%%%%%%%%%%%%%%%%%%%
N=5.51*10^24;
bbap=5*10^(-25);bbas=0.1*10^(-25);bbes=2.5*10^(-25);
%Taup=0.8964;
%Taup=0.01;
Taus=0.7520;
ap=0.00035;as=0.00035;

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

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

Pin=1;
P(1)=1;

for i=1:n
    
    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;
        %P(i)=Pin*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.1;
        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.1;
        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
    end
%figure(1);plot(z,P,'-*',z,S1,'->',z,S2,'-<')  
hold on
%plot(z,P,'-*')
figure(2);plot(z,N1,z,N2,z,N4)