ouheqi.rar_2×2耦合器_coupler_fiber_光纤_光纤耦合

%基本参数定义 m=0; n=1; n_clad=1.4955; wavelength=1.55*1e-06; k0=2*pi/wavelength; z0=122*pi; %波阻抗 a1_0=1; %输入边界条件 a2_0=0; eption_0=1/36/pi*1e-9; %真空介电常数 omega=2*pi/wavelength*3*1e8; %光波角频率 %定义1号光纤的参数 r1_core=5*1e-06; r1_clad=3*r1_core; nc1=1.5; delta1=(nc1-n_clad)/nc1; V1=nc1*k0*r1_core*sqrt(2*delta1); %判断1号光纤是否满足单模条件 if V1>2.40483,error('1号光纤不是单模光纤');end %由特征方程求1号光纤的归一化径向相位常数和归一化径向衰减常数 fun=@(U1) U1*besselj(m-1,U1)/besselj(m,U1)+sqrt(V1^2-U1^2)*besselk(m-1,sqrt(V1^2-U1^2))/besselk(m,sqrt(V1^2-U1^2)); U1=fzero(fun,1); W1=sqrt(V1^2-U1^2); %求解1号光纤的传播常数beta beta1=sqrt((U1/r1_clad)^2-nc1^2*k0^2); %求解1号光纤的A A1=U1/V1*besselk(m,W1)/sqrt(besselk(m-1,W1)*besselk(m+1,W1))*sqrt(4*z0/n_clad/pi/r1_core^2); %定义2号光纤的参数 r2_core=5*1e-06; r2_clad=3*r2_core; nc2=1.5; delta2=(nc2-n_clad)/nc2; V2=nc2*k0*r2_core*sqrt(2*delta2); %判断1号光纤是否满足单模条件 if V2>2.40483,error('2号光纤不是单模光纤');end %由特征方程求2号光纤的归一化径向相位常数和归一化径向衰减常数 fun=@(U2) U2*besselj(m+1,U2)/besselj(m,U2)-sqrt(V2^2-U2^2)*besselk(m+1,sqrt(V2^2-U2^2))/besselk(m,sqrt(V2^2-U2^2)); U2=fzero(fun,1); W2=sqrt(V2^2-U2^2); %求解2号光纤的传播常数beta beta2=sqrt((U2/r2_clad)^2-nc2^2*k0^2); %求解2号光纤的A A2=U2/V2*besselk(m,W2)/sqrt(besselk(m-1,W2)*besselk(m+1,W2))*sqrt(4*z0/n_clad/pi/r2_core^2); %求解耦合系数k12,k21 d=r1_clad+r2_clad; %两光纤轴向间隔 betacha=beta1-beta2; %相位失配常数 F=@(r,theta) conj(A1*cos(m*theta).*besselj(m,U1*r/r1_core)./besselj(m,U1)).*(A2*cos(m*theta).*besselk(m,W2*sqrt(d^2+r.^2-2*d*r.*cos(theta))/r2_core)/besselk(m,W2))*(nc1^2-n_clad^2).*r*omega*eption_0/4; k12=dblquad(F,0,r1_core,0,2*pi); %积分求得k12 k21=conj(k12); %二者共轭 c2=((betacha+sqrt(betacha^2+4*abs(k21)^2))*a1_0+2*k12*a2_0)/(2*sqrt(betacha^2+4*abs(k21)^2)); c1=a1_0-c2; %积分常数c1,c2 %耦合求解 %1号光纤 r1=linspace(-r1_core,r1_core,400); %rz平面上纤芯内某点的纵坐标 z1=linspace(0,8*1e-0,1000/1); %耦合距离 [Z1,R1]=meshgrid(z1,r1); theta0=0; %沿零度方位角取光纤的轴向剖面（rz平面） rr1=linspace(r1_core,(d+r2_clad),600); %rz包层内某点的纵坐标 [ZZ1,RR1]=meshgrid(z1,rr1); a1=(c1*exp(1i*(betacha+sqrt(betacha^2+4*abs(k21)^2))/2*z1)+c2*exp(1i*(betacha-sqrt(betacha^2+4*abs(k21)^2))/2*z1)); %1号光纤中场振幅 Score1=A1^2/2/z0*(cos(m*theta0))^2*nc1.*(besselj(m,U1*abs(r1)'/r1_core)./besselj(m,U1)).^2*abs(a1).^2; %纤芯功率密度 Sclad1=A1^2/2/z0*(cos(m*theta0))^2*n_clad.*(besselk(m,W1*rr1'/r1_core)/besselk(m,W1)).^2*abs(a1).^2; %包层功率密度 %2号光纤 r2=linspace(d-r2_core,d+r2_core,400); z2=linspace(0,8*1e-0,1000/1); [Z2,R2]=meshgrid(z2,r2); rr2=linspace(0,d-r2_core,400); [ZZ2,RR2]=meshgrid(z2,rr2); a2=1i/k12*(c1*1i*(betacha+sqrt(betacha^2+4*abs(k21)^2))/2*exp(1i*(-betacha+sqrt(betacha^2+4*abs(k21)^2))/2*z2)+c2*1i*(betacha-sqrt(betacha^2+4*abs(k21)^2))/2*exp(-1i*(betacha+sqrt(betacha^2+4*abs(k21)^2))/2*z2)); %2号光纤中场振幅 Score2=A2^2/2/z0*(cos(m*theta0))^2*nc2.*(besselj(m,U2*abs(r2-d)'/r2_core)./besselj(m,U2)).^2*abs(a2).^2; %纤芯功率密度 Sclad2=A2^2/2/z0*(cos(m*theta0))^2*n_clad.*(besselk(m,W2*(d-rr2)'/r2_core)/besselk(m,W2)).^2*abs(a2).^2; %包层功率密度 %耦合功率绘图 cmap=[linspace(0,1,256);zeros(1,256);zeros(1,256)]';colormap(cmap); colormap(cmap); %color map('default') %color map(gray(256)) subplot(2,1,2); pcolor(Z1,R1,Score1); hold on; pcolor(ZZ1,RR1,Sclad1); shading flat;colorbar;axis tight; xlabel('耦合距离z(m)');ylabel('r','Fontsize',13,'Fontname','Times'); title('1号光纤中的功率密度分布'); subplot(2,1,1); pcolor(Z2,R2,Score2); hold on; pcolor(ZZ2,RR2,Sclad2); shading flat;colorbar;axis tight; xlabel('耦合距离z（m）');ylabel('r','Fontsize',13,'FontName','Times'); title('2号光纤的功率密度分布'); %耦合效率 figure Pout_1=abs(a1).^2; Pout_2=abs(a2).^2; yita=Pout_2; subplot(3,1,1); plot(z1,Pout_1);xlabel('耦合距离z(m)');ylabel('直通臂功率Pout_1'); subplot(3,1,2); plot(z2,Pout_2);xlabel('耦合距离z(m)');ylabel('耦合臂功率Pout_2'); subplot(3,1,3); plot(z1,yita);xlabel('耦合距离z(m)');ylabel('耦合效率');