% ashrot.m
% This script calculates the heading correction <gyro-ashtech>
% angle file for rotating the data. Gaps are filled by linear
% interpolation or extrapolation. Besides the angle file,
% output consists of a plot (optional) of time vs ashtech-gyro
% with filled values plotted as (+), and a log file with
% statistics and gap information.
%
% To run manually, edit the cruise specific variables below.
% It is used by the perl automated processing routine newday.prl,
% so be very careful about changing it. Many parameters are set
% externally if the appropriate variables exist, otherwise they
% are set internally; see "yearbase" below, for example.
%
% -- EF --
% 98/01/28 replaced "pause" with "paused" Jules
%>>>>>>>>>>>>>>>>>>>>> edit the following >>>>>>>>>>>>>>>>>>>>>>>>

ndays = 2.0;         % no. of days to plot per panel
pause_seconds = 0;   % seconds to pause after each figure
no_plot = 0;         % 1 to suppress plotting, 0 otherwise


% Editing criteria:
min_n_used     = 20;
max_dh_span    =  5;
max_dh         = 10;
max_mean_roll  =  7;
max_mean_pitch =  5;


if exist('yearbase') == 0,
        yearbase = 1996;
end
if exist('cruise') == 0,
        cruise = 'MW9603';
end

if exist('dh_offset') == 0
   dh_offset = 0    %
end

if exist('fnstring') == 0,     %file name base.
  fnstring = 'a9603';
end

if exist('ash_mat_prefix') == 0,
   ash_mat_prefix = '';
end

if exist('ash_mat_path') == 0,
   ash_mat_path = '../../nav/';    %keep trailing slash
end

log_file = 'ashrot.log';
ps_file =  [fnstring '.ps'];
out_prefix = '';                   % for angle file; normally ''

%<<<<<<<<<<<<<<<<<<<<<<<<< end editing <<<<<<<<<<<<<<<<<<<<<<<<<<<
ang_file = [out_prefix fnstring '.ang'];

ashfile = [ash_mat_path ash_mat_prefix fnstring '.mat'];

if ~exist(ashfile, 'file')
   error(['\nAshtech attitude file from ubprint,\n   '...
       ashfile '\n   is not found.\n'])
end

load(ashfile, '-mat');

fid_log = fopen(log_file, 'a');
fprintf(fid_log, ['\n--> ashrot ' pctime '\n']);
fprintf(fid_log, ['cruise: ' cruise '\n']);
fprintf(fid_log, ['fnstring: ' fnstring '\n']);
fprintf(fid_log, ['input file:   '  ashfile  '\n']);
fprintf(fid_log, ['output file:  '  ang_file '\n']);
fprintf(fid_log, 'yearbase:         %6d\n', yearbase);
fprintf(fid_log, 'dh_offset:           %6.2f\n', dh_offset);
fprintf(fid_log, 'min_n_used:       %6d\n',      min_n_used);
fprintf(fid_log, 'max_dh_span:        %6.1f\n',  max_dh_span);
fprintf(fid_log, 'max_dh:             %6.1f\n',  max_dh);
fprintf(fid_log, 'max_mean_roll:      %6.1f\n',  max_mean_roll);
fprintf(fid_log, 'max_mean_pitch:     %6.1f\n',  max_mean_pitch);

% Calculate means where data are present, for use in editing.
igood = find(~isnan(att(i_r_mean,:)));

if isempty(igood)
   disp('Warning: No valid attitude measurements were found.')
   good_att_mask = zeros(1, size(att, 2));
else
   dhm = mean(att(i_dh_mean, igood));
   rm = mean(att(i_r_mean, igood));
   pm = mean(att(i_p_mean, igood));

   % Editing mask:
   good_att_mask = att(i_n_used,:) > min_n_used                      ...
                 & att(i_dh_max,:)-att(i_dh_min,:) < max_dh_span     ...
                 & abs(att(i_dh_mean,:) - dhm)     < max_dh          ...
                 & abs(att(i_r_mean,:)  - rm )     < max_mean_roll   ...
                 & abs(att(i_p_mean,:)  - pm )     < max_mean_pitch;
end

i_good_att = find(good_att_mask);

dh = att(i_dh_mean,:).';
tt = att(i_day,:).';
n_att = length(dh);

if isempty(i_good_att)
   disp('Warning: No attitude measurements have passed the editing checks.')
   disp('         A value of 0 will be substituted.')
   dh = zeros(size(tt));
else

   % Fill missing values at start of array with first good value.
   i_0 = i_good_att(1);
   if i_0 ~= 1,
      dh(1:(i_0-1)) = dh(i_0) * ones(i_0-1,1);
   end

   % Fill missing values at end of array with last good value.
   i_n = i_good_att(length(i_good_att));
   if i_n ~= n_att,
      dh((i_n+1):n_att) = dh(i_n) * ones(n_att-i_n, 1);
   end

   % Next: linearly interpolate any gaps in the middle.
   i_bad_att = i_0 - 1 + find(~good_att_mask(i_0:i_n));
   if ~isempty(i_bad_att),
     dh(i_bad_att) = table1([tt(i_good_att) dh(i_good_att)], tt(i_bad_att));
   end
end

%------Save the ascii "angle file" for use by ROTATE.------
day_ang = [tt (-dh + dh_offset)];
eval(['save -ascii ' ang_file ' day_ang'])

%------Calculations are finished; describe the results.-----
if isempty(i_good_att)
   fprintf(fid_log, 'total number of ensembles: %6d\n', n_att);
   fprintf(fid_log, 'Attitude samples prior to editing: %6d\n', length(igood));
   fprintf(fid_log, 'No attitude measurements passed the editing\n');
   fprintf(fid_log, 'Value of 0 was substituted for the Ashtech dh\n');
   fclose(fid_log);
   return
else
   fprintf(fid_log, 'mean angle:         %5.3f\n', mean(day_ang(:,2)));
   fprintf(fid_log, 'std angle:          %5.3f\n', std(day_ang(:,2)));
   fprintf(fid_log, 'total number of ensembles: %6d\n', n_att);
   fprintf(fid_log, 'extrapolated at the start: %6d\n', i_0 - 1);
   fprintf(fid_log, 'extrapolated at the end:   %6d\n', n_att - i_n);
   fprintf(fid_log, 'interpolated:              %6d\n', length(i_bad_att));
end



% Find significant interpolation ranges:
range_mask = [1 good_att_mask];
edge_mask = diff(range_mask);
i_start = find(edge_mask == -1);

range_mask = [good_att_mask 1];
edge_mask = diff(range_mask);
i_end   = find(edge_mask == 1);

n_ranges = length(i_end);
if (n_ranges == length(i_start))
   time_ranges =  [tt(i_start) tt(i_end)];  % two columns
   range_lengths = diff(time_ranges.').';  % one column
   i_long_range = find(range_lengths > 1/(3*24) ); % 20 minutes
   n_long_ranges = length(i_long_range);
   fprintf(fid_log, 'ranges longer than about 25 minutes: %d\n', ...
                        n_long_ranges);
   if n_long_ranges ~= 0
      fprintf(fid_log, 'DD ranges, lengths in hours:\n');
      fprintf(fid_log, '%3d  %9.5f to %9.5f  %5.2f\n',...
         [ (1:n_long_ranges).' ...
            time_ranges(i_long_range,:) ...
            range_lengths(i_long_range)*24].');

      num_strings   = sprintf('%3d  ', 1:n_long_ranges);
      num_strings   = reshape(num_strings, 5, n_long_ranges).';
      start_strings = to_date(yearbase, time_ranges(i_long_range,1));
      end_strings   = to_date(yearbase, time_ranges(i_long_range,2));
      mid_strings   = ' to ';
      mid_strings   = mid_strings(ones(1,n_long_ranges),:);
      lf_strings    = ones(n_long_ranges,1) * '\n';

      fprintf(fid_log, '\nDate ranges:\n');
      fprintf(fid_log, [num_strings start_strings ...
                   mid_strings end_strings lf_strings].');
   end
 else
   fprintf(' ERROR: Number of starts and ends of interp ranges not equal')
end

fprintf(fid_log, '\n');
fclose(fid_log);
%-------------end of calculations and logging-------------

if no_plot == 1, return, end

%----------------------------plot--------------------------

if exist(ps_file)
  eval(['delete ' ps_file]);
end

axleft = [floor(tt(1)):ndays:tt(n_att)];
nplots = length(axleft);
npanels = min(6, nplots); % no. of panels per page
npages = ceil(nplots/npanels);

yrange = 1.2 * [min(dh) max(dh)];
ym = 0.1;       % top & bottom margin
ysep = 0.03;    % separation between panels
yp = (1 - 2 * ym - (npanels-1) * ysep) / npanels;  % panel height
xm = 0.15;                                         % side margin
xp = 1 - 2 * xm;                                   % panel width

for ipage = [1:npages]
  figure
  y0 = 1 - ym - yp;
  for ipanel = [1:npanels]
    if (ipage-1) * npanels + ipanel <= nplots
      axes('position', [xm y0 xp yp])
      plot(tt(i_good_att),     dh(i_good_att)    , '.', ...
        tt(~good_att_mask), dh(~good_att_mask), '+')
      axl = axleft(npanels*(ipage-1) + ipanel);
      axis([axl axl+ndays yrange])
      if ipanel == 1
        title('Ashtech heading - gyro heading')
        text(axl+ndays, yrange(2), cruise, ...
            'horizontalalignment', 'right', 'verticalalignment', 'bottom')
      end
      ylabel('deg')
      y0 = y0 - yp - ysep;
    end
  end
  xlabel([int2str(yearbase) ' decimal days'])
  paused(pause_seconds)
  if (npanels > 2)
    orient('tall')
  else
    orient('landscape')
  end
  date_it('ashrot.m')
  eval(['print -dpsc -append ' ps_file]);
end

