% clkrate estimates the PC clockspeed error rate by a least-squares fit
% of the difference between the PC time and satellite time.  It
% automatically weeds out outliers prior to estimation and detects
% breaks in the trend (most probably due to the PC clock being reset
% at mid-cruise), fitting each detected trend separately.
%
% Input:  output of scanping
%         This can be the .trs output of scanping, if Transit fixes were
%         recorded with the ADCP ping files; or the .scn output of
%         scanping, if GPS fixes were recorded with the ADCP ping files.
% Output: .clk file with estimate of PC clockspeed vs. satellite clock,
%         the PC time, and the correct time for each segment, in format
%         expected by loadping.
% NOTE: the times chosen as start of a segment are estimates; they must be
% compared with transitions (such as a change in header within a pingdata file
% or ensemble lenghts longer than normal, which might indicate a reset of the
% PC clock) seen in the scan output file and adjusted as needed.
%
% Usage:  1) Edit the file clkrate.m to specify the input parameters:
%            a) Examine the (PCtime - Fixtime) column of the input file
%          for outliers.
%            b) Set the abs_max_bad parameter below to exclude those outliers.
%            c) Examine the column for breaks in the general trend.
%            d) Set the abs_max_gap parameter below to catch those breaks, if any.
%         2) In Matlab, type 'clkrate'.
%            The script will first attempt to fit all the included points.
%            If breaks are detected, it will then fit each seegment separately and
%            display an estimate of the PC clock rate, starting PC time and
%            correct time for each section.
%         3) Use these for time_correction entries in the loadping control file.
% Note:   It may be necessary to manually edit the scan file
%         and delete specific lines that defy underlying trends
%         but cannot be caught by suitable thresholds.

% 91/05/21 Willa Zhu
% 97/05/12 EF, JF modified to set time when pc=fix time (intercept) to first time
%  in segment for calculating rate in that segment; output now gives for each
%  segment the start PC time and the correct time, as well as the
%  clock_rate, as needed for time correction during load; see NOTE
%  above before using these "as is" in loadping.cnt.
%  PC time and pc-fix now read from scanping output file by function
%  readscn.m
%
% 98/01/28 replaces "pause" with "paused" Jules
% 2001/11/25 EF, various updates and bug fixes;
%    incl. 4-digit yr, which is OK for loadping.
%    Plots are now saved as postscript.
%    Results are appended to output file.
%

% >>>>>>>>>>>>>>>>>>>> user must edit the following lines >>>>>>>>>>>>>>>>>>>>

infile = '9309.scn';     % input file
outfile = '9309.clk';    % output file (diary)
yearbase = 1993;            % base year for decimal days
abs_max_gap = 7;          % break-in-trend threshold (absolute value, seconds)
abs_max_bad = 1000;       % outlier threshold (absolute value, seconds)

% <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< end editing <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
psfile = 'clkrate.ps';
fid_out = fopen(outfile, 'a');
fprintf(fid_out, '\n%% clkrate run on %s\n', datestr(now, 31));
[fix_t, dt] = readscn(infile);

ind_bad = find(abs(dt) >= abs_max_bad);
if(length(ind_bad) ~= 0)
   fix_t(ind_bad) = [];
   dt(ind_bad) = [];
end

t0 = fix_t(1);
tt = fix_t - t0;

title_str = 'Fitting one trend line to ALL data';
c = polyfit(tt,dt,1);
fit_t = polyval(c,tt);
figure
plot(fix_t,dt,'*',fix_t,fit_t)
xlabel('fix time')
ylabel('PC time - fix time')
title(title_str)
date_it('clkrate.m');
print( '-dps', psfile);

% ginput
clock_rate = 1/(c(1)/86400+1);
tc = t0 - c(2)/86400;      % True (fix) time
tcd = to_date(yearbase,tc);
t0d = to_date(yearbase,t0);     % PC time as date

ftcd = [tcd(1:10),'-',tcd(12:19)];  % reformat with  "-"
ft0d = [t0d(1:10),'-',t0d(12:19)];
fprintf(fid_out, ...
      '%%  abs_max_gap = %.0f  abs_max_bad = %.0f\n',abs_max_gap, abs_max_bad);
fprintf(fid_out, ['\n%% ', title_str, '\n']);
fprintf(fid_out, ' correct_time: %s\n',ftcd);
fprintf(fid_out, ' PC_time:      %s\n', ft0d);
fprintf(fid_out, ' clock rate:   %f\n', clock_rate);

% find segments to fit; loop through, calculating rate and correct time

difft = diff(dt);
ind_gap = find(abs(difft) >= abs_max_gap);
if (length(ind_gap) ~= 0)
   fprintf(fid_out, ...
         '\n\% %2.0f trends detected as follows: \n',length(ind_gap)+1);
   ind_gap = [0; ind_gap; length(dt)];
   for ii = 2:length(ind_gap)
      N = ind_gap(ii) - ind_gap(ii-1);
      title_str = sprintf(...
            'Trend %2.0f (number of points used: %.0f)',ii-1,N);
      x0 = fix_t(ind_gap(ii-1)+1);                   % start segment time
      x = fix_t((ind_gap(ii-1)+1):ind_gap(ii));
      xx = x - x0;
      y = dt((ind_gap(ii-1)+1):ind_gap(ii));
      c = polyfit(xx,y,1);
      yfit = polyval(c,xx);
      figure
      plot(x,y,'*',x,yfit)
      xlabel('fix time')
      ylabel('PC time - fix time')
      title(title_str)
      date_it('clkrate.m');
      print( '-dps', 'append', psfile);
      clock_rate = 1/(c(1)/86400+1);
      xc = x0 - c(2)/86400;
      xcd = to_date(yearbase,xc);
      x0d = to_date(yearbase,x0);
      % reformat with "-"
      fxcd = [xcd(1:10),'-',xcd(12:19)];
      fx0d = [x0d(1:10),'-',x0d(12:19)];
      fprintf(fid_out, ['\n%% ', title_str, '\n']);
      fprintf(fid_out,' correct_time: %s\n', fxcd);
      fprintf(fid_out,' PC_time:      %s\n', fx0d);
      fprintf(fid_out,' clock_rate:   %f\n', clock_rate);

   end
end
fprintf(fid_out, '\n');
fclose(fid_out);

