%function [] = hw
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%
%% 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: %%
%% hw.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 %%
%% %%
%% %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
global hw_num_pos
global mmhg_pa
global hw_p_inf
global rho_hg
global g
global nu
global hw_folder_name
global precal_fl_name
global postcal_fl_name
global dat_fl_name
global zpt_fl_name
global fl_name
global hw_p_cal
global hw_t_cal
global hw_e_cal
global hw_t_inf
global hw_u_cal
global hw_e_fit
global hw_u_fit
global hw_ubar
global hw_ustd
global hw_uskew
global hw_ukurt
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
global re_delta99
global re_dstar
global re_theta
global rho_hg g rho_air nu z_dat hw_ubar
% Define Conastants
hw_num_pos = 50;
mmhg_pa = 133.322368421; % Conversion Factor, mm HG to Pa
hw_p_inf = 750.900024*mmhg_pa; % Change atmospheric pressure here
rho_hg = 13579.04; % Density of Mercury, kg/m3
g = 9.8; % Acceleration due to gravity, m/s^2
nu = 1.51e-5; % Kinematic Viscosity, m^2/s (computed by interpolating for a given temperature, from the web
% Preston tube data for u_tau computations, usually taken from the HW data
dP_pres = 0.0353; % Preston tube reading for utau computatissons
diapres = 0.001; % Preston Tube diameter, m
%hw_folder_name = 'C:\Documents and Settings\Jazz\My Documents\my research\blp091706_6ms\';
hw_folder_name= 'E:\Anurag\Laptop\Research\blp091706_6ms\';
precal_fl_name = 'precal091706.dat';
postcal_fl_name = 'postcal091706.dat';
dat_fl_name = 'blp091706_6ms_pos'; % note that its only the partial file name and ch & position are added to it
zpt_fl_name = 'zpoints091706.dat';
fl_name = [hw_folder_name,'calibration\',precal_fl_name];
fid = fopen(fl_name,'r');
disp('Computing calibration data ...');
for i=1:15 % Skip frist 15 lines
fgetl(fid);
end
temp = fscanf(fid,'%g %g %g %g %g %g %g %g %g %g',[10,inf]);
hw_p_cal = temp(3,:); % pressure in mmHg
hw_t_cal = mean(temp(5,:).*10); % mean temperature
hw_e_cal = temp(7,:); % hot wire voltage corresponding to the velocity
hw_t_inf = hw_t_cal + 273.15; % temperature in deg Kelvin
rho_air = hw_p_inf/(287.3*hw_t_inf); % Density of air @ , kg/m^3
hw_u_cal = sqrt(2*hw_p_cal*rho_hg*g*1.0e-3/rho_air);
hw_e_fit = polyfit(hw_e_cal, hw_u_cal, 3);
hw_u_fit = polyval(hw_e_fit,hw_e_cal);
%%% --------------------------------------------------------------------------------------------------------------
%%% Following peice of commented out code needs some modification for correct post calibration computations.
%%% --------------------------------------------------------------------------------------------------------------
%hw_t_inf_pre = hw_t_cal + 273.15; % temperature in deg Kelvin
%
%rho_air = hw_p_inf/(287.3*hw_t_inf_pre); % Density of air @ , kg/m^3
%
%
%hw_u_cal = sqrt(2*hw_p_cal*rho_hg*g*1.0e-3/rho_air);
%
%
%hw_e_fit_pre = polyfit(hw_e_cal, hw_u_cal, 3);
%hw_u_fit_pre = polyval(hw_e_fit_pre,hw_e_cal);
%
%
%
%% post calibration data
%
%clear hw_p_cal hw_t_cal hw_e_cal hw_u_cal;
%fl_name = [hw_folder_name,'calibration\',postcal_fl_name];
%fid = fopen(fl_name,'r');
%for i=1:15 % Skip frist 15 lines
% fgetl(fid);
%end
%
%temp = fscanf(fid,'%g %g %g %g %g %g %g %g %g %g',[10,inf]);
%
%hw_p_cal = temp(3,:); % pressure in mmHg
%hw_t_cal = mean(temp(5,:).*10); % mean temperature
%
%hw_e_cal = temp(7,:); % hot wire voltage corresponding to the velocity
%
%hw_t_inf_post = hw_t_cal + 273.15; % temperature in deg Kelvin
%
%rho_air = hw_p_inf/(287.3*hw_t_inf_post); % Density of air @ , kg/m^3
%
%
%hw_u_cal = sqrt(2*hw_p_cal*rho_hg*g*1.0e-3/rho_air);
%
%
%hw_e_fit_post = polyfit(hw_e_cal, hw_u_cal, 3);
%hw_u_fit_post = polyval(hw_e_fit_post,hw_e_cal);
%
%
%
%--------------- end calibration -----------------
disp('Reading position data');
fl_name = [hw_folder_name,'data\',zpt_fl_name];
fid = fopen(fl_name,'r');
z_dat = fscanf(fid,'%g',[1,inf]);
fclose(fid);
% Begin reading the hotwire data and compute the profile
% Read hot wire data files
hw_ubar = squeeze(zeros(1,hw_num_pos));
hw_ustd = squeeze(zeros(1,hw_num_pos));
hw_uskew = squeeze(zeros(1,hw_num_pos));
hw_ukurt = squeeze(zeros(1,hw_num_pos));
disp('----------');
for i=1:hw_num_pos
disp(['Processing data, position #',int2str(i)]);
fl_name = [hw_folder_name,'data\',dat_fl_name,int2str(i),'_ch2.bin'];
fid = fopen(fl_name,'r');
[a] = fread(fid,inf,'short');
e = a.*20./(2.^16); % convert the numbers to flaoting point precision (20 is p2p voltage setting, 2^16 is because of 16 bit processor on batchelor
t(i)= mean(e.*10);
fclose(fid);
fl_name = [hw_folder_name,'data\',dat_fl_name,int2str(i),'_ch3.bin'];
fid = fopen(fl_name,'r');
[a] = fread(fid,inf,'short');
e = a.*20./(2.^16); % convert the numbers to flaoting point precision (20 is p2p voltage setting, 2^16 is because of 16 bit processor on batchelor
utemp = polyval(hw_e_fit,e);
%utemp_pre = polyval(hw_e_fit_pre ,e);
%utemp_post = polyval(hw_e_fit_post,e);
%utemp = (utemp_pre - utemp_post)*(t(i) - hw_t_inf_post)/(hw_t_inf_pre - hw_t_inf_post) + utemp_post;
hw_ubar(i) = mean(utemp);
hw_ustd(i) = std(utemp);
hw_uskew(i) = mean((utemp - hw_ubar(i)).^3)./(hw_ustd(i).^3);
hw_ukurt(i) = mean((utemp - hw_ubar(i)).^4)./(hw_ustd(i).^4);
fclose(fid);
end
