function [uh,Xlag] = read_uh_HiWinGS_2013(dfl,ddd,hhh,w_psd)
%{
reads Univ Hawaii R2A sonic wind and motion data specified by dfl
and returns array with the following 11 columns:

1   jd_ref      10 Hz timestamp
2   U           x axis wind velocity, m/s
3   V           y axis wind velocity, m/s
4   W           z axis wind velocity, m/s
5   Tsonic      sonic temperature, C
6   accx        x axis acceleration, m/s2
7   accy        y axis acceleration, m/s2
8   accz        z axis acceleration, m/s2
9   ratex       x axis rotation rate, deg/s
10  ratey       y axis rotation rate, deg/s
11  ratez       z axis rotation rate, deg/s

input parameters: dfl = string path to scs file
                  ddd = julian date
                  hhh = hour
                  w_psd = psd vertical wind variable for lag estimate
%}
%% initialize
disp(['Reading uh file for hour ',int2str(hhh)]);
flist = fopen(dfl,'r');

st1 = [];               % data array
temp = {};              % cell array for textscan output

%% read son data
while feof(flist)==0,   % read entire file into cell array
    try
        % temp will be 1x11 cell array, cells will be ~36000x1 arrays
        temp =textscan(flist,'%f %*f %*f %*f %*f %*f %f %f %f %f %f %f %f %f %f %f %*[^\n]', 'delimiter', ',', 'headerlines', 1,  'emptyvalue', NaN);
        st1 = [st1 cell2mat(temp)'];  % 11x36000 array
    catch
        for ii=1:11  % length of last cell reflects missing values from any/all fields
            if length(temp{1,ii})~=length(temp{1,11})
                temp{1,ii}(length(temp{1,11})+1) = [];  %truncate length of all cells to length(temp{1,11})
            end;
        end;
        if ~isempty(temp{1,1});
            st1 = [st1 cell2mat(temp)']; % append truncated cells to st1 and return for more
        end;
    end;
end;
fclose (flist);
tsons = 1 + st1(1,:)/86400;    % time in decimal yearday

%% scale wind to physical units and rotate
u = (st1(2,:)-2.5)*12;	%  m/s
v = (st1(3,:)-2.5)*12; 
W = (st1(4,:)-2.5)*12;

% rotate 30 deg
uv = [u;v];
theta = (-30*(pi/180));
Tm = [(cos(theta)) (-sin(theta));(sin(theta)) (cos(theta))];
UVrot = Tm*uv;
U = -UVrot(1,:);
V = UVrot(2,:);

% Replace unrealistic values with means
U(U>35) = nanmean1(U);
V(V>35) = nanmean1(V);
W(W>35) = nanmean1(W);

%% scale motion variables and despike

%ratex : port up is positive phi
%ratey : bow down is positive theta
%ratez : bow to port is positive psi (right-handed)
%accx  : forward is positive 
%accy  : to port is positive
%accz  : up is positive
ratex = st1(7,:)*(1000/50.134)*(pi/180);   % rad/sec
ratey = st1(6,:)*(1000/50.034)*(pi/180);
ratez = -st1(8,:)*(1000/49.975)*(pi/180);
accx = st1(10,:)*(9.80665/7.510);          % m/s2
accy = st1(9,:)*(9.80665/7.492);
accz = -st1(11,:)*(9.80665/3.747);

[ratex,~] = despike2(ratex);
[ratey,~] = despike2(ratey);
[ratez,~] = despike2(ratez);
[accx,~] = despike2(accx);
[accy,~] = despike2(accy);
[accz,~] = despike2(accz);

%% compute sonic temperature, NOAA R2 formula
SpdSnd = st1(5,:)*16 + 290;  % speed of sound in m/s
Tsonic = (SpdSnd.^2 + (U.^2 + W.^2)/2. + V.^2)/403 - 273.15;

%% gain adjustment for UH motion pak based on comparison with psd - 6/27/2014 values
accx = accx/0.94992;    % accx, accy, accz are the accelerations in the 
accy = accy/0.91555;    % bow-stern, port-starboard, and up-down directions,
accz = accz/0.94844;    % respectively
ratex = ratex/1.02158;
ratey = ratey/1.02782;
ratez = ratez/1.04830;

%% Interpolate to same timestamp as psd data,just to be sure
start = (ddd+hhh/24);
delta = double(1.0/864000);
last = start + 35999*delta;
tref = [start:delta:last];
% [~,ii,~] = unique(tsons);
% su = interp1(tsons(ii),U(ii),tref);      su(isnan(su)) = nanmean1(su);
% sv = interp1(tsons(ii),V(ii),tref);      sv(isnan(sv)) = nanmean1(sv);
% sw = interp1(tsons(ii),W(ii),tref);      sw(isnan(sw)) = nanmean1(sw);
% st = interp1(tsons(ii),Tsonic(ii),tref); st(isnan(st)) = nanmean1(st);
% srx = interp1(tsons(ii),ratex(ii),tref); srx(isnan(srx)) = nanmean1(srx);
% sry = interp1(tsons(ii),ratey(ii),tref); sry(isnan(sry)) = nanmean1(sry);
% srz = interp1(tsons(ii),ratez(ii),tref); srz(isnan(srz)) = nanmean1(srz);
% sax = interp1(tsons(ii),accx(ii),tref);  sax(isnan(sax)) = nanmean1(sax);
% say = interp1(tsons(ii),accy(ii),tref);  say(isnan(say)) = nanmean1(say);
% saz = interp1(tsons(ii),accz(ii),tref);  saz(isnan(saz)) = nanmean1(saz);

% UH is exactly 36000 pnts, so don't interp
su = U; sv = V; sw = W; st = Tsonic;
srx = ratex; sry = ratey; srz = ratez;
sax = accx; say = accy; saz = accz;

%% remove sonic spikes in 10-min segments - set rejection criteria in despike2.m
for ii = 1:6000:30001
    [su(ii:ii+5999),~] = despike2(su(ii:ii+5999));         
    [sv(ii:ii+5999),~] = despike2(sv(ii:ii+5999));
    [sw(ii:ii+5999),~] = despike2(sw(ii:ii+5999));
    [st(ii:ii+5999),~] = despike2(st(ii:ii+5999));
end;

%% fix lag between PSD and UH data systems
if all(isnan(w_psd))
    Xlag = 0;
else
    Xlag = findlagPositive(w_psd,sw',150,1); % over 15 sec lag interval
end;
    
% Ajdust UH onto PSD timestamp: shift time series by Xlag points
if (Xlag>=0) % positive Xlag: UH lagged w/respect to PSD
    sax(1:36000-Xlag) = sax(Xlag+1:36000);
    say(1:36000-Xlag) = say(Xlag+1:36000);
    saz(1:36000-Xlag) = saz(Xlag+1:36000);
    srx(1:36000-Xlag) = srx(Xlag+1:36000);
    sry(1:36000-Xlag) = sry(Xlag+1:36000);
    srz(1:36000-Xlag) = srz(Xlag+1:36000);
    su(1:36000-Xlag) = su(Xlag+1:36000);
    sv(1:36000-Xlag) = sv(Xlag+1:36000);
    sw(1:36000-Xlag) = sw(Xlag+1:36000);
    st(1:36000-Xlag) = st(Xlag+1:36000);
else        % negative Xlag: PSD lagged w/respect to UH
    sax(-Xlag+1:36000) = sax(1:36000+Xlag);
    say(-Xlag+1:36000) = say(1:36000+Xlag);
    saz(-Xlag+1:36000) = saz(1:36000+Xlag);
    srx(-Xlag+1:36000) = srx(1:36000+Xlag);
    sry(-Xlag+1:36000) = sry(1:36000+Xlag);
    srz(-Xlag+1:36000) = srz(1:36000+Xlag);
    su(-Xlag+1:36000) = su(1:36000+Xlag);
    sv(-Xlag+1:36000) = sv(1:36000+Xlag);
    sw(-Xlag+1:36000) = sw(1:36000+Xlag);
    st(-Xlag+1:36000) = st(1:36000+Xlag);
end;

uh = [tref',su',sv',sw',st',sax',say',saz',srx',sry',srz'];

