mat3.m

基于MATLAB

h=6.626e-34;
vp=3.0612e14;
vs=1.9595e14;
p=0.58;
s=0.53;
%n0=2e25;
g(1)=1;
ac=3e-12;
ps=1e-3;
pl=1e-3;
t=10e-3;
ap=2.58e-25;
as=5.36e-25;
es=5.41e-25;
b=p*ap/(as+es)/s;
%c=b*s*es*n0;
d=vp*ps/vs;
%e=s*as*n0;
m=ac*h*vp/t/s/(as+es);

hhh=1;
for pp=0:0.001:0.15
n=1;
n01=8e25
z1=0.3;
c=b*s*es*n01
e=s*as*n01;
fg(n)=g(n).^b*exp(-c*z1)*pp-pp+d*(g(n)-1)+(log(g(n))+e*z1)*m;
fg1(n)=b*g(n).^(b-1)*exp(-c*z1)*pp+d+m/g(n);
g(n+1)=g(n)-fg(n)/fg1(n);
gg=g(n);
while g(n+1)-g(n)>0.0001
n=n+1;
fg(n)=g(n).^b*exp(-c*z1)*pp-pp+d*(g(n)-1)+(log(g(n))+e*z1)*m;
fg1(n)=b*g(n).^(b-1)*exp(-c*z1)*pp+d+m/g(n);
g(n+1)=g(n)-fg(n)/fg1(n);
gg=g(n);
end
n02=8e25;
z2=0.5;
c=b*s*es*n02;
e=s*as*n02;
fg(n)=g(n).^b*exp(-c*z2)*pp-pp+d*(g(n)-1)+(log(g(n))+e*z2)*m;
fg1(n)=b*g(n).^(b-1)*exp(-c*z2)*pp+d+m/g(n);
g(n+1)=g(n)-fg(n)/fg1(n);
yy=g(n);
while g(n+1)-g(n)>0.0001
n=n+1;
fg(n)=g(n).^b*exp(-c*z2)*pp-pp+d*(g(n)-1)+(log(g(n))+e*z2)*m;
fg1(n)=b*g(n).^(b-1)*exp(-c*z2)*pp+d+m/g(n);
g(n+1)=g(n)-fg(n)/fg1(n);
yy=g(n);
end
n03=8e25;
z3=0.8;
c=b*s*es*n03;
e=s*as*n03;
fg(n)=g(n).^b*exp(-c*z3)*pp-pp+d*(g(n)-1)+(log(g(n))+e*z3)*m;
fg1(n)=b*g(n).^(b-1)*exp(-c*z3)*pp+d+m/g(n);
g(n+1)=g(n)-fg(n)/fg1(n);
zz=g(n);
while g(n+1)-g(n)>0.0001
n=n+1;
fg(n)=g(n).^b*exp(-c*z3)*pp-pp+d*(g(n)-1)+(log(g(n))+e*z3)*m;
fg1(n)=b*g(n).^(b-1)*exp(-c*z3)*pp+d+m/g(n);
g(n+1)=g(n)-fg(n)/fg1(n);
zz=g(n);
end
ggg(hhh)=10*log10(gg);
yyy(hhh)=10*log10(yy);
zzz(hhh)=10*log10(zz);
P=ps*exp(ggg*z1/10);
P=ps*exp(yyy*z2/10);
P=ps*exp(zzz*z3/10)
hhh=hhh+1;
end
pp=0:0.001:0.15
plot(pp,ggg,pp,yyy,pp,zzz)
xlabel('泵浦功率(W)','FontSize',12)
ylabel('增益(dB)','FontSize',12)