%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% impulse propagation through a panel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear

%%% INPUT PARAMETERS  %%%
K        = 4;                   % Oversampling Ratio
N_samples= 2048*K;              % Time Samples
v_rat    = 1;                   % Velocity Ratio (=c/ugrid)
T_max    = 1024*10^(-9);        % Time Window
dT       = T_max/N_samples;     % Time Step
dL       = 0.15/K;              % Spatial Step
N_nodes  = K*101;               % Number of E-nodes
L        = (N_nodes-1)*dL;      % Spatial Domain Length
T0       = K*16*dT;             % Truncated Gaussian Pulse Duration
T_shift  = T0/2;                % Gaussian Pulse Time Shift
B_eff    = 10^9;                % Effective Bandwidth
f_max    = 1/(2*dT);            % Nyquist Frequency
N_obs    = round(N_nodes/(K*2));% Observation Point, Total Field
c        = 3e8;
eta0     = 377;
ugrid    = dL/dT;
%%% ASSIGN INPUT VARIABLES %%%
eps_rel  = 4;               % Relative permittivity is fixed
mu_rel   = 1;               % Relative permeability is fixed
sigma    = 0;               % Conductivity inside of panel
diel_start = 6;             % Left edge of dielectric
diel_stop  = 7;             % Right edge of dielectric

%%% SET UP PLOT RANGE AND DRAW THE WALL &&&
	Lmax=15;
	Amax=0;
	Amin=-40;
	Xw1 = [diel_start, diel_start];
    Xw2 = [diel_stop, diel_stop];
    Y = [Amin, Amax];

%%% AVERAGE PERMITTIVITY AND CONDUCTIVITY (mu_rel=1 ASSUMED) %%%
B = ones(1,N_nodes-1);
C = ones(1,N_nodes);

%%% HOMOGENEOUS PANEL - A IS THE SAME FOR ALL NODES %%%
ugrid=dL/dT;

if rem(diel_start,dL) <= dL/2
   left_diel_node = ceil(diel_start/dL);
else
   left_diel_node = ceil(diel_start/dL) + 1;
end

if rem(diel_stop,dL) >= dL/2
   right_diel_node = floor(diel_stop/dL);
else
   right_diel_node = floor(diel_stop/dL) - 1;
end

for k=left_diel_node:1:right_diel_node
	B(1,k) = (1/eps_rel)/(1+v_rat*dL*sigma/2/eps_rel);
	C(1,k) = (1-v_rat*dL*sigma/2/eps_rel)/(1+v_rat*dL*sigma/2/eps_rel);
end

%%% TIME BASE %%%
time_base = linspace(0,T_max,N_samples);

%%% FREQUENCY BASE %%%
freq_base = linspace(0,(1/K)*f_max/10^9, round(N_samples/(2*K)));

%%% SPACE %%%
space_base = linspace(0,L,N_nodes);

%%% INITIAL CONDITIONS %%%
E_inc = 0;
old_E_inc = 0;
Ex = zeros(1,N_nodes);
Hy = zeros(1,N_nodes-1);
Ex_scat_0   = 0; % scattered field at z=0
Ex_scat_1   = 0; % scattered field at z=dl
Ex_scat_N   = 0; % scattered field at z=L
Ex_scat_N_1 = 0; % scattered field at z=L-1

%%% GAUSSIAN PULSE INCIDENT AT Z=0 %%%
B_eff = 10^9;

%%% TIME LOOP %%%
kk=0;
 for n=1:N_samples
  t = (n-1)*dT;

% Incident Field
  E_inc = exp((-1)*(B_eff*(t-T0/2))^2);
  E_inc_wave(n) = E_inc;

% Update Scattered Field at z=0 and z=L (grid's edges)
  Ex_scat_0 = Ex_scat_1;
  Ex_scat_N = Ex_scat_N_1;

% H-Field Update at half-time step (really = eta0*Hy)
  Hy = Hy+v_rat*(Ex(1,1:N_nodes-1)-Ex(1,2:N_nodes));

% Total Field at z=0 (adding the incident field to the grid's edge at z=0)
  Ex(1,1) = Ex_scat_0+E_inc;

% Field at z=L
  Ex(1,N_nodes) = Ex(1,N_nodes-1);

% E-Field Update at full time step, non-edge nodes
  Ex(1,2:N_nodes-1) = Ex(1,2:N_nodes-1) .* C(1,2:N_nodes-1) + ...
        v_rat*(Hy(1,1:N_nodes-2)-Hy(1,2:N_nodes-1)) .* B(1,2:N_nodes-1);

% Scattered Fields Just Inside
  Ex_scat_1  = Ex(1,2)-old_E_inc;

% Saving Old Values
  old_E_inc = E_inc;

% Observation Points Fields
  E_obs(n) = Ex(1,N_obs);
  E_obs_scat(n) = Ex_scat_0;

  if n/6==floor(n/6) & n<700
%   if n==426
    plot(space_base, 20*log10(abs(Ex)+1e-20), '-', Xw1, Y, Xw2, Y)
    axis([0, Lmax, Amin, Amax ]), grid 
    title('Gaussian Incident Field in 1-D'), xlabel('Position, m'), 
    ylabel('Amplitude, dB')
    pause(0.001)
    kk=kk+1;
    M(:,kk)=getframe;
end

end  %%% END OF TIME LOOP %%%