duichenmochang.m

用MATLAB编写的光子晶体光纤中模场的分布是准二维的

clc;
clear;
warning off;

mu0=4*pi*10^(-7);  %真空磁导率
epsilon0=8.85*10^(-12);    %真空介电常数
n_slab=1.45; 
n_air=1;
n_slab2=n_slab^2;
n_air2=n_air^1;
n_average=(n_slab2+n_air2)/2;
epsilon_slab=n_slab2*epsilon0;
epsilon_air=n_air2*epsilon0;
epsilon_average=(n_slab2+n_air2)*epsilon0/2;
%z_slab2=mu0/epsilon_slab;
z_air2=mu0/epsilon0;

pitch=1.8e-6;         % 空洞中心的间距 
ddx=0.7e-6/13;               %空间步长
kz=15/pitch;    %  定义纵向传播常数
dt=0.95/(3e8*sqrt((2/ddx^2)+(kz/2)^2)); %定义时间步长

a=ceil(pitch/ddx);    % 空洞中心的间距的网格数
%pitchs=1.5e-6;          %内层缺陷孔周期
%as=ceil(pitchs/ddx);  
%ds=1.4e-6;
d1=pitch*0.3;   
d2=pitch*0.3; 
d3=pitch*0.3; 
d4=pitch*0.3;
%rs=ceil(ds/2/ddx); 
r1=ceil(d1/2/ddx);   % x/y/z 相当于x/(y*z)
r2=ceil(d2/2/ddx);
r3=ceil(d3/2/ddx);
r4=ceil(d4/2/ddx);

layer=4;       %调节包层空气孔层数
npml=8;              % PML层数

I=(layer+1)*2*a+2*npml; % 总网格数
J=I;
ic=I/2;jc=J/2;        % I,J必须是偶数，ic,jc,kc才符合下标要求

spread=round(0.8e-6/6e8/dt);            % 信源宽度 对时间归一化
t0=round(4*spread+0.5)-1;                %脉冲中心  对时间归一化
r0=100;

%---------设置网格电磁参数

para=zeros(I,J,6);  %第一个介电常数 第二个磁导系数 第三、四个电导率(x,y) 第五、六个导磁率(x,y)
xn_max=1/(30*n_slab*pi*ddx);  % xn相当于电导率  介电常数取硅的数值
for i=1:I;
    for j=1:J;
        para(i,j,2)=mu0;
        para(i,j,1)=epsilon_slab; %先全部置成硅的介电常数
    end
end

%设置PML边界
for i=1:npml;
    xnum=npml-i+1;
    xxn=xnum/npml;
    xn=xn_max*(xxn^4); % 定义电导率 指数为4
    xm=z_air2*xn;  % 定义导磁率                         
    for j=1:J;
        para(i,j,3)=xn;
        para(i,j,5)=xm;
        para(I-i+1,j,3)=xn;
        para(I-i+1,j,5)=xm;
        para(j,i,4)=xn;
        para(j,i,6)=xm;
        para(j,I-i+1,4)=xn;
        para(j,I-i+1,6)=xm;
    end
end
        
%设置光子晶体光纤周期参数

%第一列：wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
for j=jc;
  for i=jc;   %for i=jc:a:jc+a;
    for i1=-r1:r1;
      for j1=-r1:r1;
        if sqrt(i1^2+j1^2)<r1;
               para(round(i+i1),round(j+j1),1)=epsilon_air;
           end
       end
   end
     for i1=-r1-1:r1+1;
      for j1=-r1-1:r1+1;
          if para(i+i1,j+j1,1)==epsilon_air & (para(i+i1-1,j+j1,1)==epsilon_slab | para(i+i1+1,j+j1,1)==epsilon_slab | para(i+i1,j+j1-1,1)==epsilon_slab | para(i+i1,j+j1+1,1)==epsilon_slab);
                 para(i+i1,j+j1,1)=epsilon_average;
             end
         end
     end  %设置阶梯近似
 end
end

for j=jc;%r2=3;
      for i=jc+2*a;
        for i1=-r2:r2;
         for j1=-r2:r2;
            if(sqrt(i1^2+j1^2)<r1);
                para(round(i+i1),round(j+j1),1)=n_air2*epsilon0;
            end
         end
      end   
      for i1=-r2:r2;
         for j1=-r2:r2;
            if para(i+i1,j+j1,1)==epsilon_air & (para(i+i1-1,j+j1,1)==epsilon_slab | para(i+i1+1,j+j1,1)==epsilon_slab | para(i+i1,j+j1-1,1)==epsilon_slab | para(i+i1,j+j1+1,1)==epsilon_slab);
                 para(round(i+i1),round(j+j1),1)=epsilon_average;
         end
      end 
   end  %设置阶梯近似   
end
end

for j=jc;%r3=4;
    for i=jc+3*a;
       for i1=-r3:r3;
         for j1=-r3:r3;
            if(sqrt(i1^2+j1^2)<r2);
               para(round(i+i1),round(j+j1),1)=n_air2*epsilon0;
            end
         end
      end
      for i1=-r3:r3;
        for j1=-r3:r3;
      if para(i+i1,j+j1,1)==epsilon_air & (para(i+i1-1,j+j1,1)==epsilon_slab | para(i+i1+1,j+j1,1)==epsilon_slab | para(i+i1,j+j1-1,1)==epsilon_slab | para(i+i1,j+j1+1,1)==epsilon_slab);
          para(round(i+i1),round(j+j1),1)=epsilon_average;
        end
      end 
   end %设置阶梯近似   
end
end

for j=jc;%r4=5
   for i=jc+4*a;
      for i1=-r4:r4;
         for j1=-r4:r4;
            if(sqrt(i1^2+j1^2)<r2);
                para(round(i+i1),round(j+j1),1)=n_air2*epsilon0;
            end
         end
       end  
      for i1=-r4:r4;
         for j1=-r4:r4;
          if para(i+i1,j+j1,1)==epsilon_air & (para(i+i1-1,j+j1,1)==epsilon_slab | para(i+i1+1,j+j1,1)==epsilon_slab | para(i+i1,j+j1-1,1)==epsilon_slab | para(i+i1,j+j1+1,1)==epsilon_slab);
               para(round(i+i1),round(j+j1),1)=epsilon_average;
          end
      end 
   end %设置阶梯近似  
end
end

%第二列：wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
w1=sqrt(3)/2*a;w2=a/2;
for j=jc+w1;%r1=2;
   for i=jc+w2;
      for i1=-r1:r1;
         for j1=-r1:r1;
            if(sqrt(i1^2+j1^2)<r1);
                para(round(i+i1),round(j+j1),1)=n_air2*epsilon0;
            end
         end
      end   
      for i1=-r1-1:r1+1;
         for j1=-r1-1:r1+1;
           if para(round(i+i1),round(j+j1),1)==epsilon_air & (para(round(i+i1-1),round(j+j1),1)==epsilon_slab | para(round(i+i1+1),round(j+j1),1)==epsilon_slab | para(round(i+i1),round(j+j1-1),1)==epsilon_slab | para(round(i+i1),round(j+j1+1),1)==epsilon_slab);
               para(round(i+i1),round(j+j1),1)=epsilon_average;
           end
      end 
  end %设置阶梯近似
end
end

for j=jc+w1;%r2=3;
   for i=jc+w2+a;
      for i1=-r2:r2;
         for j1=-r2:r2;
            if(sqrt(i1^2+j1^2)<r1);
                para(round(i+i1),round(j+j1),1)=n_air2*epsilon0;
            end
         end
       end  
      for i1=-r2:r2;
         for j1=-r2:r2;
            if  para(round(i+i1),round(j+j1),1)==epsilon_air & (para(round(i+i1-1),round(j+j1),1)==epsilon_slab | para(round(i+i1+1),round(j+j1),1)==epsilon_slab | para(round(i+i1),round(j+j1-1),1)==epsilon_slab | para(round(i+i1),round(j+j1+1),1)==epsilon_slab);
               para(round(i+i1),round(j+j1),1)=epsilon_average;
            end
      end
    end %设置阶梯近似  
end
end

for j=jc+w1;%r3=4;
   for i=jc+w2+2*a;
      for i1=-r3:r3;
         for j1=-r3:r3;
            if(sqrt(i1^2+j1^2)<r2);
              para(round(i+i1),round(j+j1),1)=n_air2*epsilon0;
            end
         end
      end  
      for i1=-r3:r3;
         for j1=-r3:r3;
           if  para(round(i+i1),round(j+j1),1)==epsilon_air & (para(round(i+i1-1),round(j+j1),1)==epsilon_slab | para(round(i+i1+1),round(j+j1),1)==epsilon_slab | para(round(i+i1),round(j+j1-1),1)==epsilon_slab | para(round(i+i1),round(j+j1+1),1)==epsilon_slab);
               para(round(i+i1),round(j+j1),1)=epsilon_average;
         end
      end
    end %设置阶梯近似
end
end

for j=jc+w1;%r4=5;
   for i=jc+w2+3*a;
      for i1=-r4:r4;
         for j1=-r4:r4;
            if(sqrt(i1^2+j1^2)<r2);
             para(round(i+i1),round(j+j1),1)=n_air2*epsilon0;
            end
         end
      end   
      for i1=-r4:r4;
         for j1=-r4:r4;