Need to extract boundary layer thickness from 2D planar PIV data.

This is 2D Planar PIV data at a freestream velocity of ~ 13 m/s.
My code currently selects a velocity profile where max(profile_u) = 12.41 m/s. I want to use a velocity profile where max(profile_u) is closer to 13 m/s.
Scatter(u,y) will show all the velocity profiles I am working with.
clc;
close all;
clear all;

%% BASELINE
% Load the data
opts = detectImportOptions(‘FlatPlate_Baseline_avg0001(in).csv’);
opts.VariableNamingRule = ‘preserve’;
data = readtable(‘FlatPlate_Baseline_avg0001(in).csv’, opts);

% Extract variables
x = data.x / 1000; % Convert to meters
y = data.y / 1000; % Convert to meters
u = data.u;
v = data.v;

%%
% Filter data for x values between 0.05m to 0.08m
filter_indices = (x >= 0.05) & (x <= 0.08);
x_filtered = x(filter_indices);
y_filtered = y(filter_indices);
u_filtered = u(filter_indices);
v_filtered = v(filter_indices);

% Find unique x and y values within the filtered range
unique_x = unique(x_filtered);
unique_y = unique(y_filtered);

nu = 1.5e-5;
kappa = 0.41;
B = 5.0;
profile_spacing = 2 / 1000;

profile_x_positions = unique_x(unique_x >= 0.05);

wall_shear_stress = zeros(size(profile_x_positions));
u_tau_values = zeros(size(profile_x_positions));
delta_99_values = zeros(size(profile_x_positions));
delta_star_values = zeros(size(profile_x_positions)); % Displacement thickness
theta_values = zeros(size(profile_x_positions)); % Momentum thickness

for i = 1:length(profile_x_positions)
profile_indices = find(abs(x – profile_x_positions(i)) < 1e-6);
profile_y = y(profile_indices);
profile_u = u(profile_indices);

valid_indices = profile_y > 0 & profile_y * max(profile_u) / nu > 1e-6;
profile_y = profile_y(valid_indices);
profile_u = profile_u(valid_indices);

if isempty(profile_y) || isempty(profile_u)
fprintf(‘No valid data points at x = %f meters. Skipping…n’, profile_x_positions(i));
continue;
end

% Estimate u_inf (freestream velocity) as the max value of the profile
u_inf = max(profile_u);

% f(eta) = u / u_inf
f_eta = profile_u / u_inf;
[f_eta_unique, unique_indices] = unique(f_eta);
profile_y_unique = profile_y(unique_indices);

% Find y_99 where f(eta) = 0.99
y_99 = interp1(f_eta_unique, profile_y_unique, 0.99, ‘linear’, ‘extrap’);
eta_99 = 0.99;

% delta_99
delta_99 = y_99 / eta_99;
delta_99_values(i) = delta_99;

% displacement thickness (delta_star)
delta_star = abs(trapz(profile_y, (1 – f_eta)));
delta_star_values(i) = delta_star;

% momentum thickness (theta)
theta = abs(trapz(profile_y, f_eta .* (1 – f_eta)));
theta_values(i) = theta;
Re_theta = (u_inf * theta)/nu;

% Calculate wall shear stress and friction velocity
fun_u_tau = @(u_tau) mean(profile_u – (u_tau / kappa) * log(profile_y * u_tau / nu + 1) – B * u_tau);
try
options = optimset(‘Display’, ‘off’);
u_tau = fzero(fun_u_tau, [0.01, 10], options);
u_tau_values(i) = u_tau;
catch err
fprintf(‘Cannot converge at x = %f meters: %sn’, profile_x_positions(i), err.message);
continue;
end

rho = 1.225;
tau_wall = rho * u_tau^2;

wall_shear_stress(i) = tau_wall;

% % Plot y/delta vs u/u_inf
% figure;
% y_over_delta = profile_y / delta_99;
% u_over_u_inf = profile_u / u_inf;
% plot(u_over_u_inf, y_over_delta, ‘-o’);
% title(sprintf(‘Velocity Profile at x = %.3f meters’, profile_x_positions(i)));
% xlabel(‘u / u_{infty}’);
% ylabel(‘y / delta’);
% grid on;
end

% Display the results
results_table = table(profile_x_positions, delta_99_values, delta_star_values, theta_values);
disp(results_table);
disp([‘Momentum Thickness Reynolds Number (Re_theta): ‘, num2str(Re_theta)]);This is 2D Planar PIV data at a freestream velocity of ~ 13 m/s.
My code currently selects a velocity profile where max(profile_u) = 12.41 m/s. I want to use a velocity profile where max(profile_u) is closer to 13 m/s.
Scatter(u,y) will show all the velocity profiles I am working with.
clc;
close all;
clear all;