%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% example.m                                                               %
% Script to open hdf5 dataset file in MATLAB and access data              %
% Author: Thomas Haas, KULeuven                                           %
% Date: 19/12/2018                                                        %
% Contact: thomas.haas@kuleuven.be                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc; clear; close all;

% ================================================ %
% USER INPUT: Case settings
%   - Specify path to dataset
%   - Specify wind class ('c0', 'c1' or 'c2')
%   - Specify time serie ('t1', 't2' or 't3')
%   - Turn on/off visualization functions
%   > Retrieve file name
path = '...';
wcase = 'c0';
tserie = 't1';
fname = [path, 'pdbl_', wcase, '_tserie', tserie(end),'.h5'];
% ================================================= %

% ------------------------------------------------- %
% Display file structure:
% h5disp(fname)

% ------------------------------------------------- %
% Scaling:
% The computed velocity fields are non-dimensional. The dimensional field
% is obtained by multiplying the velocity components with the friction
% velocity. The friction velocity is computed from logarithmic velocity 
% profile u = utau/kappa * log( z/z0 ) with kappa=0.4 the von Kàrman 
% constant and z0 the roughness length from the considered class. The time
% step also needs to be scaled with respect to reference friction velocity
% speed and length (ABL height). 
% Tuning parameters:
% u: Reference velocity in [m/s] at reference height
% z: Reference height in [m]
% z0: Roughness length in [m]
% h: Boundary layer height in [m]
% utau: Friction velocity in [m/s]
% Default values:
if (wcase=='c0')
   u=10.0; z=100.0; z0=0.0002; h=1000.0; utau=0.30482;
elseif (wcase=='c1')
   u=9.0; z=100.0; z0=0.03; h=1000.0; utau=0.44380;
elseif (wcase=='c2')
   u=8.0; z=100.0; z0=0.1; h=1000.0; utau=0.46325;
end

% ------------------------------------------------- %
% Metadata: (time and position vectors)
time = (h/utau)*h5read(fname,['/time_array']);
xvec = h*h5read(fname,['/xvec']);
yvec = h*h5read(fname,['/yvec']);
zvec = h*h5read(fname,['/zvec']);
% Number of time steps and grid points
Nt = length(time);
Nx = length(xvec);
Ny = length(yvec);
Nz = length(zvec);
% Temporal and spatial resolution
dt = time(2) - time(1);
dx = xvec(2) - xvec(1);
dy = yvec(2) - yvec(1);
dz = zvec(2) - zvec(1);
% Time period and domain size
Tp = time(end);
Lx = xvec(end);
Ly = yvec(end);
Lz = zvec(end);
% Compute grid
[X,Y,Z] = meshgrid(xvec,yvec,zvec);

% ------------------------------------------------- %
% Load velocity field:
% Each velocity component is located in a specific group:
% - The axial flow velocity component in '/u/'
% - The spanwise flow velocity component in '/v/'
% - The vertical flow velocity component in '/w/'
% Each time step is associated to a specific time stamp
% Example: Retrieve the axial velocity field at t=600s (1st time serie).
nt = 601;
grpname = ['/u/'];
tstamp = ['field_t_',num2str(time(nt)/(h/utau),'%.8f')];
u = h5read(fname,[grpname,tstamp]);

% ------------------------------------------------- %
% Plot slices of the axial velocity field at t=600s. 
% Scale data
U = utau*permute(u,[2,1,3,4]);

% Rearrange arrays (plot 1) 
Xp1 = squeeze(X(31,:,:));
Zp1 = squeeze(Z(31,:,:));
Up1 = squeeze(U(31,:,:));
% Rearrange arrays (plot 2)
Yp2 = squeeze(Y(:,51,:));
Zp2 = squeeze(Z(:,51,:));
Up2 = squeeze(U(:,51,:));
% Rearrange arrays (plot 3)
Xp3 = squeeze(X(:,:,41));
Yp3 = squeeze(Y(:,:,41));
Up3 = squeeze(U(:,:,41));

% Plot 1: XZ-plane (y = 0)
% Create figure and axes
fig = figure(1);
fig.Units = 'centimeters';
fig.Position(1:4) = [1 1 16 10];
ax = axes('units', 'centimeters', 'position', [2.5 2 10 7]);
% Plot data
pcolor( Xp1, Zp1, Up1)
% Axes layout
xlabel('X [m]')
ylabel('Z [m]')
set(ax,'XLim',[xvec(1) xvec(end)],'XTick',[linspace(xvec(1),xvec(end),5)])
set(ax,'YLim',[zvec(1) zvec(end)],'YTick',[linspace(zvec(1),zvec(end),5)])
set(ax,'FontSize', 12)
% Appearance and colorbar
shading interp
caxis([5 13])
cb = colorbar;    
cb.Location = 'manual';
cb.Units = 'centimeters';
cb.Position = [13 2 0.5 7];
cb.Label.String = 'U [m/s]';
cb.FontSize = 12;
% Title
t_str = ['Axial velocity u at time t=',num2str(time(nt),'%.0f'),'s in XZ-plane (Side view)'];
title(t_str)

% Plot 2: YZ-plane (x = 500)
% Create figure and axes
fig = figure(2);
fig.Units = 'centimeters';
fig.Position(1:4) = [1 1 16 10];
ax = axes('units', 'centimeters', 'position', [2.5 2 10 7]);
% Plot data
pcolor( Yp2, Zp2, Up2)
% Axes layout
xlabel('Y [m]')
ylabel('Z [m]')
set(ax,'XLim',[yvec(1) yvec(end)],'XTick',[linspace(yvec(1),yvec(end),5)])
set(ax,'YLim',[zvec(1) zvec(end)],'YTick',[linspace(zvec(1),zvec(end),5)])
set(ax,'FontSize', 12)
% Appearance and colorbar
shading interp
caxis([5 13])
cb = colorbar;    
cb.Location = 'manual';
cb.Units = 'centimeters';
cb.Position = [13 2 0.5 7];
cb.Label.String = 'U [m/s]';
cb.FontSize = 12;
% Title
t_str = ['Axial velocity u at time t=',num2str(time(nt),'%.0f'),'s in YZ-plane (Front view)'];
title(t_str)

% Plot 3: XY-plane (z = 402.5)
% Create figure and axes
fig = figure(3);
fig.Units = 'centimeters';
fig.Position(1:4) = [1 1 16 10];
ax = axes('units', 'centimeters', 'position', [2.5 2 10 7]);
% Plot data
pcolor( Xp3, Yp3, Up3)
% Axes layout
xlabel('X [m]')
ylabel('Y [m]')
set(ax,'XLim',[xvec(1) xvec(end)],'XTick',[linspace(xvec(1),xvec(end),5)])
set(ax,'YLim',[yvec(1) yvec(end)],'YTick',[linspace(yvec(1),yvec(end),5)])
set(ax,'FontSize', 12)
% Appearance and colorbar
shading interp
caxis([5 13])
cb = colorbar;    
cb.Location = 'manual';
cb.Units = 'centimeters';
cb.Position = [13 2 0.5 7];
cb.Label.String = 'U [m/s]';
cb.FontSize = 12;
% Title
t_str = ['Axial velocity u at time t=',num2str(time(nt),'%.0f'),'s in XY-plane (Top view)'];
title(t_str)