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 %  保存結果