S5_Sig1000_RS.m

%% ************************************************************************
%***** Nortek Signature 1000 Processing - Step 5                      *****
%*****                                                                *****
%***** S5_Sig1000_RS.m                                                *****
%*****                                                                *****
%***** A script to perform the Reynolds decomposition of the radial   *****
%***** beam velocity measurements into a mean and a fluctuating part. *****
%***** The period over which the mean is calculated has to be chosen  *****
%***** according to the flow conditions at hand to remove any trends  *****
%***** from the tide or similar large scale influences.               *****
%***** The variance technique by Lohrmann is then performed to        *****
%***** calculate Reynolds stress components from the variance of      *****
%***** beam velocitiy measurements of opposing beams. The approach    *****
%***** for a 5-beam ADCP by Dewey & Stringer 2007 is followed.        *****
%*****                                                                *****
%***** Input: .mat files generated by S3_Sig1000_QC.m                 *****
%***** Output: .mat files of ADCP data and calculated Reynolds        *****
%*****         stresses                                               *****
%*****                                                                *****
%*****                                                                *****
%***** Prepared by J. Tiede, Januar 2018                              *****
%**************************************************************************

%% Clean workspace
clc; clear; close all;

%% Add tools and functions to MATLAB path.
addpath('C:\Users\j_tie\OneDrive\LuFI\Projekte\LuggagePoint_JournalPaper\70_Matlab\ADCP\NortekSig1000\Functions')
addpath('C:\Users\j_tie\OneDrive\LuFI\Projekte\LuggagePoint_JournalPaper\70_Matlab\tools')

%% Define ADCP base folder.
basefolder = 'C:\Users\j_tie\OneDrive\LuFI\Projekte\LuggagePoint_JournalPaper\60_Data\Nortek_061117\'

%% Check if folder exists.
DataFolder = [basefolder '05_QCed\'];
if ~isdir(DataFolder)
  errorMessage = sprintf('Error: Folder does not exist\n%s', DataFolder);
  uiwait(warndlg(errorMessage));
  return;
end

%% Define folder where files are saved to. 
savepath = [basefolder '06_RS\'];

%% Loop over all .mat files
for p = 1
%% User specified variables.
no = 5 % number of beams
th = 25; % Angle of beams to instrument axis
fs = 8; % Sampling frequency
avt = 10; % Period in minutes for calculating mean velocity for Reynolds decomposition
Ravt1 = 1; % Period in minutes for averaging Reynolds stresses + TKE
Ravt10 = 10; % Period in minutes for averaging Reynolds stresses + TKE
Ravt20 = 20; % Period in minutes for averaging Reynolds stresses + TKE
Ravt60 = 55; % Period in minutes for averaging Reynolds stresses + TKE
Ravt100 = 100; % Period in minutes for averaging Reynolds stresses + TKE
Vavt1 = 1; % Period in minutes for averaging horizontal velocity
Vavt10 = 10; % Period in minutes for averaging horizontal velocity
Vavt20 = 20; % Period in minutes for averaging horizontal velocity
Vavt60 = 55; % Period in minutes for averaging horizontal velocity
Vavt90 = 90; % Period in minutes for averaging horizontal velocity
Vavt300 = 300 % Period in minutes for averaging horizontal velocity
npoints1 = fs * 60 * Ravt1; % 8 Hz * 60 s * Ravt min
npoints10 = fs * 60 * Ravt10; % 8 Hz * 60 s * Ravt min
npoints20 = fs * 60 * Ravt20; % 8 Hz * 60 s * Ravt min
npoints60 = fs * 60 * Ravt60; % 8 Hz * 60 s * Ravt min
npoints100 = fs * 60 * Ravt100; % 8 Hz * 60 s * Ravt min
vpoints1 = fs * 60 * Vavt1; % 8 Hz * 60 s * Ravt min
vpoints10 = fs * 60 * Vavt10; % 8 Hz * 60 s * Ravt min
vpoints20 = fs * 60 * Vavt20; % 8 Hz * 60 s * Vavt min
vpoints60 = fs * 60 * Vavt60; % 8 Hz * 60 s * Vavt min
vpoints90 = fs * 60 * Vavt90; % 8 Hz * 60 s * Vavt min
vpoints300 = fs * 60 * Vavt300 % 8 Hz * 60 s * Vavt min

%% Load quality controlled data from Step 4.
load([DataFolder,'PNCU' num2str(p) '_p_qc.mat'])

%% Orientation to radians
ADCP.phi1 = ADCP.phi1*pi/180
ADCP.phi2 = ADCP.phi2*pi/180
ADCP.phi3 = ADCP.phi3*pi/180

%% Correct compass for magnetic declination
ADCP.phi1 = ADCP.phi1 + 10.59
ADCP.VelDirection = ADCP.VelDirection + 10.59

 %% Find ebb flood etc
[max idmax] = max(ADCP.Press)
[min idmin] = min(ADCP.Press)
 

%% Utilise function RS5 to estimate Reynolds stresses and TKE budget.
[RST.uw,RST.vw,RST.u2,RST.v2,RST.w2,TKE.q,Bf.B1f,Bf.B1m,Bf.B2f,Bf.B2m,Bf.B3f,Bf.B3m,Bf.B4f,Bf.B4m,Bf.B5f,Bf.B5m,Bf.b10f,Bf.Buf,Bf.Bvf,RST.uv,Bf.um,Bf.vm] = RS5(no,avt,th,ADCP.phi1,ADCP.phi2,ADCP.phi3,ADCP.b1,ADCP.b2,ADCP.b3,ADCP.b4,ADCP.b5,ADCP.VelX,ADCP.VelY);


%% Filter for waves and other
% f = 8;
% T = 1/f;
% frequencies = [0.1 0.9];
% [b,a] = butter(5, frequencies./(f/2));
% filterButter = tf(b,a,T);
% figure;
% bode(filterButter);
% 
% t = 1:length(Bf.B1f);
% for i = 1:26
% y = double(Bf.B1f(:,i));
% y(isnan(y)) = 0;
% x(:,i) = y;
% yFiltered(:,i) = filtfilt(b,a,x(:,i));
% end
% figure;
% plot(t,x(:,10))
% hold on
% plot(t,yFiltered(:,20),'r')
% legend('unfiltered','filtered')
%%
figure
plot(Bf.b10f(:,10),'b')
hold on
plot(Bf.B1f(:,10),'r')

%% Signal analysis
fs = 8
t = 0:1/fs:10-1/fs;
y = fft(Bf.B1f(:,10));

n = length(y);          % number of samples
f = (0:n-1)*(fs/n);     % frequency range
power = abs(y).^2/n;    % power of the DFT

plot(f,power)
xlabel('Frequency')
ylabel('Power')
hold on

y0 = fftshift(y);         % shift y values
f0 = (-n/2:n/2-1)*(fs/n); % 0-centered frequency range
power0 = abs(y0).^2/n;    % 0-centered power

plot(f0,power0)
xlabel('Frequency')
ylabel('Power')

%% Rotate stresses to earth coordinates.
RSTenu.uw = RST.uw.*cos(mean(ADCP.phi1)) + RST.vw.*sin(mean(ADCP.phi1));
RSTenu.vw = RST.vw.*cos(mean(ADCP.phi1)) - RST.uw.*sin(mean(ADCP.phi1));

%% Rotate stresses to true velocity vector.
flooddir = (nanmean(nanmean(ADCP.VelDirection)))

U = [ADCP.VelX, ADCP.VelY];                             %create velocity matrix
U(isnan(U)) = 0;
U(isinf(U)) = 0;
U = U - repmat(mean(U,1),length(U),1);  %convert matrix to deviate form
R = U'*U/(length(U)-1);                 %compute covariance matrix (alternatively - cov(U))
[V,lambda]=eig(R);                      %calculate eigenvalues and eigenvectors for covariance matrix

%sort eignvalues in descending order so that major axis is given by first eigenvector
[lambda, ilambda]=sort(diag(lambda),'descend');     %sort in descending order with indices
lambda=diag(lambda);                                %reconstruct the eigenvalue matrix
V=V(:,ilambda);                                     %reorder the eigenvectors
ra = atan2(V(2,1),V(2,2));   %rotation angle of major axis in radians relative to cartesian coordiantes
PA = -ra*180/pi+90;                         %express principal axis in compass coordinates
varxp_PA = diag(lambda(1))/trace(lambda);   %variance captured by principal axis

ebbdir = PA
% flooddir = PA + 180

flooddir = flooddir-90
flooddir = flooddir*pi/180

% RSTuvw.uw = RST.uw*cos(flooddir)-RST.vw*sin(flooddir);
% RSTuvw.vw = RST.uw*sin(flooddir)+RST.vw*cos(flooddir);

RSTuvw.uw = RSTenu.vw.*cos(flooddir)-RSTenu.uw.*sin(flooddir);
RSTuvw.vw = RSTenu.vw.*sin(flooddir)+RSTenu.uw.*cos(flooddir);

RSTuvw.u2 = RST.u2 * cos(flooddir)^2 - RST.v2 * cos(flooddir)^2 - 2 * RST.uv * sin(flooddir) * cos(flooddir)
RSTuvw.v2 = RST.u2 * sin(flooddir)^2 + RST.v2 * cos(flooddir)^2 - 2 * RST.uv * sin(flooddir) * cos(flooddir)
RSTuvw.w2 = RST.w2
% ADCP.u=ADCP.VelX.*cos(flooddir)-ADCP.VelY.*sin(flooddir);
% ADCP.v=ADCP.VelX.*sin(flooddir)+ADCP.VelY.*cos(flooddir);

ADCP.u = ADCP.VelE.*cos(flooddir)-ADCP.VelN.*sin(flooddir);
ADCP.v = ADCP.VelE.*sin(flooddir)+ADCP.VelN.*cos(flooddir);

%% Sign Speed
    U=sqrt(ADCP.VelX.^2+ADCP.VelY.^2);
    Udir=atan2(ADCP.VelX,ADCP.VelY)*180/pi;
flood_heading=146.68; %FROM NORTH
flood_heading=flood_heading + [-90, +90];
flood_heading = flood_heading*pi/180

ADCP.U_sign1 = sign_speed(ADCP.VelX,ADCP.VelY,U,Udir,flood_heading);

%% Calculate friction velocity from lowest bin of streamwise Reynolds stresses (Nezu & Nakagawa).
ADCP.ufric = sqrt((RSTuvw.uw(:,1)*-1))

%% Averaging of Reynolds stresses over Ravt minutes.
% 1 min
for i = 1: size(RSTuvw.uw,2)
uwr = RSTuvw.uw(:,i);
M  = size(uwr, 1) - mod(size(uwr,1), npoints1);
y  = reshape(uwr(1:M), npoints1, []);
RSTuvw.uw_avg1(:,i) = transpose(sum(y, 1) / npoints1);
end
for i = 1: size(RSTuvw.vw,2)
vwr = RSTuvw.vw(:,i);
M  = size(vwr, 1) - mod(size(vwr,1), npoints1);
y  = reshape(vwr(1:M), npoints1, []);
RSTuvw.vw_avg1(:,i) = transpose(sum(y, 1) / npoints1);
end
% 10 min
for i = 1: size(RSTuvw.uw,2)
uwr = RSTuvw.uw(:,i);
M  = size(uwr, 1) - mod(size(uwr,1), npoints10);
y  = reshape(uwr(1:M), npoints10, []);
RSTuvw.uw_avg10(:,i) = transpose(sum(y, 1) / npoints10);
end
for i = 1: size(RSTuvw.vw,2)
vwr = RSTuvw.vw(:,i);
M  = size(vwr, 1) - mod(size(vwr,1), npoints10);
y  = reshape(vwr(1:M), npoints10, []);
RSTuvw.vw_avg10(:,i) = transpose(sum(y, 1) / npoints10);
end
% 20 min
for i = 1: size(RSTuvw.uw,2)
uwr = RSTuvw.uw(:,i);
M  = size(uwr, 1) - mod(size(uwr,1), npoints20);
y  = reshape(uwr(1:M), npoints20, []);
RSTuvw.uw_avg20(:,i) = transpose(sum(y, 1) / npoints20);
end
for i = 1: size(RSTuvw.vw,2)
vwr = RSTuvw.vw(:,i);
M  = size(vwr, 1) - mod(size(vwr,1), npoints20);
y  = reshape(vwr(1:M), npoints20, []);
RSTuvw.vw_avg20(:,i) = transpose(sum(y, 1) / npoints20);
end
% 60 min
for i = 1: size(RSTuvw.uw,2)
uwr = RSTuvw.uw(:,i);
M  = size(uwr, 1) - mod(size(uwr,1), npoints60);
y  = reshape(uwr(1:M), npoints60, []);
RSTuvw.uw_avg60(:,i) = transpose(sum(y, 1) / npoints60);
end
for i = 1: size(RSTuvw.vw,2)
vwr = RSTuvw.vw(:,i);
M  = size(vwr, 1) - mod(size(vwr,1), npoints60);
y  = reshape(vwr(1:M), npoints60, []);
RSTuvw.vw_avg60(:,i) = transpose(sum(y, 1) / npoints60);
end
% 100 min
for i = 1: size(RSTuvw.uw,2)
uwr = RSTuvw.uw(:,i);
M  = size(uwr, 1) - mod(size(uwr,1), npoints100);
y  = reshape(uwr(1:M), npoints100, []);
RSTuvw.uw_avg100(:,i) = transpose(sum(y, 1) / npoints100);
end
for i = 1: size(RSTuvw.vw,2)
vwr = RSTuvw.vw(:,i);
M  = size(vwr, 1) - mod(size(vwr,1), npoints100);
y  = reshape(vwr(1:M), npoints100, []);
RSTuvw.vw_avg100(:,i) = transpose(sum(y, 1) / npoints100);
end

%% Averaging of TKE.q.
% 1 min
for i = 1: size(TKE.q,2)
vwr = TKE.q(:,i);
M  = size(vwr, 1) - mod(size(vwr,1), npoints1);
y  = reshape(vwr(1:M), npoints1, []);
TKE.q_avg1(:,i) = transpose(sum(y, 1) / npoints1);
end
% 10 min
for i = 1: size(TKE.q,2)
vwr = TKE.q(:,i);
M  = size(vwr, 1) - mod(size(vwr,1), npoints10);
y  = reshape(vwr(1:M), npoints10, []);
TKE.q_avg10(:,i) = transpose(sum(y, 1) / npoints10);
end
% 20 min
for i = 1: size(TKE.q,2)
vwr = TKE.q(:,i);
M  = size(vwr, 1) - mod(size(vwr,1), npoints20);
y  = reshape(vwr(1:M), npoints20, []);
TKE.q_avg20(:,i) = transpose(sum(y, 1) / npoints20);
end
% 60 min
for i = 1: size(TKE.q,2)
vwr = TKE.q(:,i);
M  = size(vwr, 1) - mod(size(vwr,1), npoints60);
y  = reshape(vwr(1:M), npoints60, []);
TKE.q_avg60(:,i) = transpose(sum(y, 1) / npoints60);
end
% 90 min
for i = 1: size(TKE.q,2)
vwr = TKE.q(:,i);
M  = size(vwr, 1) - mod(size(vwr,1), npoints100);
y  = reshape(vwr(1:M), npoints100, []);
TKE.q_avg90(:,i) = transpose(sum(y, 1) / npoints100);
end
%% Averaging of u2 v2 w2.
for i = 1: size(RSTuvw.u2,2)
u2r = RSTuvw.u2(:,i);
M  = size(u2r, 1) - mod(size(u2r,1), npoints60);
y  = reshape(u2r(1:M), npoints60, []);
RSTuvw.u2_avg60(:,i) = transpose(sum(y, 1) / npoints60);
end
for i = 1: size(RSTuvw.v2,2)
v2r = RSTuvw.v2(:,i);
M  = size(v2r, 1) - mod(size(v2r,1), npoints60);
y  = reshape(v2r(1:M), npoints60, []);
RSTuvw.v2_avg60(:,i) = transpose(sum(y, 1) / npoints60);
end
for i = 1: size(RSTuvw.w2,2)
w2r = RSTuvw.w2(:,i);
M  = size(w2r, 1) - mod(size(w2r,1), npoints60);
y  = reshape(w2r(1:M), npoints60, []);
RSTuvw.w2_avg60(:,i) = transpose(sum(y, 1) / npoints60);
end

%% Average of Speed.
% 1 min
for i = 1: size(ADCP.Speed,2)
spe = ADCP.Speed(:,i);
M  = size(spe, 1) - mod(size(spe,1), vpoints1);
y  = reshape(spe(1:M), vpoints1, []);
ADCP.Speed_avg1(:,i) = transpose(sum(y, 1) / vpoints1);
end
% 10 min
for i = 1: size(ADCP.Speed,2)
spe = ADCP.Speed(:,i);
M  = size(spe, 1) - mod(size(spe,1), vpoints10);
y  = reshape(spe(1:M), vpoints10, []);
ADCP.Speed_avg10(:,i) = transpose(sum(y, 1) / vpoints10);
end
% 20 min
for i = 1: size(ADCP.Speed,2)
spe = ADCP.Speed(:,i);
M  = size(spe, 1) - mod(size(spe,1), vpoints20);
y  = reshape(spe(1:M), vpoints20, []);
ADCP.Speed_avg20(:,i) = transpose(sum(y, 1) / vpoints20);
end
% 60 min
for i = 1: size(ADCP.Speed,2)
spe = ADCP.Speed(:,i);
M  = size(spe, 1) - mod(size(spe,1), vpoints60);
y  = reshape(spe(1:M), vpoints60, []);
ADCP.Speed_avg60(:,i) = transpose(sum(y, 1) / vpoints60);
end
% 90 min
for i = 1: size(ADCP.Speed,2)
spe = ADCP.Speed(:,i);
M  = size(spe, 1) - mod(size(spe,1), vpoints90);
y  = reshape(spe(1:M), vpoints90, []);
ADCP.Speed_avg90(:,i) = transpose(sum(y, 1) / vpoints90);
end
%% Average of vertical Speed.
% 1 min
for i = 1: size(ADCP.b5,2)
spe = ADCP.b5(:,i);
M  = size(spe, 1) - mod(size(spe,1), vpoints1);
y  = reshape(spe(1:M), vpoints1, []);
ADCP.VelUp1_avg1(:,i) = transpose(sum(y, 1) / vpoints1);
end
% 10 min
for i = 1: size(ADCP.b5,2)
spe = ADCP.b5(:,i);
M  = size(spe, 1) - mod(size(spe,1), vpoints10);
y  = reshape(spe(1:M), vpoints10, []);
ADCP.VelUp1_avg10(:,i) = transpose(sum(y, 1) / vpoints10);
end
% 60 min
for i = 1: size(ADCP.b5,2)
spe = ADCP.b5(:,i);
M  = size(spe, 1) - mod(size(spe,1), vpoints60);
y  = reshape(spe(1:M), vpoints60, []);
ADCP.VelUp1_avg60(:,i) = transpose(sum(y, 1) / vpoints60);
end
%% Average of flow direction.
% 1 min
for i = 1: size(ADCP.VelDirection,2)
dir = ADCP.VelDirection(:,i);
M  = size(dir, 1) - mod(size(dir,1), vpoints1);
y  = reshape(dir(1:M), vpoints1, []);
ADCP.VelDirection_avg1(:,i) = transpose(sum(y, 1) / vpoints1);
end
% 10 min
for i = 1: size(ADCP.VelDirection,2)
dir = ADCP.VelDirection(:,i);
M  = size(dir, 1) - mod(size(dir,1), vpoints10);
y  = reshape(dir(1:M), vpoints10, []);
ADCP.VelDirection_avg10(:,i) = transpose(sum(y, 1) / vpoints10);
end
% 20 min
for i = 1: size(ADCP.VelDirection,2)
dir = ADCP.VelDirection(:,i);
M  = size(dir, 1) - mod(size(dir,1), vpoints20);
y  = reshape(dir(1:M), vpoints20, []);
ADCP.VelDirection_avg20(:,i) = transpose(sum(y, 1) / vpoints20);
end
% 30 min
for i = 1: size(ADCP.VelDirection,2)
dir = ADCP.VelDirection(:,i);
M  = size(dir, 1) - mod(size(dir,1), vpoints60);
y  = reshape(dir(1:M), vpoints60, []);
ADCP.VelDirection_avg60(:,i) = transpose(sum(y, 1) / vpoints60);
end
% 90 min
for i = 1: size(ADCP.VelDirection,2)
dir = ADCP.VelDirection(:,i);
M  = size(dir, 1) - mod(size(dir,1), vpoints90);
y  = reshape(dir(1:M), vpoints90, []);
ADCP.VelDirection_avg110(:,i) = transpose(sum(y, 1) / vpoints90);
end
%% Averaging of velocity in streamwise direction u.
% 10 min
for i = 1: size(ADCP.u,2)
dir = ADCP.u(:,i);
M  = size(dir, 1) - mod(size(dir,1), vpoints10);
y  = reshape(dir(1:M), vpoints10, []);
ADCP.u_avg10(:,i) = transpose(sum(y, 1) / vpoints10);
end
% 20 min
for i = 1: size(ADCP.u,2)
dir = ADCP.u(:,i);
M  = size(dir, 1) - mod(size(dir,1), vpoints20);
y  = reshape(dir(1:M), vpoints20, []);
ADCP.u_avg20(:,i) = transpose(sum(y, 1) / vpoints20);
end
% 300 min
for i = 1: size(ADCP.u,2)
dir = ADCP.u(:,i);
M  = size(dir, 1) - mod(size(dir,1), vpoints300);
y  = reshape(dir(1:M), vpoints300, []);
ADCP.u_avg300(:,i) = transpose(sum(y, 1) / vpoints300);
end
%% Depth averaging of velocity in streamwise direction u.
ADCP.u20_davg = nanmean(ADCP.u_avg20)
ADCP.u300_davg = nanmean(ADCP.u_avg300)
ADCP.Time_avg1 = datenum((ADCP.Time(1):minutes(Ravt1):(ADCP.Time(end))))
ADCP.Time_avg10 = datenum((ADCP.Time(1):minutes(Ravt10):(ADCP.Time(end))))
ADCP.Time_avg20 = datenum((ADCP.Time(1):minutes(Ravt20):(ADCP.Time(end))))
ADCP.Time_avg60 = datenum((ADCP.Time(1):minutes(Ravt60):(ADCP.Time(end))))
ADCP.Time_avg100 = datenum((ADCP.Time(1):minutes(Ravt100):(ADCP.Time(end))))

%% Calculate averaged friction velocity.
% 20 min
ADCP.u_star20 = nanmean(0.0025^0.5 * ADCP.u20_davg)
% 300 min
ADCP.u_star300 = nanmean(0.0025^0.5 * ADCP.u300_davg)

%% Calc Uncertainty. Check noise for every deployment.
noise_slanted4 = 0.0734
noise_vert = 0.046
noise = sqrt(0.0734^2+0.046^2)
RSTuvw.unc = sqrt(noise^4/((8*60*60)*sin(2*20)^2))

%% Calculate theoretical value for Reynolds stress profile
z = 14
H = 14
delta_v = 0.5
delta_z = 0.5
visc = 1000
x = ADCP.u_star300 * (1-(z/H));
y = visc .* (delta_v(:,:)./delta_z);
thRS_uw = ADCP.u_star300^2*(1-(z/H)) - visc .* (delta_v./delta_z);
thRS_uw_norm = thRS_uw - nanmin(thRS_uw)
thRS_uw_norm = thRS_uw_norm ./ nanmax(thRS_uw_norm(:))

%% Calculate theoretical TKE.q (Nezu & Nakagawa).
TKE.q_th = 9.56 * ADCP.u_star20.^2 .* exp(-2*(z/H))

%% Save .mat files.
fprintf(1, 'Now saving PNCU%s_RS.mat.\n', num2str(p));
save([savepath 'PNCU' num2str(p) '_RS.mat'],'ADCP','RST','RSTuvw','RSTenu','TKE','M','Bf','thRS_uw','thRS_uw_norm','-v7.3','-nocompression')

%% Clear
% clear ADCP remc1 remc2 remc3 remc4 remc5 rema1 rema2 rema3 rema4 rema5 percbad

end

fprintf(1, 'Done! Saved %s file(s).', num2str(p));

%% Test
figure
plot((ADCP.U_sign1(5000:6000,10)*-1))
hold on
plot((ADCP.Speed(5000:6000,10)),'r')