fone.m

计算Tm

function x=fone(P,S1,S2)
tau2=334.7*10^(-6);tau4=14.2*10^(-6);tau42=0.72*tau4;tau6=567*10^(-6);tau8=1.4*10^(-6);tau9=0.022*10^(-6);tau41=0.28*tau4;
%Taup=0.8964;tau6=567*10^(-6)
Taup=0.01;
Taus=0.7520;
h=6.63*10^(-34);c=3*10^8;A=1.39*10^(-11);Lp=790*10^(-9);Ls=2*10^(-6);
k4212=5*10^(-25);k2124=0.084*k4212;         % Tm交叉弛豫系数
k2156=3.3*10^(-21);k1265=k2156/12.54;       % Tm-Ho能量传递系数
k2168=1.0*10^(-25);k1286=0.1*k2168;         % 上转换1
k4169=1*10^(-25);k1496=0.1*k4169;           % 上转换2
bbap=1.2*10^(-24);bbep=5*10^(-25);bbas=2*10^(-24);bbes=3*10^(-24);  % 受激吸收与辐射截面

Ntm=6.18*10^(25);
%Nho=1.05*10^(25);
%Nho=6.18*10^(24);
Nho=Ntm/10;
ap=0.00035;as=0.00035;L=1;
n=1.45;
P0=2*h*(c/n)^2/(Ls^3);

R4=0.9;R3=1;
%P=1;


maxit=100;

w14=Taup*Lp*bbap*P/(h*c*A);
w56=Taus*Ls*bbas*(S1+S2)/(h*c*A);
w65=Taus*Ls*bbes*(S1+S2)/(h*c*A);
x=zeros(5,1);f=zeros(5,1);
for k=1:maxit
    f(1)=(1/tau42+2*k4212*Ntm)*x(2)-(1/tau2+k2156*Nho)*x(1)+2*k4212*x(1)*x(2)-2*k4212*x(2)^2-2*k2124*x(1)^2+(k2156-k1265-k2168)*x(1)*x(3)+(k2156-k1286)*x(1)*x(4)+k2156*x(1)*x(5)+k1265*Ntm*x(3)-k1265*x(2)*x(3)+k1286*Ntm*x(4)-k1286*x(2)*x(4);
    f(2)=w14*Ntm-w14*x(1)-(w14+1/tau4+k4212*Ntm)*x(2)+k4212*x(1)*x(2)+k4212*x(2)^2+k2124*x(1)^2-k4169*x(2)*x(3)+k1496*Ntm*x(5)-k1496*x(1)*x(5)-k1496*x(2)*x(5);
    f(3)=(-1/tau6-k1265*Ntm-w56-w65)*x(3)+(k1286*Ntm-w56)*x(4)+(k1496*Ntm-w56)*x(5)+k2156*Nho*x(1)+(k1265-k2156-k2168)*x(1)*x(3)-(k2156+k2168)*x(1)*x(4)-(k2156+k1496)*x(1)*x(5)+(k1265-k4169)*x(2)*x(3)-k1286*x(2)*x(4)-k1496*x(2)*x(5)+w56*Nho;
    f(4)=-(1/tau8+k1286*Ntm)*x(4)+k2168*x(1)*x(3)+k1286*x(1)*x(4)+k1286*x(2)*x(4);
    f(5)=-(1/tau9+k1496*Ntm)*x(5)+k4169*x(2)*x(3)+k1496*x(1)*x(5)+k1496*x(2)*x(5);
    J=[-(1/tau2+k2156*Nho)+2*k4212*x(2)-4*k2124*x(1)+(k2156-k1265-k2168)*x(3)+(k2156-k1286)*x(4)+k2156*x(5) 1/tau42+2*k4212*Ntm+2*k4212*x(1)-4*k4212*x(2)-k1265*x(3)-k1286*x(4) (k2156-k1265-k2168)*x(1)+k1265*Ntm-k1265*x(2) (k2156-k1286)*x(1)+k1286*Ntm-k1286*x(2) k2156*x(1)
        -w14+k4212*x(2)+2*k2124*x(1)-k1496*x(5) -(w14+1/tau4+k4212*Ntm)+k4212*x(1)+2*k4212*x(2)-k4169*x(3)-k1496*x(5) -k4169*x(2) 0 k1496*(Ntm-x(1)-x(2))
        k2156*Nho+(k1265-k2156-k2168)*x(3)-(k2156+k1286)*x(4)-(k2156+k1496)*x(5) (k1265-k4169)*x(3)-k1286*x(4)-k1496*x(5) (-1/tau6-k1265*Ntm-w56-w65)+(k1265-k2156-k2168)*x(1)+(k1265-k4169)*x(2) (k1286*Ntm-w56)-(k2156+k1286)*x(1)-k1286*x(2) (k1496*Ntm-w56)-(k2156+k1496)*x(1)-k1496*x(2)
        k2168*x(3)+k1286*x(4) k1286*x(4) k2168*x(1) -(1/tau8+k1286*Ntm)+k1286*(x(1)+x(2)) 0
        k1496*x(5) k4169*x(3)+k1496*x(5) k4169*x(2) 0 -(1/tau9+k1496*Ntm)+k1496*(x(1)+x(2))];
    dx=-J\f;
    x=x+dx;
end