《FDTD electromagnetic field using MATLAB》讀書筆記之 Figure 1.14

標籤：sig   tla   field   程式   sof   seconds   電磁波   產生   fonts   

       背景：

        基於公式1.42（Ez分量）、1.43（Hy分量）的1D FDTD實現。

        計算電場和磁場分量，該分量由z方向的電流片Jz產生，Jz位於兩個理想導體極板中間，兩個極板平行且向y和z方向無限延伸。

          平行極板相距1m，差分網格Δx=1mm。

        電流面密度導致分介面（電流薄層）磁場分量的不連續，在兩側產生Hy的波，每個強度為5×10^（-4）A/m。因為電磁波在自由空間中傳播，本徵阻抗為η0。

          顯示了從左右極板反射前後的傳播過程。本例中是PEC邊界，正切電場分量（Ez）在PEC表面消失。

        觀察圖中step650到700的場的變換情況，經過PEC板的入射波和反射波的傳播特徵。經過PEC反射後，Ez的極性變反，這是因為反射係數等於-1；

而磁場分量Hy沒有反向，反射係數等於1。

下面是書中的代碼（幾乎沒動）：

       第1個是主程式：

%% ------------------------------------------------------------------------------%%            Output Info about this m-filefprintf(‘\n****************************************************************\n‘);fprintf(‘\n   <FDTD 4 ElectroMagnetics with MATLAB Simulations>     \n‘);fprintf(‘\n                                             Listing A.1     \n\n‘);time_stamp = datestr(now, 31);[wkd1, wkd2] = weekday(today, ‘long‘);fprintf(‘           Now is %20s, and it is %7s  \n\n‘, time_stamp, wkd2);%% ------------------------------------------------------------------------------% This program demonstrates a one-dimensional FDTD simulation.% The problem geometry is composed of two PEC plates extending to% infinity in y and z dimensions, parallel to each other with 1 meter% separation. The space between the PEC plates is filled with air.% A sheet of current source paralle to the PEC plates is placed% at the center of the problem space. The current source excites fields% in the problem space due to a z-directed current density Jz,% which has a Gaussian waveform in time.% Define initial constantseps_0 = 8.854187817e-12;                  % permittivity of free spacemu_0  = 4*pi*1e-7;                        % permeability of free spacec     = 1/sqrt(mu_0*eps_0);               % speed of light% Define problem geometry and parametersdomain_size = 1;                          % 1D problem space length in metersdx = 1e-3;                                % cell size in meters, Δx=0.001mdt = 3e-12;                               % duration of time step in secondsnumber_of_time_steps = 2000;              % number of iterationsnx = round(domain_size/dx);               % number of cells in 1D problem spacesource_position = 0.5;                    % position of the current source Jz% Initialize field and material arraysCeze       = zeros(nx+1, 1);Cezhy      = zeros(nx+1, 1);Cezj       = zeros(nx+1, 1);Ez         = zeros(nx+1, 1);Jz         = zeros(nx+1, 1);eps_r_z    = ones (nx+1, 1);              % free spacesigma_e_z  = zeros(nx+1, 1);              % free spaceChyh       = zeros(nx, 1);Chyez      = zeros(nx, 1);Chym       = zeros(nx, 1);Hy         = zeros(nx, 1);My         = zeros(nx, 1);mu_r_y     = ones (nx, 1);                % free spacesigma_m_y  = zeros(nx, 1);                % free space% Calculate FDTD updating coefficientsCeze       = (2 * eps_r_z * eps_0 - dt * sigma_e_z) ...           ./(2 * eps_r_z * eps_0 + dt * sigma_e_z); Cezhy      = (2 * dt / dx) ...           ./(2 * eps_r_z * eps_0 + dt * sigma_e_z);Cezj       = (-2 * dt) ...           ./(2 * eps_r_z * eps_0 + dt * sigma_e_z);Chyh       = (2 * mu_r_y * mu_0 - dt * sigma_m_y) ...           ./(2 * mu_r_y * mu_0 + dt * sigma_m_y);Chyez      = (2 * dt / dx) ...           ./(2 * mu_r_y * mu_0 + dt * sigma_m_y);Chym       = (-2 *dt) ...           ./(2 * mu_r_y * mu_0 + dt * sigma_m_y);% Define the Gaussian source waveformtime                  = dt * [0:number_of_time_steps-1].‘;Jz_waveform           = exp(-((time-2e-10)/5e-11).^2)*1e-3/dx;source_position_index = round(nx * source_position/domain_size)+1; % Subroutine to initialize plottinginitialize_plotting_parameters;% FDTD loopfor time_step = 1:number_of_time_steps% Update Jz for the current time stepJz(source_position_index) = Jz_waveform(time_step);% Update magnetic fieldHy(1:nx) = Chyh(1:nx) .* Hy(1:nx) ...     + Chyez(1:nx) .* (Ez(2:nx+1) - Ez(1:nx)) ...     + Chym(1:nx) .* My(1:nx);% Update electric fieldEz(2:nx) = Ceze (2:nx) .* Ez(2:nx) ...         + Cezhy(2:nx) .* (Hy(2:nx) - Hy(1:nx-1)) ...         + Cezj (2:nx) .* Jz(2:nx);         Ez(1)    = 0;       % Apply PEC boundary condition at x = 0 mEz(nx+1) = 0;       % Apply PEC boundary condition at x = 1 m% Subroutine to plot the current state of the fieldsplot_fields;end

　　第2個是initialize_plotting_parameters，看名字就知道是初始化參數：

% Subroutine used to initialize 1D plot Ez_positions = [0:nx]*dx;Hy_positions = ([0:nx-1]+0.5)*dx;v = [0 -0.1 -0.1; 0 -0.1 0.1; 0 0.1 0.1; 0 0.1 -0.1; ...     1 -0.1 -0.1; 1 -0.1 0.1; 1 0.1 0.1; 1 0.1 -0.1];f = [1 2 3 4; 5 6 7 8];axis([0 1 -0.2 0.2 -0.2 0.2]);lez = line(Ez_positions, Ez*0, Ez, ‘Color‘, ‘b‘, ‘linewidth‘, 1.5);lhy = line(Hy_positions, 377*Hy, Hy*0, ‘Color‘, ‘r‘, ‘LineWidth‘, 1.5, ‘LineStyle‘,‘-.‘);set(gca, ‘fontsize‘, 12, ‘fontweight‘, ‘bold‘);set(gcf,‘Color‘,‘white‘); axis square;legend(‘E_{z}‘, ‘H_{y} \times 377‘, ‘location‘, ‘northeast‘);xlabel(‘x [m]‘);ylabel(‘[A/m]‘);zlabel(‘[V/m]‘);grid on;p = patch(‘vertices‘, v, ‘faces‘, f, ‘facecolor‘, ‘g‘, ‘facealpha‘, 0.2);text(0, 1, 1.1, ‘PEC‘, ‘horizontalalignment‘, ‘center‘, ‘fontweight‘, ‘bold‘);text(1, 1, 1.1, ‘PEC‘, ‘horizontalalignment‘, ‘center‘, ‘fontweight‘, ‘bold‘);

　　第3個就是畫圖：

% Subroutine used to plot 1D transient field delete(lez);delete(lhy);lez = line(Ez_positions, Ez*0, Ez, ‘Color‘, ‘b‘, ‘LineWidth‘, 1.5);lhy = line(Hy_positions, 377*Hy, Hy*0, ‘Color‘, ‘r‘, ‘LineWidth‘, 1.5, ‘LineStyle‘, ‘-.‘);ts  = num2str(time_step);ti  = num2str(dt*time_step*1e9);title([‘time step = ‘ ts ‘ , time = ‘ ti  ‘ ns‘]);drawnow;

　　運行結果：

 

        ，從PEC板反射後，電場分量極性變反，磁場分量極性不變。

        反射後，電場分量極性再次改變；

 

《FDTD electromagnetic field using MATLAB》讀書筆記之 Figure 1.14