%//************************************************************************
%//    fig_2.m: MATLAB visualization program
%//***********************************************************************
%//    Program author: Dr.-Ing. René Marklein
%//                    University of Kassel, Department of Electrical and 
%//                    Computer Engineering, Electromagnetic Field Theory
%//                    Wilhelmshöher Allee 71
%//                    D-34121 Kassel
%//                    Tel. 0561 804 6426, marklein@uni-kassel.de
%//
%//    Date of this version:  November 15, 2001
%//
%//    MATLAB program to visualize the output file Ex.bin and Hy.bin of the
%//    1D FDTD code with PEC boundary conditions. The main output figure 
%//    contains 4 subfigures displaying each snapshots of the Ex and Hy
%//    component.
%//************************************************************************

close all;
clear all;

Nt = 500;
Nz = 200;

Dt = 2.5e-12;
Dx = 1.5e-3;

T = Dt*Nt;
Z = Dx*Nz;

Ex = zeros(Nt,Nz);
Hy = zeros(Nt,Nz);

X = zeros(Nt*Nz,1);

figure(1);

fprintf('Read binary file: Ex.bin\n');
fp=fopen('Ex.bin','r');
[X,nread] = fread(fp,Nt*Nz,'float32');
fclose(fp);

for nt=1:Nt
  for nz=1:Nz    
    Ex(nt,nz) = X((nt-1)*Nz+nz);
  end
end

fprintf('Determine min and max of Ex\n');
Ex_max = max(max(Ex));
Ex_min = min(min(Ex));

fprintf('Ex_min = %g Ex_max = %g\n', Ex_min, Ex_max);

fprintf('Read binary file: Hy.bin\n');
fp=fopen('Hy.bin','r');
[X,nread] = fread(fp,Nt*Nz,'float32');
fclose(fp);

for nt=1:Nt
  for nz=1:Nz    
    Hy(nt,nz) = X((nt-1)*Nz+nz);
  end
end

fprintf('Determine min and max of Hy\n');
Hy_max = max(X);
Hy_min = min(X);

fprintf('Hy_min = %g Hy_max = %g\n', Hy_min, Hy_max);

subplot(2,2,1)
title('E_x component');

subplot(2,2,3)
title('H_y component');

for nt=1:6:Nt/5
    subplot(2,2,1)
    plot(Ex(nt,:),'b')    
    set(gca, 'Ytick',[-2 -1 0 1 2],'FontSize',10);
    set(gca, 'YTickLabel',[-2 -1 0 1 2],'FontSize',10);
    xlabel('Amplitude of E_x');
    axis([1 Nz -2 2]);
    set(gca, 'XTick',[1 Nz/4 Nz/2 3*Nz/4 Nz],'FontSize',10);
    set(gca, 'XTickLabel',[0 Z/4 Z/2 3*Z/4 Z],'FontSize',10);
    xlabel('z in m');
    title(['E_x component at nt = ',num2str(nt)],'FontSize',11)
    grid on

    subplot(2,2,3)
    plot(Hy(nt,:),'r')
    set(gca, 'Ytick',[-2 -1 0 1 2],'FontSize',10);
    set(gca, 'YTickLabel',[-2 -1 0 1 2],'FontSize',10);
    xlabel('Amplitude of H_y');
    axis([1 Nz -2 2]);
    set(gca, 'XTick',[1 Nz/4 Nz/2 3*Nz/4 Nz],'FontSize',10);
    set(gca, 'XTickLabel',[0 Z/4 Z/2 3*Z/4 Z],'FontSize',10);
    xlabel('z in m');
    title(['H_y component at nt = ',num2str(nt)],'FontSize',11)
    grid on
    pause(0.005)
end

subplot(2,2,2)
pcolor(Ex)
shading flat
title('E_x component','FontSize',11);    
axis([1 Nz 1 Nt]);
set(gca, 'XTick',[1 Nz/4 Nz/2 3*Nz/4 Nz],'FontSize',10);
set(gca, 'XTickLabel',[0 Z/4 Z/2 3*Z/4 Z],'FontSize',10);
xlabel('z in m');
set(gca, 'Ytick',[1    0.25*Nt 0.5*Nt 0.75*Nt Nt],'FontSize',10);
set(gca, 'YTickLabel',[0 T/4/1e-12 T/2/1e-12 3*T/4/1e-12 T/1e-12],'FontSize',10);
ylabel('t in ps');

colorbar;
grid on;

subplot(2,2,4)
pcolor(Hy)
shading flat
title('H_y component','FontSize',10);
axis([1 Nz 1 Nt]);

set(gca, 'XTick',[1 Nz/4 Nz/2 3*Nz/4 Nz],'FontSize',10);
set(gca, 'XTickLabel',[0 Z/4 Z/2 3*Z/4 Z],'FontSize',10);
xlabel('z in m');

set(gca, 'Ytick',[1    0.25*Nt 0.5*Nt 0.75*Nt Nt],'FontSize',10);
set(gca, 'YTickLabel',[0 T/4/1e-12 T/2/1e-12 3*T/4/1e-12 T/1e-12],'FontSize',10);
ylabel('t in ps');
grid on;

colorbar;
grid on;

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

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

figure(1)
colormap(gray)
fprintf('Generate EPS file: fdtd1d_fig_2_sw.eps\n');
print -deps2 -r150 fdtd1d_fig_2_sw.eps
fprintf('Generate EPS file: fdtd1d_fig_2_sw.jpg\n');
print -djpeg -r150 fdtd1d_fig_2_sw.jpg


