try22.m

二维光子晶体弯曲波导的fdtd程序,matlab语言 ,TE波, PML边界条件,平面波光源

clear %  清空内存
clc
opengl neverselect %  系统显示设置

nm=1e-9;
um=1e-6;  %  单位设定 

L=600*nm;  %  入射波长

global I J c dx dy dt %  I,J为网格数，c是真空光速，dx、dy，dt分别是空间和时间步长

e0=(1e-9)/(36*pi);    %  真空的介电常数
u0=pi*4e-7;           %  真空的磁导率 
g0=0;                 %  电导率  

c=1/sqrt(e0*u0);      %  真空光速
w=2*pi*c/L;           %  圆频率
K=2*pi/L;             %  波矢

% 计算空间参数
e1=2.3^2*e0;
u1=u0;
g1=0;          %  玻璃介质参数

e2=e0;
u2=u0;
g2=getcon(L);  %  银膜介质参数

ds=20*nm;
dx=ds;
dy=ds;
dt=ds/(2*c);   %  设置计算步长

I=200;
J=100;%  设置计算空间的尺寸

Is=100;
Js=50; %  设置波源位置

%  电磁场分量定义(TM波)
Hz=zeros(I,J,3);
Ex=zeros(I,J-1,3);
Ey=zeros(I-1,J,3);

%-----------------------------------------
%  设置空间参数

setcon %  设置材料参数子程序
% break
%  M=Model(imread('crystal.bmp')); %  计算模型保存为bmp格式的位图(黑白)文件
%                                  %  黑色代表介质，白色代表空气
%                                  %  Model函数将位图转化为数值矩阵 
%  C(:,:)=M*C1+~M*C0;
%  D(:,:)=M*D1+~M*D0;
 
% 除用作图法得到模型形状外，此处另可添加设置模型的其它子程序

%设置光源
 Sd=length(Js);   %  波源长度
 SL=1:length(Js);
 SW=20;           %  高斯光源半高宽
 S_center=(1+J)/2;%  波源中心位置 
%-----------------------------------------
for n=0:1000
    
if n<pi/w/dt
   U=.5*(1-cos(pi*n/(pi/w/dt)));
else 
    U=1;
end   %设置升余弦开关函数,减少波动。

 S=exp(-j*w*n*dt)*U;  %  高斯波源

 movon %  FDTD计算子程序
  

if mod(n,4)==0
    %  每隔20步显示计算结果
    
fprintf('第 %d 时间步已完成\n',n);

% EI=abs((Ex(1:I-1,1:J-1,3))).^2+abs((Ey(1:I-1,1:J-1,3))).^2; %  光强分布
EI=real(Hz(1:I-1,1:J-1));
% look  %  数值结果可视化子程序
 pcolor(EI')
 shading interp
 axis image
 caxis([-.1,.1])
disp('正在计算......');
end
pause(0.001);
end
% save result Ez Hx Hy EI I J n dx dy ds %  保存结果