fdplo.m

研究FDTD模拟一维光子晶体中的光波传输

%程序目的
% 一维FDTD模拟（行波吸收边界条件）
%计算了一个平面波入射到光子晶体上的透射率频谱
%%%%%%%%%%%%%%%%%%%%%%%%

close all
clear all%清除工作空间变量

%%%%%%%%%%%%%%%%%%%%%%%%
% 程序3 个初始参数
ke=400;            % ke是对模拟区域所取的节点（计算点）的数目
T=0;                %时间置零
nsteps=4000;        % 循环的总次数  即模拟的总时间
%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%
%在边界处放置脉冲源
%此为脉冲源的3个参数
kc=100;      %光子晶体的起始点
E0=1;         %平面波的振幅
%omiga0=0.07295;    %中心频率 
%omiga0=0.08241;
%pai=3.1415;
%t0=2*pai/omiga0;
%%%%%%%%%%%%%%%%%%%%%%%%%%

%设置介质参数
kstart=kc+1;                        %介质在kc+1处开始 
epsilon1=5.5225;                          %一种介质的相对介电常数为4
epsilon2=1.9044;                          %另一种介质的相对介电常数为2
n1=sqrt(epsilon1);
n2=sqrt(epsilon2);
cb=zeros(1,ke);                 %设置参数cb，全部置零 (其含义公式1.16a)
cb(1:kstart-1)=0.5;             %真空中的cb参数为0.5
%介质中的cb参数为(0.5/epsilon)  (参见公式1.17)
a=10;                 %周期数
for l=1:a;            %循环计算介质层中的参数cb     
cb(kstart+(l-1)*11:kstart+(l-1)*11+3)=0.5/epsilon1;
cb(kstart+(l-1)*11+4:kstart+(l-1)*11+10)=0.5/epsilon2;
end
cb(kstart+a*11:kstart+a*11+3)=0.5/epsilon1;
cb(kstart+a*11+4:ke)=0.5;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%光子晶体初始化
omiga0=pi/(2*(n1*4+n2*7));
t0=2*pi/omiga0;


for m=1:200
    omiga(m)=omiga0*m/100;    %扫频
      
%%%%%%%%%%%%%%%%%%%%%%%%
%对ex和hy进行初始化，全部置0
ex=zeros(1,ke); 
hy=zeros(1,ke);
ss=zeros(1,ke);
%%%%%%%%%%%%%%%%%%%%%%%%



%为吸收边界条件准备零时参数并置零
ex_low_m1=0;
ex_low_m2=0;
ex_high_m1=0;
ex_high_m2=0;




%%%%%%%%%%%%%%%%%%%%%%%%
% 程序主体！！！！
    %%%%%%%%%%%%%%%%%%%%%%%%%%
%     % %计算部分

for n=1:nsteps   % 循环的总次数  即模拟的总时间
     T=T+1;   %记录循环的总次数（其实和n是一样的）
     
     %在low end 附近放置一个波源ex    
     if T<t0
         pulse=E0*exp(i*omiga(m)*T)*0.5*(1-cos(pi*T/t0)); 
         %pulse=E0*sin(omiga*T)*0.5*(1-cos(pi*T/t0));
     else
         pulse=E0*exp(i*omiga(m)*T);
         %pulse=E0*sin(omiga*T);
     end
         ex(6)=ex(6)+pulse; 
     
    %计算电场部分
    for k=2:ke
        ex(k)=ex(k)+cb(k)*(hy(k-1)-hy(k));  %书上公式1.16a  要理解这个公式要看书
    end
    E1=ex(6);
    E3=ex(kstart+221);
    
   
    
    
    
    %吸收边界条件
    ex(1)=ex_low_m2;
    ex_low_m2=ex_low_m1;
    ex_low_m1=ex(2);
    
    ex(ke)=ex_high_m2;
    ex_high_m2=ex_high_m1;
    ex_high_m1=ex(ke-1);
    
    
    %计算磁场部分
    for k=1:ke-1
        hy(k)=hy(k)+0.5*(ex(k)-ex(k+1));    %书上公式1.9b 要理解这个公式要看书
         ss(k)=real(ex(k)*conj(hy(k)))/2;
    end
    H1=hy(6);
    H3=hy(kstart+221);
    
%     %%%%%计算光源和透射点处的能流密度；透射率
%     SS1(T)=real(E1*conj(H1))/2;
%     SS3(T)=real(E3*conj(H3))/2;
%     TT(T)=SS3(T)/SS1(T);
%     subplot(3,1,1)
%     plot(SS1)
%     subplot(3,1,2)
%     plot(SS3)
%     subplot(3,1,3)
%     plot(TT)
    
    
      
  if m==100
      
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % 实时画出电场Ex和磁场Hy随时间的变化的情况
    figure(1)
  plot(ss,'-b')
        %X=[101 222 222 101 101];
        %Y=[-2 -2 2 2 -2];
        %fill(X,Y,'r');
        %hold on
       % plot(ss,'-b')
        
        axis([0 ke -2 2])  
        grid on
        title('平均能流密度')
         xlabel('z')
  %         ylabel('振幅')
        drawnow
  %  subplot(2,1,2)
     %    G=[101 222 222 101 101];
     %    Q=[-2 -2 2 2 -2];
      %  fill(G,Q,'r');
      %  hold on
      %   plot(hy,'b')
         
     %   axis([0 ke -2 2]) 
      %  grid on
      %  title('磁场Hy')
        % ylabel('振幅')
      %  xlabel('z')
       % drawnow
    %%%%%%%%%%%%%%%%%%%%%%%%
    pause(0.0001)
  end
end
% s1(m)=real(E1*conj(H1))/2;
s3(m)=real(E3*conj(H3))/2;
TT(m)=s3(m)/0.5;
end
hold off
figure(2)
plot(omiga/omiga0,TT);
title('光子晶体透射谱')
 xlabel('归一化频率')