%//************************************************************************
%//    fig_3.m: MATLAB visualization program
%//***********************************************************************
%//    Program author: Dr.-Ing. René Marklein
%//                    University of Kassel, Department of Electrical and 
%//                    Computer Engineering, Electromagnetic Theory
%//                    Wilhelmshöher Allee 71
%//                    D-34121 Kassel
%//                    Tel. 0561 804 6426, marklein@uni-kassel.de
%//
%//    Date of this version:  November 15, 2001
%//
%//    This program reads the binary data of the time history of the excitation
%//    pulsem plots the time history, computes the temporal frequency spectrum,
%//    and plots the magnitude of the temporal frequency spectrum.
%//************************************************************************

clear all;
close all;

fprintf('\nVisualize the time history and its temporal frequency spectrum of the\n');
fprintf('exciation function read from binary data file.\n');

Nt = 250;
fprintf('Number of time steps: Nt = %d\n', Nt);

Dt = 2.5e-12;
fprintf('Sampling time rate  : Dt = %g s\n', Dt);
fprintf('Sampling frequency  :  f = %g GHz\n', 1/Dt/1e9);

f=zeros(Nt);
e=zeros(Nt);

fp=fopen('Pulse.bin','rb');
[e,nread] = fread(fp,Nt,'float32');
fprintf('Number of time samples read: Nt = %d\n', nread)

% Undersampling in the time-domain
undersampling_rate = 4;
fprintf('Perform an undersampling in time-domain: factor = %d\n', undersampling_rate);
j=1;
for i=1:undersampling_rate:Nt
    f(j) = e(i);
    j=j+1;
end

Dt=undersampling_rate*Dt;
T= Dt*Nt;

omega_min = -pi/Dt;
Domega    = -2.0*omega_min/Nt;
omega_max = -omega_min+Domega;

fprintf('omega_min = %g GHz\n', omega_min/1e9);
fprintf('omega_max = %g GHz\n', omega_max/1e9);
fprintf('Domega    = %g MHz\n', Domega/1e6);

f_min = -1/Dt/2;
Df    = -omega_min/Nt/pi;
f_max = -f_min+Df;

fprintf('f_min = %g GHz\n', f_min/1e9);
fprintf('f_max = %g GHz\n', f_max/1e9);
fprintf('Df    = %g MHz\n', Df/1e6);

figure(1);

% Plot the time history of the pulse
subplot(2,1,1)
plot(f,'b')
%set(gca, 'Ytick',[pulse_min pulse_max],'FontSize',6);
%axis([1 Nz pulse_min pulse_max]);
grid on 

% Compute the temporal frequency spectrum by perform fast Fourier transform (FFT
F=fft(f);

% Compute the magnitude of the temporal frequency spectrum
Fabs = abs(F);

% Plot the temporal frequency spectrum of the pulse
subplot(2,1,2)
plot(Fabs,'b')
grid on 

subplot(2,1,1)
title('Excitation Pulse: RC2(t) --- (Time-Domain)','FontSize',11)
axis([1 Nt -1.2 2]);   
set(gca, 'XTick',[1 200e-12/T*Nt Nt/4 Nt/2 3*Nt/4 Nt],'FontSize',10);
set(gca, 'XTickLabel',[0 200 T/4/1e-12 T/2/1e-12 3*T/4/1e-12 T/1e-12],'FontSize',10);
xlabel('t in ps','FontSize',10);    
set(gca, 'YTick',[-1 0 1 2],'FontSize',10);
set(gca, 'YTickLabel',[-1 0 1 2],'FontSize',10);
ylabel('Amplitude of RCn(t)','FontSize',10);

subplot(2,1,2)
title('Excitation Pulse: |RC2(f)| --- (Frequency-Domain)','FontSize',11)

axis([1 Nt 0 12]);   
set(gca, 'XTick',[1 pi*10e9/omega_max*Nt pi*20e9/omega_max*Nt Nt/4 Nt/2 3*Nt/4 Nt],'FontSize',10);
set(gca, 'XTickLabel',[0 10e9/1e9 20e9/1e9 f_max/2/1e9 f_max/1e9 f_min/2/1e9 -Df/1e9],'FontSize',10);
xlabel('f in GHz');    
set(gca, 'YTick',[0 5 10],'FontSize',10);
set(gca, 'YTickLabel',[0 5 10],'FontSize',10);
ylabel('Magnitude of |RCn(f)|');


set(gcf,'PaperType','a4');

fprintf('Generate EPS file: fdtd1d_fig_3.eps\n');
print -depsc2 -r300 fdtd1d_fig_3.eps
fprintf('Generate EPS file: fdtd1d_fig_3.jpg\n');
print -djpeg -r150 fdtd1d_fig_3.jpg

figure(1)
colormap(gray)

fprintf('Generate EPS file: fdtd1d_fig_3_sw.eps\n');
print -deps2 -r300 fdtd1d_fig_3_sw.eps
fprintf('Generate EPS file: fdtd1d_fig_3_sw.jpg\n');
print -djpeg -r150 fdtd1d_fig_3_sw.jpg


