[讨论]高手来看下这个一维的FDTD的MATLAB仿真

function FDTD_debug 77bZ  
;hRpAN  
%constants />j+7ts  
c_0 = 3E8;                  % Speed of light in free space \kJt@ [w%  
Nz = 100;                    % Number of cells in z-direction 0f}Q~d=QL  
Nt = 200;                    % Number of time steps i!+3uHWu`)  
N = Nz;                      % Global number of cells VA&OI;=ri  
5Z>pa`_$2  
E = zeros(Nz,1);            % E component _F$t#.o  
E2 = zeros(Nz,1); 3x;y}:wQa  
E1 = zeros(Nz,1); m^u&g&^  
Es = zeros(Nz,Nz);          % Es is the output for surf function `]
 dx%  
%H = zeros(Nz,1);            % H component F8r455_W"  
pulse = 0;        % Pulse Z=5}17kA  
%dWFg<< |  
%to be set QIz N
#;g  
mu_0 = 4.0*pi*1.0e-7;      % Permeability of free space CR8r|+(8  
eps_0 = 1.0/(c_0*c_0*mu_0);  % Permittivty of free space >*Z{@1*h  
QT&Ws+@ s{  
c_ref = c_0;                % Reference velocity FSZoT!  
eps_ref = eps_0; t-gNG!B  
mu_ref = mu_0; hm} :Me$[)  
EvardUB)  
f_0 = 10e9;                  % Excitation frequency 2BU)qv-  
f_ref = f_0;                % Reference frequency $(mdz)Cfy  
]TZWFL-  
omega_0 = 2.0*pi*f_0;        % Excitation circular frequency GWE0 UO}  
omega_ref = omega_0; {O,M
}0Eg  
~FrkLP  
lambda_ref = c_ref/f_ref;    % Reference wavelength (,9cCnvmYU  
{{)[Ap)  
Dx_ref = lambda_ref/20;      % Reference cells width ii]=C(e9  
Dz = Dx_ref; #WmAkzvq  
nDz = Dz/Dx_ref; 3usA  
5Yibv6:3a  
Z = Nz*Dz; YH+\rb_  
r = 1;                      % Normalized time step "Ohpb!J9  
Dt_ref = r*Dx_ref/c_ref;    % Reference time step [q{Txe  
Dt = Dt_ref; #1hz=~YO  
pj-HLuZR  
% Source position R~c vml  
Nz_Source = 0.5*Nz; H5MAN,`  
N_Source = Nz_Source; TOF62,
 
<XcMc<h~  
for n = 1:Nt-1 = (h;L$  
    t = Dt_ref*r*(n-1);    % Actual time b0x0CMf  
    ^9f`3~!#bc  
    %Source function series vrO$8* sy  
    Source_type = 1; V+<AG*[  
    switch Source_type U.sPFt  
        case 1  % modified source function m~IWazj;A  
            ncycles = 2; ?g#t3j>zoF  
            if t < ncycles*2.0*pi/(omega_0)  GyM%vGl 3  
                pulse = -0.5*( 1.0 - cos(omega_0*t/ncycles) ) * sin(omega_0*t); v2X0Px_  
            else v^I%Wm  
                pulse = 0; {gHscj;SM  
            end adPU)k_j:  
        case 2  % sigle cos source function \kGtYkctZ  
            if t < 2.0*pi/(omega_0) 2_~XjwKE  
                pulse = 8*c_0^2*Dt_ref^2*mu_ref*omega_0*cos(omega_0*t); U;dt-3?=.h  
            else csA.3|rv  
                pulse = 0; E/-Kd!|"  
            end ph (k2cb  
        case 3  % Gaussian pulse |a:Vp
M  
            if t < Dt_ref*r*50 ;Sl0kSu  
                pulse = -40*c_0*(t-t*25/(n-1))*exp(-(t-t*25/(n-1))^2/2/(50/2.3548)^2/(t/(n-1))^2); g<tTZD\g  
            else "X}F%:HL  
                pulse = 0; fU7:3"|s8  
            end      .wOLi Ms  
    end bc}OmPE  
    E2(N_Source) = E2(N_Source) - r*pulse; ~g6[ [  
    E1 = E2; naCI55Wx  
    E2 = E; {%~Ec4r  
      f]65iE?x  
    m = 2:Nz-1; &fhurzzAm  
    E(m) = r^2*(E2(m+1) - 2*E2(m) + E2(m-1)) + 2*E2(m) - E1(m); x
9 L\"  
        ;m:GUp^[  
    %Boundary n\al}KG  
    E(1) = 0; E(Nz) = 0;                                % Dirichlet .wn_e=lT  
    %E(1) = E(2);E(Nz) = E(Nz-1);                      % Neumann _llaH  
    %E(Nz) = E2(Nz-1) + (r-1)/(r+1)*(E(Nz-1)-E2(Nz));    % Mur ABC @ right side z = Nz _q}%!#4  
    %E(1) = E2(2) + (r-1)/(r+1)*(E(2)-E2(1));            % Mur ABC @ left side z = 0 ;4#8#
;  
      {+mkXp])R  
    %display (G!J=
=  
    %choice*********************************************** b'"%  
    display = 3; % choice: '1' = line plot; '2' = imagesc; '3' = surf MQMy Z:  
    switch display BQ)43Rr>  
        case 1 CKw)J}z  
            i = 1:Nz; }ucg!i3C  
            plot(i, E(i)); zk~rKQ,  
            axis([0 Nz -4 4]); BCB/cBE  
        case 2 aT1W]i  
            A = rot90(E); CPE F,,\  
            imagesc(A); 3t6'5{  
        case 3 `l#$l3v+  
            i = 1:Nz; fB}5,22  
            for j = 1:Nz r>@/XYK&\  
                Es(i,j) = E(i); 69[k ?')LM  
  &nbs .. 7%}}m&A7h  
{!bJ.O l  

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

一维FDTD的最多数据是一维的，还真不能绘制3D的图。 H;ib3
?  
不过可以同时绘制电场、磁场分量的历史曲线，他们在垂直的两个平面，这样可以是3D的

不像1D的FDTD啊，像一维时域波动方程