请教一个用MATLAB编的一维FDTD程序（有程序）

要仿真一个一维的FDTD传播问题，想根据这个改动一下，在距离高斯点源15cm的地方加入一个相对介电常数为9，厚度为5CM的介质，边界条件由PEC改为吸收边界，请教大虾给条明路.... [txE .7p  
% Define initial constants K;?+8(H  
eps_0 = 8.854187817e-12;          % permittivity of free space  e'~3oqSvR  
mu_0  = 4*pi*1e-7;                % permeability of free space    zhQJy?>'m  
c     = 1/sqrt(mu_0*eps_0);       % speed of light I7onX,U+  
="+#W6bZT  
% Define problem geometry and parameters A.SvA Yn  
domain_size = 0.6;                  % 1D problem space length in meters ?,z}%p  
dx = 1e-3;                        % cell size in meters   j2k"cmsKh  
dt = 3e-12;                       % duration of time step in seconds  y29m/i:  
number_of_time_steps = 2000;      % number of iterations { 6il`>=C  
nx = round(domain_size/dx);       % number of cells in 1D problem space kiEa<-]  
source_position = 0.2;            % position of the current source Jz w)f#V s  
2y4bwi  
% Initialize field and material arrays 
$'vU2L  
Ceze      = zeros(nx+1,1); Gc?a+T  
Cezhy     = zeros(nx+1,1); i-1op> Y  
Cezj      = zeros(nx+1,1); `5*}p#G  
Ez        = zeros(nx+1,1); %{W6PrY{  
Jz        = zeros(nx+1,1); 1MFbQs^  
eps_r_z   = ones (nx+1,1); % free space 00
(\ZUj  
sigma_e_z = zeros(nx+1,1); % free space /ZX}Nc g  
\bXa&Lq  
Chyh      = zeros(nx,1); \fOEqe*5SM  
Chyez     = zeros(nx,1); pa+hL,w{6  
Chym      = zeros(nx,1); :
O
T&  
Hy        = zeros(nx,1); pglVR </  
My        = zeros(nx,1); E.h*g8bXe  
mu_r_y    = ones (nx,1); % free space W/N7vAx X  
sigma_m_y = zeros(nx,1); % free space z{q`GwW  
U{mYTN*:j$  
% Calculate FDTD updating coefficients KI.unP%  
Ceze = (2 * eps_r_z * eps_0 - dt * sigma_e_z) ... }e1ZbmW  
     ./(2 * eps_r_z * eps_0 + dt * sigma_e_z); w0. u\  
+{]j]OP  
Cezhy = (2 * dt / dx) ... iZmcI;?u  
     ./(2 * eps_r_z * eps_0 + dt * sigma_e_z); >P(.:_^p  
kh<2BOV  
Cezj  = (-2 * dt) ... ?,/ }`3Vw  
     ./(2 * eps_r_z * eps_0 + dt * sigma_e_z); h[ ZN+M  
i8p6Xht  
Chyh  = (2 * mu_r_y * mu_0 - dt * sigma_m_y) ... Wwo0%<2y  
     ./(2 * mu_r_y * mu_0 + dt * sigma_m_y); X}]-*T|a  
+`4A$#$+y  
Chyez = (2 * dt / dx) ... sOY:e/_F  
     ./(2 * mu_r_y * mu_0 + dt * sigma_m_y); l/D} X  
;uW FHc5@B  
Chym  = (-2 * dt) ... ?dTD\)%A  
     ./(2 * mu_r_y * mu_0 + dt * sigma_m_y); Z+SRXKQ  
\U0Q<ot/7  
% Define the Gaussian source waveform :RYTL'hes  
time       = dt*[0:number_of_time_steps-1].'; x`s>*^  
Jz_waveform = exp(-((time-8e-11)/2e-11).^2); , gHDx  
source_position_index = round(nx*source_position/domain_size)+1; _1^'(5f$  
[g,}gyeS(  
\V:^h[ad  
Ez_positions = [0:nx]*dx; *8q.YuZ  
Hy_positions = ([0:nx-1]+0.5)*dx; +ZYn? #IQ  
v = [0 -0.1 -0.1; 0 -0.1 0.1; 0 0.1 0.1; 0 0.1 -0.1; ... XppOU  
     1 -0.1 -0.1; 1 -0.1 0.1; 1 0.1 0.1; 1 0.1 -0.1]; ZCw]m#lS  
f = [1 2 3 4; 5 6 7 8]; =4!mAo}  
axis([0 1 -0.2 0.2 -0.2 0.2]); f$( e\++  
lez = line(Ez_positions,Ez*0,Ez,'Color','b','LineWidth',1.5); gw(z1L5 n  
lhy = line(Hy_positions,377*Hy,Hy*0,'Color','r', ... AA_%<zK  
    'LineWidth',1.5,'linestyle','-.'); {g6%(X\r.r  
set(gca,'fontsize',12,'FontWeight','bold'); Cx"sw }  
axis square; al0L&z\  
legend('E_{z}', 'H_{y} \times 377','Location','NorthEast'); M|-)GvR$J  
xlabel('x [m]'); }Z>)DN=+  
ylabel('[A/m]'); A&{Nh` q
 
zlabel('[V/m]'); ~&O%N  
grid on; zs;JJk^  
p = patch('vertices', v, 'faces', f, 'facecolor', 'g', 'facealpha',0.2); pAEx#ck  
~[: 2I  

Dq xs+
 
% FDTD loop 1YA% -~  
for time_step = 1:number_of_time_steps
@HW*09TG  
Ac6=(B  
    if (0.35=<Ez_positions<=0.4)%%%%%这样改肯定不对...额 X&zis1A<  
eps_r_z =9.*ones (nx+1,1); Vl]>u+YqE  
else }u|q0>^8  
eps_r_z=ones (nx+1,1); 'qi}|I  
end 9uY'E'm*  
    % Update Jz for the current time step :gT4K-Oj  
    Jz(source_position_index) = Jz_waveform(time_step); E7hhew  
    E^PB)D(.  
    % Update magnetic field )jj0^f1!j  
        Hy(1:nx) =   Chyh(1:nx) .* Hy(1:nx) ... a.'*G6~Qgw  
         + Chyez(1:nx) .* (Ez(2:nx+1) - Ez(1:nx))  ... J4utIGF  
         + Chym(1:nx) .* My(1:nx); |qLh5Ty  
0x7'^Z>-oe  
    % Update electric field }G=M2V<L  
    Ez(2:nx) = Ceze (2:nx) .*  Ez(2:nx) ... N!32 wJ  
             + Cezhy(2:nx) .* (Hy(2:nx) - Hy(1:nx-1))  ... /?!u{(h}  
             + Cezj(2:nx)  .*  Jz(2:nx); l6B@qYLZ  
R]dg_Da  
    Ez(1)    .. )"LJ hLg  
:,^gj  

未注册仅能浏览部分内容，查看全部内容及附件请先 登录 或 注册

第78-82是源文件没有的，要删除 6,uX,X5