function []=hwproc()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%
%% Copyright (C) 2005-2007 Anurag Singh %%
%% %%
%% This program/code snippet/file (hence forth refered as "library") %%
%% is free software; you can redistribute it and/or %%
%% modify it under the terms of the GNU Lesser General Public %%
%% License as published by the Free Software Foundation; either %%
%% version 2.1 of the License, or (at your option) any later version. %%
%% %%
%% This library is distributed in the hope that it will be useful, %%
%% but WITHOUT ANY WARRANTY; without even the implied warranty of %%
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU %%
%% Lesser General Public License for more details. %%
%% %%
%% You should have received a copy of the GNU Lesser General Public %%
%% License along with this library; if not, write to the Free Software %%
%% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA %%
%% %%
%% FILENAME: %%
%% hwproc.m %%
%% %%
%% AUTHOR: %%
%% %%
%% Anurag Singh, %%
%% MS, 2007 %%
%% Aerospace Engineering & Mechanics %%
%% University of Minnesota - Twin Cities. %%
%% Minneapolis, MN 55455 (USA) %%
%% %%
%% (currently working at: Exa Corporation, Burlington, MA 01803) %%
%% %%
%% CONTACT/EMAIL: %%
%% %%
%% anurag@aem.umn.edu %%
%% anurag9@gmail.com %%
%% %%
%% SOURCE CONTROL INFORMATION: %%
%% None (since I was planning on putting it under source control since it has %%
%% reached the critical file system size. Would be a good thing to put it under %%
%% source control while making changes. %%
%% %%
%% DESCRIPTION: %%
%% %%
%% Need to add %%
%% %%
%% %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global dP_pres diapres u_tau tau_w kappa const_a hw_z_pos hw_z_plus hw_u_plus
global hw_u_plus_calc hw_u_inf u_99 delta_99 delta_star theta_star shape_f mixed_s colew
global re_tau re_delta99 re_dstar re_theta
global rho_hg g rho_air nu z_dat hw_ubar hw_num_pos hw_u_cal hw_e_cal hw_u_fit hw_ustd
% Note preston tube measurements seem to be wrong so using clausen chart method
% Preston tube and other calculations
dp = rho_hg*g*0.001*dP_pres; % 0.0353 is mmHg which was measured
% Patel(JFM) Method - Doesn't work for wind tunnel for some reason
%y_star = x_star/2 + 0.037;
x_star = log10(dp*diapres^2/(4*rho_air*nu*nu));
y_star = 0.8287 - 0.1381*x_star + 0.1437*x_star^2 - 0.0060*x_star^3;
%str = sprintf('y+2*log(1.95*y + 4.10) - %g',x_star); y_star = solve(str)
tau_w = 4*rho_air*nu^2*10^(y_star)/diapres^2;
u_tau = sqrt(tau_w/rho_air);
disp(['u_tau computed from from Preston tube measurements = ',num2str(u_tau)]);
% Law of the wall computations
kappa = 0.41;
const_a = 5.0;
epsilon = 0.0002;
hw_z_pos = z_dat*0.001 + epsilon; % hw_z_pos corrected for the initial distance from the wall
% Clausen Chart Method (Works better in general
f = inline('u_tau*( log(hw_z_pos*u_tau/1.5100e-005)/0.41 + 5.0)','u_tau','hw_z_pos');
u_tau = lsqcurvefit(f,0.25,hw_z_pos(16:29),hw_ubar(16:29));
tau_w = u_tau*u_tau*rho_air;
disp(['u_tau computed from from Clausen Chart method (this is the value being used) = ',num2str(u_tau)]);
hw_z_plus = hw_z_pos*u_tau/nu;
hw_u_plus = hw_ubar/u_tau;
hw_u_plus_calc = log(hw_z_plus)/kappa + const_a;
% Compute Boundary Layer Parameters
% BL Thickness, delta_99
hw_u_inf = hw_ubar(hw_num_pos);
delta_99 = 0;
u_99 = 0.99*hw_ubar(hw_num_pos);
for i=1:hw_num_pos
if(hw_ubar(i) > u_99)
break
end
end
delta_99 = hw_z_pos(i-1) + (hw_z_pos(i) - hw_z_pos(i-1)) * (u_99 - hw_ubar(i-1))/(hw_ubar(i) - hw_ubar(i-1));
%Displacement Thickness
delta_star = trapz(hw_z_pos, (1 - hw_ubar/hw_u_inf));
%Momentum Thickness
theta_star = trapz(hw_z_pos, (hw_ubar/hw_u_inf).*(1 - hw_ubar/hw_u_inf));
% Shape Factor
shape_f = delta_star/theta_star;
% Coles Wake Parameter
colew = 0;
for i=2:hw_num_pos
colew = colew + (hw_ubar(i)/hw_u_inf) * (1 - hw_ubar(i)/hw_u_inf) * (hw_z_pos(i) - hw_z_pos(i-1));
end
% Mixed Scaling
mixed_s = hw_u_inf/u_tau;
% Reynolds Number
re_delta99 = hw_u_inf*delta_99/nu;
re_dstar = hw_u_inf*delta_star/nu;
re_theta = hw_u_inf*theta_star/nu;
re_tau = u_tau*delta_99/nu;
% Print values
disp('-------------------------------------------------------------------------');
disp(['Free Stream Velocity hw_u_inf (m/s) = ' num2str(hw_u_inf)]);
disp(['Displacement Thickness, delta_star (mm) = ' num2str(delta_star*1000)]);
disp(['Momentum Thickness, theta (mm) = ' num2str(theta_star*1000)]);
disp(['Shape Factor, delta_star/theta = ' num2str(shape_f)]);
disp(['Skin Friction velocity, u_tau = ' num2str(u_tau)]);
disp(['Mixed Scaling(hw_u_inf/u_tau) = ' num2str(mixed_s)]);
disp(['Boundary Layer Thickness, delta_99 (mm) = ' num2str(delta_99*1000)]);
disp(['Re_theta = ' num2str(re_theta)]);
disp(['Re_tau = ' num2str(re_tau)]);
disp('-------------------------------------------------------------------------');
% Write all the boundary layer data to a .mat file
save ../piv/pivblp.mat hw_u_inf delta_star theta_star shape_f u_tau mixed_s delta_99 re_theta re_tau hw_num_pos hw_z_plus hw_u_plus hw_u_plus_calc hw_u_inf hw_ustd
