www.pudn.com > EDFA-980nm.zip > main_980q_dan_ns.m, change:2014-11-04,size:4083b


clear 
clc 
% 所有未注明的值的单位为标准国际单位 
% 调取函数 VUA01 VUA11 afs axs 
% 常用参数 国际标准单位 
a=2.2e-6;NA=0.24;w=500;h0=6.62606896e-34;v=4e12; 
c=2.99793e8;r=linspace(0,a,w); 
bp=980e-9;vp=c./bp; 
% 各种光的输入功率 
ps01=0.01e-3;pa01=0; 
pp01=100e-3;pp11=0e-3; 
Ps=ps01;Pa=pa01;Pp=[pp01 pp11]; 
% 光纤长度和粒子浓度 
L=30; 
n0=zeros(1,w); 
for i=1:w 
%     n0(i)=5e24*(r(i)<=5e-6)+23e24*(r(i)>5e-6);%%%二分阶跃型浓度分布 
%     n0(i)=36*1e24*exp(-((r(i)-5.5e-6)./ 1e-6).^2)+1e24;%%%高斯型浓度分布 
    n0(i)=1e24;%%%均匀型浓度分布 
end 
% % % % % % 调取函数 VUA01 VUA11 
%计算泵浦光的归一化光场分布 
[vp01,Up01,Ap01]= VUA01(a,bp,NA); 
[vp11,Up11,Ap11]= VUA11(a,bp,NA); 
Fp01= Ap01*(besselj(0,Up01*r/a)/besselj(0,Up01)).^2; 
Fp11= Ap11*(besselj(1,Up11*r/a)/besselj(1,Up11)).^2; 
Fp=[Fp01;Fp11]; 
% % % % % % 调取函数 VUA01 
% 信号光的归一化光场分布 
sl=1525;su=1625; 
bs=(sl:1:su).*1e-9;bbs=bs*1e9;%%%bbs单位为nm (方便画图) 
vs=c./bs;oo=length(bs); 
vs01=zeros(1,oo);Us01=zeros(1,oo); 
As01=zeros(1,oo);Fs01=zeros(oo,w); 
for j=1:oo 
    [vs01(j),Us01(j),As01(j)]= VUA01(a,bs(j),NA); 
    Fs01(j,:)= As01(j)*(besselj(0,Us01(j)*r/a)/besselj(0,Us01(j))).^2; 
end 
% % % % % % 调用函数 afs axs 
% 吸收和辐射截面 
oap=2.86e-25; 
oes=afs(bs); 
oas=axs(bs); 
 
% 计算开始 
nc=100;dh=L/nc; 
t=100;%上能级粒子寿命的倒数,单位hz 
yps01=zeros(oo,nc+1);ypa01=zeros(oo,nc+1); 
ypp01=zeros(1,nc+1);ypp11=zeros(1,nc+1); 
yps01(:,1)=Ps;ypa01(:,1)=Pa; 
ypp01(1)=Pp(1);ypp11(1)=Pp(2); 
n1=zeros(nc,w);n2=zeros(nc,w); 
as=0;ap=0; 
mp=zeros(1,w);ms0=zeros(oo,w);ms1=zeros(oo,w);ms2=zeros(oo,w); 
for x=1:nc 
    mp=(ypp01(1,x)*Fp(1,:)+ypp11(1,x)*Fp(2,:))*(oap./(h0.*vp)); 
    for y=1:oo 
        ms0(y,:)=(yps01(y,x)+ypa01(y,x))*Fs01(y,:); 
        ms1(y,:)=(oes(y)/(h0*vs(y))).*ms0(y,:); 
        ms2(y,:)=((oes(y)+oas(y))/(h0*vs(y))).*ms0(y,:); 
    end 
    ms3=sum(ms1,1);ms4=sum(ms2,1); 
    n1(x,:)=n0.*(t+ms3)./(t+ms4+mp);n2(x,:)=n0-n1(x,:); 
    for z=1:oo 
        %计算信号的第一个模式LP01s的功率 
        sss01=fs01(a,w,n1(x,:),n2(x,:),oes(z),oas(z),Fs01(z,:));        %调用子函数fs01 
        s1s1=yps01(z,x).*(sss01-as); 
        s1s2=(yps01(z,x)+(dh/2).*s1s1).*(sss01-as); 
        s1s3=(yps01(z,x)+(dh/2).*s1s2).*(sss01-as); 
        s1s4=(yps01(z,x)+dh.*s1s3).*(sss01-as); 
        yps01(z,x+1)=yps01(z,x)+(dh/6).*(s1s1+2*s1s2+2*s1s3+s1s4); 
        
        %计算ASE的第1个模式LP01a的功率 
        aaa01=fa01(a,w,n2(x,:),oes(z),Fs01(z,:),vs(z));        %调用子函数fa01 
        a1a1=ypa01(z,x).*(sss01-as)+aaa01; 
        a1a2=(ypa01(z,x)+(dh/2).*a1a1).*(sss01-as)+aaa01; 
        a1a3=(ypa01(z,x)+(dh/2).*a1a2).*(sss01-as)+aaa01; 
        a1a4=(ypa01(z,x)+dh.*a1a3).*(sss01-as)+aaa01; 
        ypa01(z,x+1)=ypa01(z,x)+(dh/6).*(a1a1+2*a1a2+2*a1a3+a1a4); 
         
    end 
    %计算泵浦的第1个模式LP01p的功率 
    ppp01=fp01(a,w,n1(x,:),Fp(1,:),oap);        %调用子函数fp01 
    p1p1=ypp01(x).*(ppp01-ap); 
    p1p2=(ypp01(x)+(dh/2).*p1p1).*(ppp01-ap); 
    p1p3=(ypp01(x)+(dh/2).*p1p2).*(ppp01-ap); 
    p1p4=(ypp01(x)+(dh).*p1p3).*(ppp01-ap); 
    ypp01(x+1)=ypp01(x)+(dh/6).*(p1p1+2*p1p2+2*p1p3+p1p4); 
    
    %计算泵浦的第2个模式LP11p的功率 
    ppp11=fp11(a,w,n1(x,:),Fp(2,:),oap);        %调用子函数fp11 
    p4p1=ypp11(x).*(ppp11-ap); 
    p4p2=(ypp11(x)+(dh/2).*p4p1).*(ppp11-ap); 
    p4p3=(ypp11(x)+(dh/2).*p4p2).*(ppp11-ap); 
    p4p4=(ypp11(x)+dh.*p4p3).*(ppp11-ap); 
    ypp11(x+1)=ypp11(x)+(dh/6).*(p4p1+2*p4p2+2*p4p3+p4p4); 
end 
 
Gs01=(yps01(:,end)./yps01(:,1))'; 
ya01=(ypa01(:,end))'; 
s01=10*log10(Gs01); 
NF01=10*log10(1./Gs01+ya01./(h0.*vs.*Gs01.*v)); 
 
% 画图 
figure(1) 
set(gca,'FontName','Times New Roman','FontSize',14); 
plot(bbs,s01,bbs,NF01,'LineWidth',2); 
dd=legend('信号LP01','信号NF01'); 
set(dd,'FontSize',10); 
xlim([sl su]) 
set(gca,'xtick',sl:10:su) 
title(['掺杂分布, L=' num2str(L) 'm'],'FontName','Times New Roman','fontsize',18,'FontWeight','demi'); 
xlabel('波长(nm)','FontName','Times New Roman','fontsize',18,'FontWeight','demi'); 
ylabel('增益(dB)','FontName','Times New Roman','fontsize',18,'FontWeight','demi');