function [son,badTime,badPnts,missing] = read_sonic_WHOTS_2014(dfl,ddd,hhh,sonicmodel,rotationsonic)
%{
reads son file specified by dfl and returns array with the
following columns:

1   jd_ref      10 Hz timestamp
2   U           x axis acceleration, m/s
3   V           y axis acceleration, m/s
4   W           z axis acceleration, m/s
5   Tsonic      x axis rotation rate, rad/s

also returns the number and times for bad data points and the number
of missing data lines in 10-min data segments

input parameters: dfl = string path to scs file
                  ddd = julian date
                  hhh = hour
                  sonicmodel = string parameter
                  rotationsonic = boolean
%}
%% initialize
disp(['Reading son file for hour ',int2str(hhh)]);
flist = fopen(dfl,'r');

st1 = [];               % data array
temp = {};              % cell array for textscan output
badPnts = zeros(1,6);   % initial value for # bad points in each 10-min seg
badTime = zeros(1,6);   % 10-min timestamp
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 ~36000x1 arrays
        temp = textscan(flist,'%2f%2f%3f %*2c %f %f %f %*c %f %*[^\n]','delimiter',', ','headerlines',1,'emptyvalue',NaN,'treatAsEmpty','Sonic');
        st1 = [st1 cell2mat(temp)'];  % 7x~36000 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;

%% count missing data lines
if (all(isfinite(tsons)))
    tdiff = diff(tsons);      % look for data gaps in tsons
    bb = find(tdiff>0.15);    % indices for data gaps
    cc = floor(bb./6000)+1;   % cc identifies which 10-min seg each gap is in 1:6
    missing(cc) = missing(cc) + round(tdiff(bb))*10;    % sum missing lines in each segment
end;                                                    % assuming each sec of gap is 10 lines

%% clean data
u = st1(4,:);       % 1x~36000, m/s
v = st1(5,:);
W = st1(6,:);
Tsonic = st1(7,:);  % 1x~36000, degC

% Count NaNs in each 10-min segment
first = (ddd+hhh/24);
last = first + 50/1440;
delta10 = 1.0/144;
st = (first:delta10:last);
cc = zeros(1,6);
dd = zeros(1,6);
if (all(isfinite(tsons)))
    for bb=1:6;     % get indices in tsons for start/end of each 10-min seg
        cc(bb) = find(tsons > st(bb),1,'first');
        dd(bb) = find(tsons < st(bb)+600,1,'last');
        badPnts(bb) = sum(isnan(u(cc(bb):dd(bb)))); % a hack - count Trues (==1) from isnan
        badTime(bb) = ddd + (hhh+(bb-1)*10/60)/24;  % julian dates for bad segmentss
    end;
end;

[u,N] = despike2(u);         % remove spikes - set rejection criteria in despike2.m
[v,N] = despike2(v);
[W,N] = despike2(W);

% replace any NaNs remain with median
W(isnan(W)) = nanmedian1(W);
u(isnan(u)) = nanmedian1(u);
v(isnan(v)) = nanmedian1(v);

%% Apply rotation as required by particular model and mounting config
%------All result arrays are 1x~36000
if rotationsonic % true when sonic coordinates are rotated by 30 deg
    switch sonicmodel
        case {'R3','WindMasterPro'}     % omnidirectional  (N symbol pointing backward!)
            U = -u*cos(30/180*pi)+v*sin(30/180*pi);     %to North
            V = -u*sin(30/180*pi)-v*cos(30/180*pi);     %to West
        case 'R3A'     % asymmetric      (N symbol pointing forward!)
            U = u*cos(30/180*pi)-v*sin(30/180*pi);      %to North
            V = u*sin(30/180*pi)+v*cos(30/180*pi);      %to West
        case 'R2'      % omnidirectional  (N symbol pointing backward!)
            U = u*cos(30/180*pi)-v*sin(30/180*pi);      %to North
            V = u*sin(30/180*pi)+v*cos(30/180*pi);      %to West
        case 'R2A'     % asymmetric      (N symbol pointing forward!)
            U = -u*cos(30/180*pi)+v*sin(30/180*pi);     %to North
            V = -u*sin(30/180*pi)-v*cos(30/180*pi);     %to West
    end;
else
    switch sonicmodel
        case 'WindMasterPro'	% (N symbol pointing backward!)
            U = -u; %to North
            V = -v; %to West
        case 'R3'               % omnidirectional (N symbol pointing backward!)
            U = -u; %to North
            V = -v; %to West
        case 'R3A'              % asymmetric (N symbol pointing forward!)
            U = u; %to North
            V = v; %to West
    end;
end;

%% Replace unrealistic values with means
% in case despike missed any...
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];
[aa,ii,bb] = 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);

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