duichenmochang.m

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

            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

%第三列：wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
ww1=sqrt(3)/2*a*2;
for j=jc+ww1;%r2=3;
   for i=jc:a:jc+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+ww1;%r3=4;
   for i=jc+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+ww1;%r4=5;
   for i=jc+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;
          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

%第四列：wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
www1=sqrt(3)/2*a*3;
for j=jc+www1;%r3=4;
   for i=jc+w2:a:jc+w2+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+www1;%r4=5;
   for i=jc+w2+2*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(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

%第五列：wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
wwww1=sqrt(3)/2*a*4;
for j=jc+wwww1;%r4=5;
   for i=jc:a:jc+2*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(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

%将以上定义的1/4结构复制到全部平面
for i=1:ic;
   for j=1:jc;
       para(ic-i+1,jc-j+1,1)=para(ic+i-1,jc+j-1,1); % 3象限
	   para(ic-i+1,jc+j-1,1)=para(ic+i-1,jc+j-1,1); % 2象限
       para(ic+i-1,jc-j+1,1)=para(ic+i-1,jc+j-1,1); % 4象限
   end
end

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

debug=1;

hz=zeros(I,J); hx=hz; hy=hz;
ex=hz; ey=hz; ez=hz;
dx=hz; dy=hz; dz=hz; bx=hz; by=hz; bz=hz;
dx_d=hz; dy_d=hz; dz_d=hz; bx_d=hz; by_d=hz; bz_d=hz;


%选取保留时域波形的点 

ey20k=zeros(I,J);
ex20k=zeros(I,J);

% 准备主循环wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww        
nsteps=100;%ceil(sqrt(2)*I*15);

for t=1:nsteps;
    %wwwwwwwwwwwwwwwwwwwwwwwwwwww定义计算所需时间wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
    if t==1
       t_start=cputime;
       disp(strcat('1 / ',int2str(nsteps),'  Calculating ,And trying to calculate the run time...'))
    else
       disp(strcat(int2str(t),' /  ',int2str(nsteps),'  Calculating...May be ','  ',num2str((nsteps-t+1)*t_elapsed),'  seconds left'))

    end
    %wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww定义完成wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww    
   
    %加入初始脉冲
    pulse=exp(-0.5*((t-t0)/spread)^2);    
    
    for i=1:I-1;   %  计算空间各处H值
        for j=1:J-1;
                  tempy_add=1/dt+para(i,j,6)/mu0/2;
                  tempx_add=1/dt+para(i,j,5)/mu0/2;
                  tempy_sub=para(i,j,6)/mu0/2-1/dt;
                  tempx_sub=para(i,j,5)/mu0/2-1/dt;
                  bx_d(i,j)=bx(i,j);
                  bx(i,j)=bx(i,j)-dt*((ez(i,j+1)-ez(i,j))/ddx-kz*ey(i,j));   
                  hx(i,j)=-tempy_sub/tempy_add*hx(i,j)+1/mu0/tempy_add*(tempx_add*bx(i,j)+tempx_sub*bx_d(i,j))+exp(-((i-ic-a)^2+(j-jc)^2)/r0^2)*pulse+exp(-((i-ic+a)^2+(j-jc)^2)/r0^2)*pulse;%在Hx(ic,jc)加激励源
                  by_d(i,j)=by(i,j);
                  by(i,j)=by(i,j)+dt*((ez(i+1,j)-ez(i,j))/ddx-kz*ex(i,j));
                  hy(i,j)=-tempx_sub/tempx_add*hy(i,j)+1/mu0/tempx_add*(tempy_add*by(i,j)+tempy_sub*by_d(i,j));%+(i==ic)*(j==jc)*pulse;%在Hy(ic,jc)加激励源
                  bz_d(i,j)=bz(i,j);
                  bz(i,j)=-tempy_sub/tempy_add*bz(i,j)+1/tempy_add*((ex(i,j+1)-ex(i,j))-(ey(i+1,j)-ey(i,j)))/ddx;
                  hz(i,j)=-tempx_sub/tempx_add*hz(i,j)+1/mu0/tempx_add/dt*(bz(i,j)-bz_d(i,j));
        end
    end      
    
    for i=2:I;  %  计算空间各处E值
        for j=2:J;
                 tempy_add=1/dt+para(i,j,4)/epsilon0/2;
                 tempx_add=1/dt+para(i,j,3)/epsilon0/2;
                 tempy_sub=para(i,j,4)/epsilon0/2-1/dt;
                 tempx_sub=para(i,j,3)/epsilon0/2-1/dt; 
                 dx_d(i,j)=dx(i,j);
                 dx(i,j)=dx(i,j)+dt*((hz(i,j)-hz(i,j-1))/ddx+kz*hy(i,j));
                 ex(i,j)=-tempy_sub/tempy_add*ex(i,j)+1/para(i,j,1)/tempy_add*(tempx_add*dx(i,j)+tempx_sub*dx_d(i,j));
                 dy_d(i,j)=dy(i,j);
                 dy(i,j)=dy(i,j)-dt*((hz(i,j)-hz(i-1,j))/ddx+kz*hx(i,j));
                 ey(i,j)=-tempx_sub/tempx_add*ey(i,j)+1/para(i,j,1)/tempx_add*(tempy_add*dy(i,j)+tempy_sub*dy_d(i,j));
                 dz_d(i,j)=dz(i,j);
                 dz(i,j)=-tempy_sub/tempy_add*dz(i,j)+1/tempy_add*((hy(i,j)-hy(i-1,j))-(hx(i,j)-hx(i,j-1)))/ddx;
                 ez(i,j)=-tempx_sub/tempx_add*ez(i,j)+1/para(i,j,1)/tempx_add/dt*(dz(i,j)-dz_d(i,j));
        end
    end
    
    %保留时域波形
    
    ey20k=ey(1:I,1:J);
    ex20k=ex(1:I,1:J);
    
    %定义剩下的时间wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
    if t==1
        t_elapsed=cputime-t_start;
    end
    %wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
end % 主循环结束。
x=-178:177;y=x;
[X,Y]=meshgrid(x,y);
Z=ey20k.^2+ex20k.^2;
surf(X,Y,Z)