function [son,bad,missing] = read_sonic_SeaState_2015(dfl,ddd,hhh)
%{
reads Metek sonic file specified by dfl and returns arrays with the
following 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

Data in file are 20Hz, output is interpolated to 10Hz.

also returns 6x1 array with the number of bad data points and the number
of missing data lines in 10-min data segments for each sonic

input parameters: dfl = string path to mtx file
                  ddd = julian date
                  hhh = hour
%}
%% initialize
disp(['Reading mtx0 file for hour ',int2str(hhh)]);
flist = fopen(dfl,'r');

st1 = [];               % data array
temp = {};              % cell array for textscan output
bad = zeros(1,6);   % initial value for # bad points in each 10-min seg
missing = zeros(1,6);   % initial value for # missing lines in each 10-min seg

%% read son data
while feof(flist)==0,   % read entire file into cell array
    try
        % temp will be 1x7 cell array, cells will be ~72000x1 arrays
        temp = textscan(flist,'%2f%2f%3f %*2c%f %f %f %f %*[^\n]','delimiter',', ','headerlines',1,'emptyvalue',NaN);
        st1 = [st1 cell2mat(temp)'];  % 7x~72000 array
    catch
        for ii=1:7  % length of last cell reflects missing values from any/all fields
            if length(temp{1,ii})~=length(temp{1,7})
                temp{1,ii}(length(temp{1,7})+1) = [];  %truncate length of all cells to length(temp{1,8})
            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 =  ddd + (hhh + (st1(1,:)+(st1(2,:)+st1(3,:)/1000)/60)/60)/24;
st1(4:7,:) = st1(4:7,:)/100;

%% count missing data points in each 10min section
start = (ddd+hhh/24);
delta = double(10/1440);
last = start + 6*delta;
tbins = start:delta:last;
for ii = 1:6
    Npts = length(find(tsons>=tbins(ii) & tsons<tbins(ii+1)));
    missing(ii) = 12000 - Npts;
end;

%% clean data
U = st1(5,:);  % y axis, positive to N
V = -st1(4,:);  % x axis, positive to W (sonic raw is + to E)
W = st1(6,:);  % z axis, positive upward
Tsonic = st1(7,:);

%% Count NaNs in each 10-min segment
if (all(isfinite(tsons)))
    for ii=1:6;     % get indices in tsons for start/end of each 10-min seg
        cc = find(tsons >= tbins(ii),1,'first');
        dd = find(tsons < tbins(ii+1),1,'last');
        bad(ii) = sum(isnan(U(cc:dd))); % a hack - count Trues (==1) from isnan
    end;
end;

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

%% Interpolate
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);

%% remove spikes in each 10-min segment - 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;

son = [tref' su' sv' sw' st'];
