auto-searching completed program.m

计算稀土离子Tm3+

%% 可以自动辨别前向跳变点，给出前向抽运阈值。
n=7.88*10^20;
m=8*10^19;
segmaTm12=1*10^-23;
segmaTm23=3.54*10^-21;
wTm21=227;
wTm31=2000;
wTm32=304;
wHo21=170;
alafa=0;
alafaTmHo=0;
gamaTmHo=1*10^-17;
gamaHoTm=3.9*10^-19;
beta=3.1*10^-17;
c=0.5;
h=6.63*10^-34;
mu=3.0*10^8/(1.06*10^-6);
X0=[n;0;m];                 % difine the initial value for iteration
i=0;

Pn=40000;            % intensity range given
dh=200;              % space step in iteration
Num=Pn/dh+2;
i=0;
i_f=Num;                    % used to identify the divergent point                 
N1f=zeros(Num-1,1);
N2f=zeros(Num-1,1);
M1f=zeros(Num-1,1);      % matrix for calculation results storage

%% circling calculation by "for"
for P=0:dh:Pn     % threshold maximum point
  
      if (i+1)<i_f     % Jump out of the circuling calculation at divergent point
        i=i+1;         % circulating number
        % difine function group matrix and Jacobi matrix  
        yita=P/(h*mu);    
        f1=sym('yita*segmaTm23*n2-beta*n1*(n-n1-n2)-(wTm32+wTm31)*(n-n1-n2)+alafa*n2^2');
        f2=sym('yita*segmaTm12*n1-yita*segmaTm23*n2+2*beta*n1*(n-n1-n2)+wTm32*(n-n1-n2)-wTm21*n2-gamaTmHo*n2*m1+gamaHoTm*n1*(m-m1)-2*alafa*n2^2-alafaTmHo*n2*(m-m1)+c*alafaTmHo*n2*(m-m1)');
        f3=sym('-wHo21*(m-m1)+gamaTmHo*n2*m1-gamaHoTm*n1*(m-m1)-c*alafaTmHo*n2*(m-m1)');      
        FM=[f1;f2;f3]; 
        DF11=diff(f1,'n1');DF12=diff(f1,'n2');DF13=diff(f1,'m1');
        DF21=diff(f2,'n1');DF22=diff(f2,'n2');DF23=diff(f2,'m1');
        DF31=diff(f3,'n1');DF32=diff(f3,'n2');DF33=diff(f3,'m1');    
        DFJ=[DF11,DF12,DF13;                        
             DF21,DF22,DF23;
             DF31,DF32,DF33];     
        n1=X0(1);
        n2=X0(2);
        m1=X0(3);     % initial value of levels concentration  
        F=subs(FM);
        DF=subs(DFJ);
        DFF=DF\F;      % with real matrix elements
        X=X0-DFF;
        a=norm(X-X0);
        b=norm(X0);
        c=a/b;
        d=0;  
        
        % Newton iteration for every specific "P"
        while c>10^-2&d<1             % judge at the given accuracy with relative error and identify the convergence or divergence of iteration sequence
             
            X0=X; 
             X1=X0;
             n1=X0(1);
             n2=X0(2);
             m1=X0(3);         %former value as the input value of later calculation
             F=subs(FM);
             DF=subs(DFJ);
             DFF=DF\F;          
             X=X0-DFF;
             X2=X;
             d1=norm(X2-X1);         % former space of the iteration sequence
              
             X0=X; 
             n1=X0(1);
             n2=X0(2);
             m1=X0(3);            %former value as the input value of later calculation
             F=subs(FM);
             DF=subs(DFJ);
             DFF=DF\F;          
             X=X0-DFF;
             X3=X;
             d2=norm(X3-X2);       % later space of the iteration sequence      
             d=norm(d2/d1);        % the ratio of later and former sequence
             a=norm(X-X0);
             b=norm(X0);
             c=a/b;   
        end
        
        if c<=10^-2
           N1f(i)=X(1);
           N2f(i)=X(2);
           M1f(i)=X(3);     % producing matrix elements for every levle concentration
           X0=X;           % former value as the input value of later calculation
       else
           disp('>>****************************************************************')
           disp('     A divergent sequence occurs as increasing pumping density ! ')
           disp('     The below corresponds to the first divergent point ')                  
           disp('>>****************************************************************')
           i_f=i;
           Pmax=P-dh;
           IfP=[i_f Pmax]              % show the first divergent point and forward pumping threshold
       end
       
   end

end

N3f=n*ones(Num-1,1)-N1f-N2f;
M2f=m*ones(Num-1,1)-M1f;
N1f;N2f;N3f;M1f;M2f;
 
%% 可以自动辨别逆向跳变点，给出逆向抽运阈值。 
n1P=1.78*10^20;
n2P=6*10^20;
m1P=0.01*10^19;              % given initial values of level concentration for power intensity much more than threshold
X0=[n1P;n2P;m1P];             % difine the initial value for iteration
Pnb=40000;
i=Pnb/dh+2;
i_b=0;
N1b=zeros(Num-1,1);
N2b=zeros(Num-1,1);
M1b=zeros(Num-1,1);      % matrix for calculation results storage

%% circling calculation by "for"
for P=Pnb:-dh:0
    
    if (i-1)>=i_b
       
        i=i-1;
       % difine function group matrix and Jacobi matrix  
       yita=P/(h*mu);    
       f1=sym('yita*segmaTm23*n2-beta*n1*(n-n1-n2)-(wTm32+wTm31)*(n-n1-n2)+alafa*n2^2');
       f2=sym('yita*segmaTm12*n1-yita*segmaTm23*n2+2*beta*n1*(n-n1-n2)+wTm32*(n-n1-n2)-wTm21*n2-gamaTmHo*n2*m1+gamaHoTm*n1*(m-m1)-2*alafa*n2^2-alafaTmHo*n2*(m-m1)+c*alafaTmHo*n2*(m-m1)');
       f3=sym('-wHo21*(m-m1)+gamaTmHo*n2*m1-gamaHoTm*n1*(m-m1)-c*alafaTmHo*n2*(m-m1)');      
       FM=[f1;f2;f3]; 
       DF11=diff(f1,'n1');DF12=diff(f1,'n2');DF13=diff(f1,'m1');
       DF21=diff(f2,'n1');DF22=diff(f2,'n2');DF23=diff(f2,'m1');
       DF31=diff(f3,'n1');DF32=diff(f3,'n2');DF33=diff(f3,'m1');    
       DFJ=[DF11,DF12,DF13;                        
            DF21,DF22,DF23;
            DF31,DF32,DF33];       
       n1=X0(1);
       n2=X0(2);
       m1=X0(3);    % initial value of levels concentration  
       F=subs(FM);
       DF=subs(DFJ);
       DFF=DF\F; % with real matrix elements
       X=X0-DFF;
       a=norm(X-X0);
       b=norm(X0);
       c=a/b;
       d=0;
   
      % Newton iteration for every specific "P"
      while c>10^(-2)&d<1             % judge at the given accuracy with relative error
            X0=X;
            X1=X0;
            n1=X0(1);
            n2=X0(2);
            m1=X0(3); %former value as the input value of later calculation
            F=subs(FM);
            DF=subs(DFJ);
            DFF=DF\F;          
            X=X0-DFF;
            X2=X;
            d1=norm(X2-X1);         % former space of the iteration sequence
                
            X0=X; 
            n1=X0(1);
            n2=X0(2);
            m1=X0(3);            %former value as the input value of later calculation
            F=subs(FM);
            DF=subs(DFJ);
            DFF=DF\F;          
            X=X0-DFF;
            X3=X;
            d2=norm(X3-X2);       % later space of the iteration sequence      
            d=norm(d2/d1);        % the ratio of later and former sequence
            a=norm(X-X0);
            b=norm(X0);
            c=a/b;    
       end
       
         if c<=10^-2
            N1b(i)=X(1);
            N2b(i)=X(2);
            M1b(i)=X(3);     % producing matrix elements for every levle concentration
            X0=X;           % former value as the input value of later calculation
         else
            disp('>>****************************************************************')
            disp('     A divergent sequence occurs as decreasing pumping density ! ')
            disp('     The below corresponds to the first divergent point ')                  
            disp('>>****************************************************************')
            i_b=i;
            Pmin=P+dh;
            IbP=[i_b Pmin]              % show the first divergent point and forward pumping threshold
         end
    end
end
N3b=n*ones(Num-1,1)-N1b-N2b;
M2b=m*ones(Num-1,1)-M1b;
N1b;N2b;N3b;M1b;M2b;
 
% 自动整合前向和逆向数据点
a=1:i_b;
N1b(a)=N1f(a);
N2b(a)=N2f(a);
M1b(a)=M1f(a);
M2b(a)=M2f(a);
b=i_f:Num-1;
N1f(b)=N1b(b);
N2f(b)=N2b(b);
M1f(b)=M1b(b);
M2f(b)=M2b(b);

% show figures
P=0:dh:Pn;
figure;
hold on;
plot(P,N1f,'k');plot(P,N2f,'r');plot(P,10*M1f,'b');plot(P,10*M2f,'g');
plot(P,N1b,'k');plot(P,N2b,'r');plot(P,10*M1b,'b');plot(P,10*M2b,'g');
 
% 保存数据
save c:\matlab6p1\work\SSL-IOB\fig3data-f N1f N2f N3f M1f M2f   % save the forward results calculated
save c:\matlab6p1\work\SSL-IOB\fig3data-b N1b N2b N3b M1b M2b   % save the backward results calculated