incTime=readmatrix('inc-8192.dat'); % incident field file                                     
dieTime=readmatrix('die-8192.dat'); % transmitted field file                                  

inc=fft(incTime); % take Fourier transforms                                                   
die=fft(dieTime);

nSteps=length(incTime); % number of time steps                                                
freqMin=1;   % minimum frequency index of interest                                            
freqMax=500; % maximum frequency of interest                                                  
freqIndex=freqMin:freqMax; % range of frequencies of interest                                 
freqSlice=freqIndex+1; % correct for offset of 1 in matlab's indexing                         

courantNumber = 1;
nLambda=nSteps./freqIndex; % points per wavelength of freq's of interest                      

clf
subplot(3,1,1)
hold on
plot(incTime(freqSlice),'-.');
plot(dieTime(freqSlice));
grid on
legend('Incident field','Transmitted field', ...
   'Interpreter', 'latex');
xlabel('Temporal Index', 'Interpreter', 'latex');
ylabel('$E_z(x_1, t)$ [V/m]', ...
   'Interpreter', 'latex');
title('(a) Temporal Fields', 'Interpreter', 'latex');

subplot(3,1,2)
hold on
plot(freqIndex,abs(inc(freqSlice)),'-.');
plot(freqIndex,abs(die(freqSlice)));
grid on
legend('Incident field','Transmitted field', ...
   'Interpreter', 'latex');
xlabel('Frequency Index', 'Interpreter', 'latex');
ylabel('$|E_z(x_1,\omega)|$', 'Interpreter', 'latex');
title('(b) Transform vs. Frequency Index', ...
   'Interpreter', 'latex');

subplot(3,1,3)
plot(freqIndex,abs(die(freqSlice)./inc(freqSlice)));
grid on
xlabel('Frequency Index', 'Interpreter', 'latex');
ylabel('$|E^t_z(x_1,\omega)/E^i_z(x_1,\omega)|$', ...
   'Interpreter', 'latex');
title('(c) Normalized Transmitted Field vs. Frequency Index', ...
   'Interpreter', 'latex');
