光纤模式本征值的计算及仿真

光纤模式本征值的计算及仿真
可以实现βlm的计算并绘出模场电场磁场的强度分布以及矢量分布
matlab代码如下：

clear;
close all;
format long
%constant
global n1 n2 a lambda k0 A B C D omega mu epsilon1 epsilon2
n1 = 1.46;
n2 = 1.44;
a  =  4.5;%um
lambda = 1.550;%um
k0 = 2*pi/lambda;
omega = 2*pi*3*10^14/1.55;%(um/s)/um
mu    = 4*pi*10^-13;%V·s/(A·um)
epsilon1 = n1^2*8.854187817 *10^-18;% F/um
epsilon2 = n2^2*8.854187817 *10^-18;% F/um
%main
m = 0;
for l = 0:1:1%取HE11和TE01模
    jugde_val = 0;
    prev_judge_val = 0;
    for beta = n1*k0-0.000001: -0.000001 : n2*k0+0.000001
        [judge_val,U, W ,A,B,C,D] = eigenEquationJudge(l,beta);
        if(judge_val * prev_judge_val < 0)
           m=m+1;
           eigen_val(m)  = 0.5*(beta + prev_beta);       
           X = []; Y = [];Z = [];Ex = [];Ey = [];Ez = [];Hx = [];Hy = [];Hz = [];
            for r = 0:0.5:6
                for phi = 0:(0.1*pi):2*pi
                    Ez   = [Ez real(cal_Ez(l,r,phi,U,W))];
                    Er   = real(cal_Er  (l,r,phi,eigen_val(m),U,W));
                    Ephi = real(cal_Ephi(l,r,phi,eigen_val(m),U,W));
                    Hz   = [Hz real(cal_Hz(l,r,phi,U,W))];
                    Hr   = real(cal_Hr  (l,r,phi,eigen_val(m),U,W));
                    Hphi = real(cal_Hphi(l,r,phi,eigen_val(m),U,W));
                    X = [X r*cos(phi)];
                    Y = [Y r*sin(phi)];
                    Z = [Z 0];
                    Ex =[Ex Er*cos(phi)-Ephi*sin(phi)];
                    Ey =[Ey Er*sin(phi)+Ephi*cos(phi)];
                    Hx =[Hx Hr*cos(phi)-Hphi*sin(phi)];
                    Hy =[Hy Hr*sin(phi)+Hphi*cos(phi)];
                end 
            end
            figure
            quiver3(X,Y,Z,Ex,Ey,Ez);
            hold on
            quiver3(X,Y,Z,Hx,Hy,Hz);
            axis equal
            for t = 1:1:length(Ex)
            IntensityE(t) = (Ex(t)^2+Ey(t)^2+Ez(t)^2)^0.5;
            end
            figure
            [x,y,z]=griddata(X',Y',IntensityE',linspace(min(X),max(X))',linspace(min(Y),max(Y)),'v4');%插值
            pcolor(x,y,z);
            shading interp
            for t = 1:1:length(Hx)
            IntensityH(t) = (Hx(t)^2+Hy(t)^2+Hz(t)^2)^0.5;
            end
            figure
            [x,y,z]=griddata(X',Y',IntensityH',linspace(min(X),max(X))',linspace(min(Y),max(Y)),'v4');%插值
            pcolor(x,y,z);
            shading interp
            %xlabel('x'),ylabel('y'),zlabel('z');
            break;
        end
        prev_judge_val = judge_val;
        prev_beta      = beta;
    end
end
figure(1);
title("TE01模电磁场矢量分布");
legend("E","H");
saveas(gca,'1.jpg');
figure(2);
title("TE01模电场强度分布");
saveas(gca,'2.jpg');
figure(3);
title("TE01模磁场强度分布");
saveas(gca,'3.jpg');
figure(4);
title("HE11模电磁场矢量分布");
legend("E","H");
saveas(gca,'4.jpg');
figure(5);
title("HE11模电场强度分布");
saveas(gca,'5.jpg');
figure(6);
title("HE11模磁场强度分布");
saveas(gca,'6.jpg');
function [result,U,W,A,B,C,D] = eigenEquationJudge(l,beta)
   global n1 n2 a lambda k0 omega mu
   U = a*(n1^2*k0^2-beta^2)^0.5;
   W = a*(beta^2-n2^2*k0^2)^0.5;  
   left1  = ( besselj(l-1,U)-besselj(l+1,U))/(2*U*besselj(l,U))           ...
          + (-besselk(l-1,W)-besselk(l+1,W))/(2*W*besselk(l,W));
   left2  = (n1^2*k0^2*( besselj(l-1,U)-besselj(l+1,U))/(2*U*besselj(l,U)) ...
           + n2^2*k0^2*(-besselk(l-1,W)-besselk(l+1,W))/(2*W*besselk(l,W)));
   if(l==0&&U>2.405)
        if(abs(left1)<abs(left2))
        %TE模本征值方程
        A = 0;
        B = 1i;%补偿相位，防止实部为零
        C = 0;
        D = B*besselj(l,U)/besselk(l,W);
        result = left1;
        else
        %TM模本征值方程
        A = 1;
        B = 0;
        C = A*besselj(l,U)/besselk(l,W);
        D = 0;
        result = left2 ;
        end
   else
       right1 = l*beta*(1/U^2+1/W^2);
       A = 1;
       B = 1i*right1/(omega*mu*left1)*A;
       C = A*besselj(l,U)/besselk(l,W);
       D = B*besselj(l,U)/besselk(l,W);
       result = left1*left2 - right1^2;
   end
  
end
function Ez = cal_Ez(l,r,phi,U,W)
global n1 n2 a lambda k0 A B C D omega mu
        if(r>=0&&r<=10)
            Ez = A*besselj(l,U/a*r)*exp(1i*l*phi);
        else
            Ez = C*besselk(l,W/a*r)*exp(1i*l*phi);
        end
end
function Er = cal_Er(l,r,phi,beta,U,W)
global n1 n2 a lambda k0 A B C D omega mu
        if(r>=0&&r<=a)
            Er = -1i*(a/U)^2*exp(1i*l*phi)*...
                ((A*beta*(U/a)     *( besselj(l-1,U*r/a)-besselj(l+1,U*r/a))/2)...
                +(B*1i*omega*mu*l/(r+eps)*  besselj(l,  U*(r+eps)/a)));                      
                
        else
            Er =  1i*(a/W)^2*exp(1i*l*phi)*...
                ((C*beta*(W/a)      *(-besselk(l-1,W*r/a)-besselk(l+1,W*r/a))/2)...
                +(D*1i*omega*mu*l/r *  besselk(l,  U*r/a)));
        end
end
function Ephi = cal_Ephi(l,r,phi,beta,U,W)
global n1 n2 a lambda k0 A B C D omega mu
        if(r>=0&&r<=a)
                 Ephi =-1i*(a/U)^2*exp(1i*l*phi)*...
                 (A*1i*beta*l/(r+eps)*besselj(l,U*(r+eps)/a)...
                 -B*omega*mu*(U/a)   *( besselj(l-1,U*r/a)-besselj(l+1,U*r/a))/2);
        else
                 Ephi = 1i*(a/W)^2*exp(1i*l*phi)*...
                 (C*beta*l/r       *besselk(l,W*r/a)...
                 -D*omega*mu*(W/a)   *(-besselk(l-1,W*r/a)-besselk(l+1,W*r/a))/2);
        end
end
function Hz = cal_Hz(l,r,phi,U,W)
global n1 n2 a lambda k0 A B C D omega epsilon1 epsilon2
        if(r>=0&&r<=10)
            Hz = B*besselj(l,U/a*r)*exp(1i*l*phi);
        else
            Hz = D*besselk(l,W/a*r)*exp(1i*l*phi);
        end
end
function Hphi = cal_Hphi(l,r,phi,beta,U,W)
global n1 n2 a lambda k0 A B C D omega epsilon1 epsilon2
        if(r>=0&&r<=a)
            Hphi = -1i*(a/U)^2*exp(1i*l*phi)*...
                ((A*omega*epsilon1*(U/a)     *( besselj(l-1,U*r/a)-besselj(l+1,U*r/a))/2)...
                +(B*1i*beta*l/(r+eps)*  besselj(l,  U*(r+eps)/a)));                      
                
        else
            Hphi =  1i*(a/W)^2*exp(1i*l*phi)*...
                ((C*omega*epsilon2*(W/a)      *(-besselk(l-1,W*r/a)-besselk(l+1,W*r/a))/2)...
                +(D*1i*beta*l/r      *  besselk(l,  U*r/a)));
        end
end
function Hr = cal_Hr(l,r,phi,beta,U,W)
global n1 n2 a lambda k0 A B C D omega epsilon1 epsilon2
        if(r>=0&&r<=a)
                 Hr =-1i*(a/U)^2*exp(1i*l*phi)*...
                (-A*1i*omega*epsilon1*l/(r+eps)*besselj(l,U*(r+eps)/a)...
                 +B*beta*(U/a)      *( besselj(l-1,U*r/a)-besselj(l+1,U*r/a))/2);
        else
                 Hr = 1i*(a/W)^2*exp(1i*l*phi)*...
                 (C*omega*epsilon2*l/r       *besselk(l,W*r/a)...
                 +D*beta*(W/a)   *(-besselk(l-1,W*r/a)-besselk(l+1,W*r/a))/2);
        end
end