%//************************************************************************
%//    snapshots.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:  December 19, 2001
%//
%//    This program reads the binary snapshot data of the H_x, H_z, and E_y 
%//    component and generates a sequence of output figures on the screen.
%//************************************************************************

clear all;
close all;

fprintf('\nVisualize time-domain snapshots of a incident plane wave diffracted\n');
fprintf('by a slit.\n');

Nx = 200;
fprintf('Number of smaplex in x-direction: Nx = %d\n', Nx);
Nz = 200;
fprintf('Number of smaplex in x-direction: Nz = %d\n', Nz);

field = zeros(Nz,Nx);
help  = zeros(Nx*Nz,1);

fp_Hx = fopen('Hx.bin','rb');
fp_Hz = fopen('Hz.bin','rb');
fp_Ey = fopen('Ey.bin','rb');

fig = figure('Position',[0 0 350 1000]);
figure(1)

for ss=1:15
       
    subplot(3,1,1);
    [help,nread] = fread(fp_Hx,Nx*Nz,'float32');
    fprintf('Number of samples read: Nx*Nz = %d\n', nread)    
    for nz=0:Nz-1
        for nx=1:Nx            
            field(Nz-nz,nx) = help(nx+nz*Nx,1);
        end
    end        
    pcolor(field);
    shading flat
    axis([1 Nx 1 Nz]);
    title('2D TM_y FDTD result: H_x component');
    xlabel('x in cells');
    ylabel('z in cells');
    colorbar;
    grid on;
    
    subplot(3,1,2);
    [help,nread] = fread(fp_Hz,Nx*Nz,'float32');
    fprintf('Number of samples read: Nx*Nz = %d\n', nread)    
    for nz=0:Nz-1
        for nx=1:Nx            
            field(Nz-nz,nx) = help(nx+nz*Nx,1);
        end
    end        
    pcolor(field);
    shading flat
    axis([1 Nx 1 Nz]);
    xlabel('x in cells');
    ylabel('z in cells');
    title('2D TM_y FDTD result: H_z component');
    colorbar;
    grid on;

    subplot(3,1,3);
    [help,nread] = fread(fp_Ey,Nx*Nz,'float32');
    fprintf('Number of samples read: Nx*Nz = %d\n', nread)    
    for nz=0:Nz-1
        for nx=1:Nx            
            field(Nz-nz,nx) = help(nx+nz*Nx,1);
        end
    end        
    pcolor(field);
    shading flat
    axis([1 Nx 1 Nz]);
    xlabel('x in cells');
    ylabel('z in cells');
    title('2D TM_y FDTD result: E_y component');
    colorbar;
    grid on;
    
    pause(0.0001);

end