function [son,dT_hist,dT_edges,badPnts,missing] = read_sonic_FATIMA_2022(dfl,rotationsonic,read_raw,fsons)
%{
    reads 20 Hz R3A son file specified by dfl and returns 10 Hz table
    with the following variables:
    
    1   dnum        matlab datenum 10 Hz timestamp
    2   U           stern-bow axis wind velocity, m/s
    3   V           stbd-port axis wind velocity, m/s
    4   W           z axis wind velocity, m/s
    5   Tsonic      sonic temperature, C
    
    also returns the number of bad data points (NaNs) and the number
    of missing data lines for each 10-min data segment of raw data
    
    input parameters: dfl: string path to son0 file
                      rotationsonic: boolean 1 or 0, 0 for FATIMA
                      read_raw: boolean flag, 1 = force read of raw txt file
                      fsons = sonic sample output rate setting, Hz
%}

fclose all;

vars = {'dnum','u','v','w','Tsonic'};
units = {'none','m/s','m/s','m/s','degC'};

[fpath,fname,~] = fileparts(dfl);
YY = fname(5:6); yyyy = 2000 + str2double(YY);
DDD = fname(7:9); ddd = str2double(DDD);
HH = fname(10:11);  hh = str2double(HH);

dflm = [fname,'.mat'];
% load son mat table if it exists, unless force read_raw is indicated
if ~read_raw && exist(fullfile(fpath,dflm),'file')
    disp(['Reading Sonic file ',fname]);
    load(fullfile(fpath,dflm),'son','dT_hist','dT_edges','badPnts','missing');
    return
end

%% 10 Hz reference timestamp and averaging bins
start = datenum(yyyy-1,12,31) + ddd + hh/24;
delta = double(1.0/864000);
last = start + 35999*delta;
dnum = start:delta:last;  dnum = dnum(:); % exact 10Hz timestamp
dnum_bin = start:delta:last + double(1.0/864000);
dnum_bin = dnum_bin(:);

%% read data if file exists, if not return NaNs
if exist(dfl,'file')==2
    disp(['Reading Sonic file ',fname]);
    %% read file
    flist = fopen(dfl,'r');
    st1 = [];
    badPnts = zeros(6,1);   % initial value for # bad points in each 10-min seg
    missing = zeros(6,1);   % initial value for # missing lines in each 10-min seg

    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 %*3c %*2c %f %f %f %f %*[^\n]','delimiter',', ',...
                'headerlines',1,'emptyvalue',NaN,'treatAsEmpty','Sonic');
            st1 = [st1 cell2mat(temp)'];
        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);
    st1 = st1';
    tsons =  start + st1(:,1)/1440 + (st1(:,2) + st1(:,3)/1000)/86400;
    % histogram of delta T for the raw data, should all be 0.05 sec
    [dT_hist,dT_edges] = histcounts(diff(tsons)*86400,-0.2:0.025:3);

    %% count missing data lines - gaps in timeseries - and bad data in each segment
    %  negative values mean a few extra points...which is not a problem
    segBins = start:10/1440:start + 1/24; % 10-min bin edges
    expected = fsons*60*10;
    for ii=1:6
        zz = tsons >= segBins(ii) & tsons < segBins(ii+1);
        missing(ii) = fix(expected - sum(zz)); % points missing in each 10min segment
        badPnts(ii) = sum(isnan(st1(zz,4)));   % data NaNs in each 10min segment
    end

    %% Apply 30° Gill rotation if required
    if rotationsonic % true when sonic coordinates are rotated by 30 deg
        utemp = st1(:,4)*cos(30/180*pi) - st1(:,5)*sin(30/180*pi); % to North
        vtemp = st1(:,4)*sin(30/180*pi) + st1(:,5)*cos(30/180*pi); % to West
        st1(:,4) = utemp; clear utemp;
        st1(:,5) = vtemp; clear vtemp;
    end

    %%  average to exactly 10 Hz and despike with 3-sigma limit
    % despike detrended wind and copy NaNs over to original wind vector
    son0 = interval_avg(tsons,st1(:,4:7),dnum_bin);
    temp = despike_3sig(detrend(son0(:,2),'omitnan'));
    U = son0(:,2); ii = isnan(temp); U(ii) = NaN;
    temp = despike_3sig(detrend(son0(:,3),'omitnan'));
    V = son0(:,3); ii = isnan(temp); V(ii) = NaN;
    temp = despike_3sig(detrend(son0(:,4),'omitnan'));
    W = son0(:,4); ii = isnan(temp); W(ii) = NaN;
    temp = despike_3sig(detrend(son0(:,5),'omitnan'));
    Tsonic = son0(:,5); ii = isnan(temp); Tsonic(ii) = NaN;

    % format output table
    temp = [dnum,U,V,W,Tsonic];
    son = array2table(temp,'VariableNames',vars);
    son.Properties.VariableUnits = units;
else
    son0 = NaN(36000,5);
    son0(:,1) = dnum;
    son = array2table(son0,'VariableNames',vars);
    son.Properties.VariableUnits = units;
    badPnts = zeros(6,1);
    missing = ones(6,1)*12000;   % all data missing
    dT_edges = [-0.2:0.025:3-0.025];
    dT_hist = NaN(size(dT_edges));
end

% save as .mat file for quicker loading next time
save(fullfile(fpath,dflm),'son','dT_hist','dT_edges','badPnts','missing');

end

