function [att, rows] = att_stats(a, varargin)
% [att, rows] = att_stats(a, varargin)
% calculates the appropriate statistics to use with ashrot
% no indexing: all values are used
%
% a is a structure with at least the following fields from adu2, all Nx1
% aux_struct =
%
%         aduheading: [2315x1 double]
%              pitch: [2315x1 double]
%               roll: [2315x1 double]
%            attdday: [2315x1 double]
%               mrms: [2315x1 double]
%               brms: [2315x1 double]
%         (all reacq fields equal to 1 have been deleted)
%
%% This version is extremely inefficient and unvectorized!


c.max_brms = 0.03;           % 0.03
c.max_mrms = 0.004;           % 0.004
c.max_dh_dev = 1;         % 1 degree
c.max_p_std_dev = 2.5;      % 2.5
c.max_r_std_dev = 2.5;      % 2.5

c = fillstruct(c, varargin);

m_reacq = isnan(a.aduheading) ;

m_brms  = (a.brms > c.max_brms); % NaN > 1 returns false, so
                                 % this mask excludes bad reacq points
m_mrms = (a.mrms > c.max_mrms);

m_good = ~(m_reacq | m_brms | m_mrms);

[m, s, n, med] = mstdgap(a.pitch(m_good));
m_pitch = abs(a.pitch - med) > c.max_p_std_dev * s;
[m, s, n, med] = mstdgap(a.roll(m_good));
m_roll = abs(a.roll - med) > c.max_r_std_dev * s;

dh = a.aduheading - a.heading;

m_dh = abs(dh - bl_filt(dh, 20)) > c.max_dh_dev;
   
m_outlier = m_pitch | m_roll | m_dh;

m_good = ~(~m_good | m_outlier);

dh_mean = mean(dh(m_good));
dh_std  = std(dh(m_good));
dh_min = min(dh(m_good));
dh_max = max(dh(m_good));

p_mean = mean(a.pitch(m_good));
p_std  = std(a.pitch(m_good));
p_min = min(a.pitch(m_good));
p_max = max(a.pitch(m_good));

r_mean = mean(a.roll(m_good));
r_std  = std(a.roll(m_good));
r_min = min(a.roll(m_good));
r_max = max(a.roll(m_good));

brms_mean = mean(a.brms(m_good));
brms_std  = std(a.brms(m_good));
brms_min = min(a.brms(m_good));
brms_max = max(a.brms(m_good));

mrms_mean = mean(a.mrms(m_good));
mrms_std  = std(a.mrms(m_good));
mrms_min = min(a.mrms(m_good));
mrms_max = max(a.mrms(m_good));


rows = {'i_day', 'i_lat', 'i_lon', ...
        'i_dh_mean', 'i_dh_std', 'i_dh_min', 'i_dh_max', ...
        'i_p_mean', 'i_p_std', 'i_p_min', 'i_p_max', ...
        'i_r_mean', 'i_r_std', 'i_r_min', 'i_r_max', ...
        'i_brms_mean', 'i_brms_std', 'i_brms_min', 'i_brms_max', ...
        'i_mrms_mean', 'i_mrms_std', 'i_mrms_min', 'i_mrms_max', ...
        'i_n_att', 'i_n_used', 'i_n_reacq', 'i_n_mrms', 'i_n_brms', ...
        'i_n_outlier'};

att = [a.dday(end), NaN, NaN, ...
       dh_mean, dh_std, dh_min, dh_max, ...
       p_mean, p_std, p_min, p_max, ...
       r_mean, r_std, r_min, r_max, ...
       brms_mean, brms_std, brms_min, brms_max, ...
       mrms_mean, mrms_std, mrms_min, mrms_max, ...
       sum(~m_reacq), sum(m_good), sum(m_reacq), sum(m_mrms), sum(m_brms), ...
       sum(m_outlier)].';

