1121204716163470.m

模拟光在二维介质中传输

function FDTD_debug

%constants
c_0 = 3E8;                   % Speed of light in free space
Nz = 100;                    % Number of cells in z-direction
Nt = 200;                    % Number of time steps
N = Nz;                      % Global number of cells

E = zeros(Nz,1);             % E component
E2 = zeros(Nz,1);
E1 = zeros(Nz,1);
Es = zeros(Nz,Nz);           % Es is the output for surf function
%H = zeros(Nz,1);             % H component
pulse = 0;         % Pulse

%to be set
mu_0 = 4.0*pi*1.0e-7;      % Permeability of free space
eps_0 = 1.0/(c_0*c_0*mu_0);  % Permittivty of free space

c_ref = c_0;                 % Reference velocity
eps_ref = eps_0;
mu_ref = mu_0;

f_0 = 10e9;                  % Excitation frequency
f_ref = f_0;                 % Reference frequency 

omega_0 = 2.0*pi*f_0;        % Excitation circular frequency
omega_ref = omega_0;

lambda_ref = c_ref/f_ref;    % Reference wavelength

Dx_ref = lambda_ref/20;      % Reference cells width
Dz = Dx_ref;
nDz = Dz/Dx_ref;

Z = Nz*Dz;
r = 1;                       % Normalized time step
Dt_ref = r*Dx_ref/c_ref;     % Reference time step
Dt = Dt_ref;

% Source position
Nz_Source = 0.5*Nz;
N_Source = Nz_Source;

for n = 1:Nt-1
    t = Dt_ref*r*(n-1);     % Actual time
    
    %Source function series
    Source_type = 1;
    switch Source_type 
        case 1  % modified source function
            ncycles = 2;
            if t < ncycles*2.0*pi/(omega_0)  
                pulse = -0.5*( 1.0 - cos(omega_0*t/ncycles) ) * sin(omega_0*t);
            else 
                pulse = 0;
            end
        case 2  % sigle cos source function
            if t < 2.0*pi/(omega_0)
                pulse = 8*c_0^2*Dt_ref^2*mu_ref*omega_0*cos(omega_0*t);
            else 
                pulse = 0;
            end
        case 3  % Gaussian pulse
            if t < Dt_ref*r*50
                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);
            else
                pulse = 0;
            end      
    end
    E2(N_Source) = E2(N_Source) - r*pulse;
    E1 = E2;
    E2 = E;
       
    m = 2:Nz-1;
    E(m) = r^2*(E2(m+1) - 2*E2(m) + E2(m-1)) + 2*E2(m) - E1(m);
        
    %Boundary
    E(1) = 0; E(Nz) = 0;                                % Dirichlet
    %E(1) = E(2);E(Nz) = E(Nz-1);                       % Neumann
    %E(Nz) = E2(Nz-1) + (r-1)/(r+1)*(E(Nz-1)-E2(Nz));    % Mur ABC @ right side z = Nz
    %E(1) = E2(2) + (r-1)/(r+1)*(E(2)-E2(1));            % Mur ABC @ left side z = 0
       
    %display
    %choice***********************************************
    display = 3; % choice: '1' = line plot; '2' = imagesc; '3' = surf
    switch display
        case 1
            i = 1:Nz;
            plot(i, E(i));
            axis([0 Nz -4 4]);
        case 2
            A = rot90(E);
            imagesc(A);
        case 3
            i = 1:Nz;
            for j = 1:Nz
                Es(i,j) = E(i);
            end
            Es = rot90(Es);
            j = 1:Nz;
            surf(i,j,Es);
            axis([0 Nz 0 Nz -4 4])    
        otherwise
            A = rot90(E);
            imagesc(A);
    end 
    pause(0.005);  
end