1 % PIVLab frame data format
 2 % This section describes the 3D data structure used in PIVLab, where the data is organized as a 3D matrix
 3 % with dimensions corresponding to rows (mRows), columns (nCols), and frames (nFrame).
 4 
 5 % Clear workspace, close figures, and initialize
 6 clc; clear; close all;
 7 
 8 % Load PIVLab data from a specified file
 9 casePath = '..';
10 data = load(fullfile(casePath, 'normalized_data', 'PIVlab.mat'));
11 
12 % Extract spatial mesh grids and display measurement domain information
13 xmesh = data.x{1, 1}; X = xmesh(1, :);
14 ymesh = data.y{1, 1}; Y = ymesh(:, 1);
15 fprintf(['Measure domain: \n'...
16     'X from %.5f to %.5f, nX = %d, dX = %.5f.\n'...
17     'Y from %.5f to %.5f, nY = %d, dY = %.5f.\n'], ...
18     X(1), X(end), length(X), mean(diff(X)), Y(end), Y(1), length(Y), mean(abs(diff(Y))));
19 
20 % Check if filtered data is available; otherwise, use original data
21 if isempty(data.u_filtered{1}) && isempty(data.v_filtered{1})
22     warning('Data not filtered in PIVLab. Mat file may contain abnormal data.')
23     u_filtered = data.u_original;
24     v_filtered = data.v_original;
25     typevector = data.typevector_original;
26 else
27     u_filtered = data.u_filtered;
28     v_filtered = data.v_filtered;
29     typevector = data.typevector_filtered;
30 end
31 
32 % Convert cell arrays to 3D matrices for further processing
33 u_filtered = cat(3, u_filtered{:});
34 v_filtered = cat(3, v_filtered{:});
35 typevector = cat(3, typevector{:});
36 
37 % Remove invalid, NaN, and infinite data based on typevector
38 % Valid vectors are identified by typevector values (0: masked, 1: valid, 2: erroneous)
39 validVector = (typevector ~= 0);
40 u_filtered = u_filtered .* validVector .* (~isnan(u_filtered)) .* (~isinf(u_filtered));
41 v_filtered = v_filtered .* validVector .* (~isnan(v_filtered)) .* (~isinf(v_filtered));
42 nValidFrame = sum(validVector, 3); % Count valid frames for each cell
43 
44 % Statistical processing
45 % Compute time-averaged velocity for each cell
46 U_t = sum(u_filtered, 3) ./ nValidFrame;
47 V_t = sum(v_filtered, 3) ./ nValidFrame;
48 
49 % Compute double-averaged velocity (spatial average in x-direction)
50 U_xt = mean(U_t, 2);
51 V_xt = mean(V_t, 2);
52 
53 % Compute turbulent velocity components (fluctuations)
54 u_pri = u_filtered - U_t; % u'
55 v_pri = v_filtered - V_t; % v'
56 
57 % Compute second moments (Reynolds shear stress and turbulent kinetic energy)
58 uv = u_pri .* v_pri; % u'v'
59 uu = u_pri.^2; % u'u'
60 vv = v_pri.^2; % v'v'
61 
62 % Compute third moments (skewness-related terms)
63 uuu = u_pri.^3;
64 vvv = v_pri.^3;
65 
66 % Compute time-averaged second and third moments
67 uv_t  = sum(uv, 3) ./ nValidFrame;
68 uu_t  = sum(uu, 3) ./ nValidFrame;
69 vv_t  = sum(vv, 3) ./ nValidFrame;
70 uuu_t = sum(uuu, 3) ./ nValidFrame;
71 vvv_t = sum(vvv, 3) ./ nValidFrame;
72 
73 % Compute double-averaged second and third moments (spatial average in x-direction)
74 uv_xt  = mean(uv_t, 2); % u'v'
75 RSS    = -1 * uv_xt; % Reynolds shear stress
76 uu_xt  = mean(uu_t, 2); % u'u'
77 vv_xt  = mean(vv_t, 2); % v'v'
78 uuu_xt = mean(uuu_t, 2); % u'u'u'
79 vvv_xt = mean(vvv_t, 2);
80 
81 % Compute root mean square (RMS) values for turbulence strength
82 u_rms = mean(sqrt(uu_t), 2); % RMS of u
83 v_rms = mean(sqrt(vv_t), 2); % RMS of v
84 
85 % Save statistical results to files
86 save(fullfile(casePath, 'figure_data', 'u_stat.mat'), 'xmesh', 'X', 'ymesh', 'Y', 'U_xt', 'V_xt', 'RSS', 'uu_xt', 'vv_xt', 'u_rms', 'v_rms', 'uuu_xt', 'vvv_xt');
87 save(fullfile(casePath, 'figure_data', 'u4pxx.mat'), '-v7.3', 'xmesh', 'X', 'ymesh', 'Y', 'U_t', 'V_t', 'u_pri', 'v_pri');
88