EDFA-980nm.zip_980单模光纤_EDFA matlab_掺铒_掺铒光纤

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');