lossy_fdtd.m

ABC_FDTD_Die(T) Implements simulation of a Gaussian Pulse over T time steps. ABC are for free spac

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lossy_FDTD FDTD in a Lossy Medium
% Lossy_FDTD(T,F) Implements a sinusiodal pulse with frequency F
% through the% problem space for T time steps.
% For a non-Lossy medium set sigma = 0.
% For Free space do the above and set dielectric_constant =1.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function Lossy_FDTD = Lossy_FDTD(time,freq)

delta_x = .01;
delta_t = delta_x/(2*3e8);
max_time = time;
max_space = 200;            
E = zeros(max_space,1);            
H = E;     
dielectric_begin = 100;
dielectric_end = max_space;
dielectic_space_a = ones(max_space,1);
dielectic_space_b = ones(max_space,1);

eta = 0.5;
dielectric_constant = 4;
etaeps = eta/dielectric_constant;

%Lossy Medium variables - for non-Lossy medium set sigma = 0;
sigma = 0.04;
epszero = 8.85419e-12;
loss_term = (delta_t*sigma) / (2*epszero*dielectric_constant);

                 
%ABS needed for LHS of problem space do not work with Lossy medium
E_low_m2 = 0;
E_low_m1 = 0;
E_high_m2 = 0;
E_high_m1 = 0;

%Initialize Pulse Variables
frequency = freq*1e6;
place_pulse = 5;
spread = 12;
center_problem_space = max_space/2;
t0 = 40.0; 

%Set up dielectric medium withing the problem space
for i = 1:max_space
    if (i>=dielectric_begin & i <= dielectric_end)
        dielectric_space_a(i)=(1-loss_term)/(1+loss_term);
        dielectric_space_b(i)=0.5/(dielectric_constant *(1+loss_term));
    else
        dielectric_space_a(i)=1;
        dielectric_space_b(i)=eta;
    end
end
       

for n = 1:max_time

    for k = 2:max_space
        E(k) =  dielectric_space_a(k)*E(k) + dielectric_space_b(k)*(H(k-1)-H(k));   
    end
        
    pulse = sin(2*pi*frequency*delta_t*n);
    E(place_pulse) = E(place_pulse) + pulse;
    
    E(1) = E_low_m2;
    E_low_m2 = E_low_m1;
    E_low_m1 = E(2);
    E(max_space) = E_high_m2;
    E_high_m2 = E_high_m1;
    E_high_m1 = E(max_space-1);
    
    for j = 1:max_space-1
        H(j) = H(j) +eta*(E(j)-E(j+1));    
    end
   
       plot(E)
       title('Lossy Medium Simulation')
       ylabel('E_x')
       xlabel('Problem Space')
       axis([0 max_space -1.2 1.2])
       drawnow
end