% Post-process vortex case

clear
clc
close all

% Read FDS output files
fds_02 = csvread('../vortex_2/vortex_2_line.csv',2,0);
fds_05 = csvread('../vortex_5/vortex_5_line.csv',2,0);
fds_10 = csvread('../vortex_10/vortex_10_line.csv',2,0);
fds_20 = csvread('../vortex_20/vortex_20_line.csv',2,0);
fds_40 = csvread('../vortex_40/vortex_40_line.csv',2,0);

%% Begin analysis of U_z (along z=0) %%
% Analytical solution (U_z along z = 0 and with finest grid)
Rc   = 0.005;           % Core radius of vortex (m)
Umax = 0.5;             % Maximum velocity (m/s)
nu   = 1.0e-5;          % Kinematic viscosity (m^2/s)
t    = 3.0;             % Time to evaluate analytical solution (s)

R = sqrt(2*nu*t+Rc^2);
U = Umax*(((2*nu*t)./(Rc^2))+1).^(-3/2);

% NB: U_z is actually evaluated at (xc,zc+dz/2)
xv   = 0;
zv   = 0;
x(:) = fds_40(:,3);
z    = 0;
U_z  = -(U*(x(:)-xv)/R).*exp( 0.5-0.5*((x(:)-xv).^2+(z-zv)^2)/(R^2) );

% Plots
% Plot Analytical Solution
figure;
plot(fds_40(:,3),U_z,'-ok','LineWidth',2, ...
     'MarkerFaceColor','k','MarkerSize',4);
hold on

% Plot FDS results
% NB: U_z is actually evaluated at (xc,zc+dz/2)=(xc,0)
% OK: curves are evaluated at same location
plot(fds_02(:,3),fds_02(:,4),'--b','LineWidth',2);
hold on
plot(fds_05(:,3),fds_05(:,4),'--g','LineWidth',2);
hold on
plot(fds_10(:,3),fds_10(:,4),'--c','LineWidth',2);
hold on
plot(fds_20(:,3),fds_20(:,4),'--m','LineWidth',2);
hold on
plot(fds_40(:,3),fds_40(:,4),'--r','LineWidth',2);

%% End analysis of U_z (along z=0) %%




