%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Exercise set 7: Properties of turbulence %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Instructions: Complete the following matlab script to accomplish the
% tasks required by the exercise set. Much of the script is done for you,
% you will need to complete lines marked with 'ReplaceMe'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;        % Clear the  command window
clear;      % Clear all stored variables
close all;  % Close all open figures


% Import the data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u = load('u.out'); % Load u-velocity data
w = load('w.out'); % Load w-velocity data
q = load('q.out'); % Load specific humidity data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   a.  Plot u,w, and q as a function of time.      %
%       Mark the mean of u w and q on each plot.    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Step 1.   Create a time vector, with the same size as the data, 
%           corresponding to a sampling frequency of 20 Hz.

fs = 20;            % Sampling frequency [Hz]
N = length(u);      % Number of samples
t = (0:1:(N-1))/fs; % Time vector

% Step 2.   Determin the mean values of the data

u_avg = mean(u);    % Find the mean value of u;
w_avg = mean(w);    % Find the mean value of w;
q_avg = mean(q);    % Find the mean value of q;

% Print the result to the command window
fprintf('Mean of u = %3.2f [m/s] \n',  u_avg)
fprintf('Mean of w = %3.2e [m/s] \n',  w_avg)
fprintf('Mean of q = %3.2f [g/m^3] \n',q_avg)
fprintf('Total sampling time = %3.2e [s] \n\n',t(N))

% Step 3.   Plot the raw data in terms of time, 

% Create a figure that is 16x16cm with a white background
% plots will appear here.
figure('color','w','unit','centimeters','position',[1 4 16 16]) 

ax{1} = subplot(311);   % Create a subplot, of a 3x1 set of subplots, 
                        % this will be the first, (hence, 311). Save a
                        % reference to the axis to ax{1}
ax{2} = subplot(312);   % Create a subplot, of a 3x1 set of subplots, 
                        % this will be the second, (hence, 312). Save a
                        % reference to the axis to ax{2}
ax{3} = subplot(313);   % Create a subplot, of a 3x1 set of subplots, 
                        % this will be the third, (hence, 313). Save a
                        % reference to the axis to ax{3}

hold(ax{1},'on')    % Allow multiple plots on the same axes  
hold(ax{2},'on')    % Allow multiple plots on the same axes  
hold(ax{3},'on')    % Allow multiple plots on the same axes                    

box(ax{1},'on')     % Display the axes as a box
box(ax{2},'on')     % Display the axes as a box
box(ax{3},'on')     % Display the axes as a box

ylabel(ax{1},'u [m/s]')   % Label the y-axis
ylabel(ax{2},'w [m/s]')   % Label the y-axis
ylabel(ax{3},'q [g/m^3]') % Label the y-axis

xlabel(ax{1},'t [s]')   % Label the x-axis
xlabel(ax{2},'t [s]')   % Label the x-axis
xlabel(ax{3},'t [s]')   % Label the x-axis

[ymin, ymax] = bounds([u w q]); % Find the upper and lower bounds of u
[xmin, xmax] = bounds(t);       % Find the upper and lower bounds of t

axis(ax{1},[xmin xmax ymin(1) ymax(1)]); % Adjust the size of the plot axes;
axis(ax{2},[xmin xmax ymin(2) ymax(2)]); % Adjust the size of the plot axes;
axis(ax{3},[xmin xmax ymin(3) ymax(3)]); % Adjust the size of the plot axes;

plot(ax{1},t,u,'b-');   % Plot u as a function of time
                        % using a solid blue line 'b-'
plot(ax{2},t,w,'b-');   % Plot w as a function of time
                        % using a solid blue line 'b-'
plot(ax{3},t,q,'b-');   % Plot q as a function of time
                        % using a solid blue line 'b-'

plot(ax{1},[t(1) t(N)],u_avg*ones(1,2),'r--'); % Plot u_avg on the same  
                                               % axes as the u-time series
plot(ax{2},[t(1) t(N)],w_avg*ones(1,2),'r--'); % Plot w_avg on the same  
                                               % axes as the w-time series
plot(ax{3},[t(1) t(N)],q_avg*ones(1,2),'r--'); % Plot q_avg on the same  
                                               % axes as the q-time series

legend(ax{1},'u','u_{avg}') % Label the time series and mean values
legend(ax{2},'u','u_{avg}') % Label the time series and mean values
legend(ax{3},'u','u_{avg}') % Label the time series and mean values

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   b.  Verify that the average of u', w', and q' = 0 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract the fluctuating component of each time series, based on
% the principle u(t) = u_avg + u', where u' is the fluctuating component
u_fluct = u - u_avg;
w_fluct = w - w_avg;
q_fluct = q - q_avg;

% Take the mean of the fluctuating component of each time series
u_fluct_avg = mean(u_fluct);
w_fluct_avg = mean(w_fluct);
q_fluct_avg = mean(q_fluct);

% Print the result to the command window
fprintf('Mean of u'' = %3.2e [m/s] \n',  u_fluct_avg)
fprintf('Mean of w'' = %3.2e [m/s] \n',  w_fluct_avg)
fprintf('Mean of q'' = %3.2e [g/m^3] \n\n',q_fluct_avg)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   c.  Plot u', w', and q' and calculate the standard deviation    %
%       of each times series                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Step 1. Plot the time series of u', w', and q'

% Create a figure that is 16x16cm with a white background
% plots will appear here.
figure('color','w','unit','centimeters','position',[1 4 16 16]) 

ax{1} = subplot(311);   % Create a subplot, of a 3x1 set of subplots, 
                        % this will be the first, (hence, 311). Save a
                        % reference to the axis to ax{1}
ax{2} = subplot(312);   % Create a subplot, of a 3x1 set of subplots, 
                        % this will be the second, (hence, 312). Save a
                        % reference to the axis to ax{2}
ax{3} = subplot(313);   % Create a subplot, of a 3x1 set of subplots, 
                        % this will be the third, (hence, 313). Save a
                        % reference to the axis to ax{3}

box(ax{1},'on')     % Display the axes as a box
box(ax{2},'on')     % Display the axes as a box
box(ax{3},'on')     % Display the axes as a box

[ymin, ymax] = bounds([u_fluct w_fluct q_fluct]); % Find the upper and lower bounds of u
[xmin, xmax] = bounds(t); % Find the upper and lower bounds of t

axis(ax{1},[xmin xmax ymin(1) ymax(1)]); % Adjust the size of the plot axes;
axis(ax{2},[xmin xmax ymin(2) ymax(2)]); % Adjust the size of the plot axes;
axis(ax{3},[xmin xmax ymin(3) ymax(3)]); % Adjust the size of the plot axes;

plot(ax{1},t,u_fluct,'b-'); % Plot u' as a function of time
                            % using a solid blue line 'b-'
plot(ax{2},t,w_fluct,'b-'); % Plot w' as a function of time
                            % using a solid blue line 'b-'
plot(ax{3},t,q_fluct,'b-'); % Plot q' as a function of time
                            % using a solid blue line 'b-'

ylabel(ax{1},'u'' [m/s]')   % Label the y-axis
ylabel(ax{2},'w'' [m/s]')   % Label the y-axis
ylabel(ax{3},'q'' [g/m^3]') % Label the y-axis

xlabel(ax{1},'t [s]')   % Label the x-axis
xlabel(ax{2},'t [s]')   % Label the x-axis
xlabel(ax{3},'t [s]')   % Label the x-axis

% Step 2. Calculate the standard deviation 

u_std = std(u); % Calculate the standard deviation of u
w_std = std(w); % Calculate the standard deviation of w
q_std = std(q); % Calculate the standard deviation of q

% Print results to the command window
fprintf('Standard deviation of u = %3.2e [m/s] \n',  u_std)
fprintf('Standard deviation of w = %3.2e [m/s] \n',  w_std)
fprintf('Standard deviation of q = %3.2e [m/s] \n\n',  q_std)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   d.  Determine the turbulent fluxes of momentum u'w'_avg         %
%       and vapor, w'q'_avg                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

uw_avg = mean(u_fluct.*w_fluct); % Determine the turbulent flux of momentum
wq_avg = mean(w_fluct.*q_fluct); % Determine the turbulent flux of vapor

% Print the result to the command window
fprintf('Mean of u''w'' = %3.2e [m^2/s^2] \n',  uw_avg)
fprintf('Mean of w''q'' = %3.2e [g m^-2 s^-1] \n\n',  wq_avg)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   e.  based on the result of (d.), what time of day were these    %
%        measurements taken?                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   f.  Knowing the vertical gradient of the mean water vapor      %
%       concentration at the measurement height (dqdz = -0.5 g m^-4)%
%       obtain the eddy diffusion coefficient for water vapor.      %
%       Compare the value you find to the molecular diffusion       %
%       coefficient Dm = 2.42e-5 m^2/s @ T=20C                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


dqdz= -0.5; % Vertical gradient of the mean water vapor concentration at 
            % the measurement height
Dt= -wq_avg/dqdz;
fprintf('Eddy diffusion coefficient for water vapor: Dt = %3.2e [m^2/s] \n\n',  Dt)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   g.  Compute and plot the aurtocorrelation function of u',       %
%       normalized so the maximum correlation is 1.0. you can use   %
%       MATLAB's built-in autocorr function for this                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u_acf = acf(u_fluct); % Compute the complete autocorrelation function.

% Step 1. Plot the time series of u', w', and q'

% Create a figure that is 16x16cm with a white background
% plots will appear here.
figure('color','w','unit','centimeters','position',[1 4 16 16]) 

ax = axes(gcf);   % Create an axes reference
box(ax,'on')      % Display the axes as a box

[ymin, ymax] = bounds(u_acf); % Find the upper and lower bounds of u
[xmin, xmax] = bounds(t);     % Find the upper and lower bounds of t

axis(ax,[xmin xmax ymin ymax]); % Adjust the size of the plot axes;

plot(ax,t,u_acf,'b-'); % Plot autocorrelation function as a function of 
                       % time-lags using a solid blue line 'b-'

ylabel(ax,'\rho')        % Label the y-axis
xlabel(ax,'\tau [s]')    % Label the x-axis
title(ax, 'Autocorrlation function') % Add title

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   h.  Calculate the integral time scale by integrating the        %
%       autocorrelation function between tau = 0 and when the acf   %
%       first becomes 0                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tau_f = min(t(u_acf<0)); % Time when the autocorrelation function first 
                         % becomes 0;
Ti = trapz(t(t<tau_f),u_acf(t<tau_f)); % Integrate over tau {0,tau_f}

% Print result to the command window
fprintf('Integral time scale Ti = %3.2f [s] \n\n', Ti);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   i.  Estimate the typical size of the largest eddies,            %
%       i.e. calculate the integral length scale.                   %
%       hint: you can use the mean velocity to calculate the        %
%       integral length scale from the integral time scale          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Li = Ti*u_avg; % Calculate the integral length scale

% Print result to the command window
fprintf('Integral length scale Li = %3.2f [m] \n\n', Li);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   j.  Estimate the Reynolds number based on the mean velocity     %
%       and integral length scale. is the low laminar or turbulent? %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Re = u_avg*Li/1.568e-5; % Calculate the Reynolds number

% Print result to the command window
fprintf('Reynolds number Re = %3.2e \n\n', Re); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ACF function                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ACF = acf(u_fluct,varargin)
% ACF = acf(u_fluct,NLags) outputs the autocorrelation function for a
% number of time-lags NLags. 
if ~isempty(varargin)
    NLags = varargin{1};
else
    NLags = length(u_fluct);
end
N = length(u_fluct);
ACF = zeros(1,NLags);
for k=1:NLags
    ACF(k)=sum(u_fluct(1:N-k+1).*u_fluct(k:N));
end
ACF=ACF/ACF(1);
end