clear;clc;  
 
% 1.初始化 
 
T=200; % 迭代次數   
Nx=100;  
Ny=100;  
npml=8; %PML的網格數量  
 
c0=3*10^8;   
f=1.5*10^(9);  
lambda=c0/f;   
 
wl=10;  
dx=lambda/wl; %取lamda/10的精度  
dy=lambda/wl;  
pi=3.14159;  
 
dt=dx/(2*c0); %由于是二維情況同，取t=dx/2c0  
t0=20*dt; %為高斯脈沖的仿真作準備  
 
epsz=1/(4*pi*9*10^9);   
epsilon=1;   
sigma=0;   
 
ic=Nx/4; % 源的X位置  
jc=Ny/4; % 源的Y位置  
 
for i=1:Nx+1;  
    for j=1:Ny+1;  
        dz(i,j)=0; %z方向電位移  
        ez(i,j)=0; % z方向電場  
        hx(i,j)=0; % x方向磁場  
        hy(i,j)=0; % y方向磁場   
        ihx(i,j)=0;  
        ihy(i,j)=0;  
        iz(i,j)=0;  % z方向求和中間參量  
    end;  
end;  
 
for i=2:Nx;   
    for j=2:Ny;  
        ga(i,j)=1;  
    end;  
end;   
%PML參數的設置 
for i=1:Nx;   
    gi2(i)=1;  
    gi3(i)=1;  
    fi1(i)=0;  
    fi2(i)=1;  
    fi3(i)=1;  
end  
 
for j=1:Ny;  
    gj2(j)=1;  
    gj3(j)=1;  
    fj1(j)=0;  
    fj2(j)=1;  
    fj3(j)=1;  
end  
 
for i=1:npml+1; %設置PML層中的參數  
    xnum=npml+1-i;  
    xn=0.33*(xnum/npml)^3; %這里的0.333并不是嚴格的計算，而是經驗值  
    gi2(i)=1.0/(1+xn);  
    gi2(Nx-1-i)=1/(1+xn);  
    gi3(i)=(1-xn)/(1+xn);  
    gi3(Nx-1-i)=(1-xn)/(1+xn);  
 
    xn=0.25*((xnum-0.5)/npml)^3; %這里的0.25和0.333也是一個道理  
    fi1(i)=xn;  
    fi1(Nx-2-i)=xn;  
    fi2(i)=1.0/(1+xn);  
    fi2(Nx-2-i)=1/(1+xn);  
    fi3(i)=(1-xn)/(1+xn);  
    fi3(Nx-2-i)=(1-xn)/(1+xn);  
end  
 
for i=1:npml+1;  
    xnum=npml+1-i;  
    xn=0.33*(xnum/npml)^3;  
    gj2(i)=1.0/(1+xn);  
    gj2(Ny-1-i)=1/(1+xn);  
    gj3(i)=(1-xn)/(1+xn);  
    gj3(Ny-1-i)=(1-xn)/(1+xn);  
 
    xn=0.25*((xnum-0.5)/npml)^3;  
    fj1(i)=xn;  
    fj1(Ny-2-i)=xn;  
 
    fj2(i)=1.0/(1+xn);  
    fj2(Ny-2-i)=1/(1+xn);  
    fj3(i)=(1-xn)/(1+xn);  
    fj3(Ny-2-i)=(1-xn)/(1+xn);  
end  
 
%2.迭代求解電場和磁場  
 
for t=1:T;  
    for i=2:Nx;  % 為了使每個電場周圍都有磁場進行數組下標處理  
        for j=2:Ny;  
      dz(i,j)=gi3(i)*gj3(j)*dz(i,j)+gi2(i)*gj2(j)*0.5*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1)); % 
        end;  
    end;  % 電場循環結束  
     
    pulse=sin(2*pi*f*t*dt);  % 正弦波源  
  %pulse=exp(-4*pi*(t*dt-t0)^2/(2/f)^2);高斯脈沖情況  
    dz(ic,jc)=dz(ic,jc)+pulse;  %軟源，就是為了防止反射做的一個處理  
     
    for i=1:Nx; %為了使每個電場周圍都有磁場進行數組下標處理  
        for j=1:Ny;  
          ez(i,j)=ga(i,j)* dz(i,j);  %反映煤質的情況都是放到這里的  
        end;  
    end; % 電荷密度循環結束  
     
    for j=1:Ny;  
        ez(1,j)=0;  
        ez(Nx,j)=0;  
    end  
 
    for i=1:Nx;  
        ez(i,1)=0;  
        ez(i,Ny)=0;  
    end  
 
    for i=1:Nx;  % 為使每個磁場周圍都有電場進行數組下標處理  
        for j=1:Ny-1;  
            curl_e=ez(i,j)-ez(i,j+1);  
          ihx(i,j)=ihx(i,j)+fi1(i)*curl_e;  
          hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5*(curl_e+ihx(i,j));  
        end;  
    end; % 磁場HX循環結束  
 
    for i=1:Nx-1;  % 為了使每個磁場周圍都有電場進行數組下標處理  
        for j=1:Ny;  
            curl_e=ez(i+1,j)-ez(i,j);  
          ihy(i,j)=ihy(i,j)+fj1(j)*curl_e;  
          hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5*(curl_e+ihy(i,j));  
        end;  
    end;  % 磁場HY循環結束  
    mesh(ez);  
    drawnow;  
end;