% 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 fix 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
%	Times output have 2-digit year at this point, as loadping is not
%	yet compatible with 4-digit year.

% 98/01/28 replaces "pause" with "paused" Jules

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

infile = 'ademo.scn';     % input file 
outfile = 'ademo.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 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

[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;

format long
fprintf(outfile,'abs_max_gap = %.0f  abs_max_bad = %.0f\n',abs_max_gap, abs_max_bad)
message = sprintf('\n Fitting one trend line to ALL data:\n')
prompt = sprintf(' RETURN to continue')
c = polyfit(tt,dt,1);
fit_t = polyval(c,tt);
plot(fix_t,dt,'*',fix_t,fit_t)
xlabel('fix time')
ylabel('PC time - fix time')
title(message)
paused

% 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(3:10),'-',tcd(12:19)];	% reformat with 2-digit year and "-"
ft0d = [t0d(3:10),'-',t0d(12:19)];
fprintf(outfile,message)
fprintf(outfile,' correct_time: %s\n',ftcd)
fprintf(outfile,' PC_time: %s\n', ft0d)
fprintf(outfile,' 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)
message=sprintf('\n %2.0f trends detected as follows: \n',length(ind_gap)+1)
fprintf(outfile,message);
  ind_gap = [0; ind_gap; length(dt)];
  for ii = 2:length(ind_gap)
      N = ind_gap(ii) - ind_gap(ii-1);
      message = sprintf('\n Trend %2.0f (number of points used: %.0f)\n ',ii-1,N)
      fprintf(outfile,message);
      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);
      plot(x,y,'*',x,yfit)
      xlabel('fix time')
      ylabel('PC time - fix time')
      title(message)
      paused
      clock_rate = 1/(c(1)/86400+1);
      xc = x0 - c(2)/86400;
      xcd = to_date(yearbase,xc);
      x0d = to_date(yearbase,x0);
% reformat with 2-digit year and "-"
      fxcd = [xcd(3:10),'-',xcd(12:19)];	
      fx0d = [x0d(3:10),'-',x0d(12:19)];
      fprintf(outfile,' correct_time: %s\n', fxcd)
      fprintf(outfile,' PC_time: %s\n', fx0d)
      fprintf(outfile,' clock rate: %f\n', clock_rate)

  end
end
format short
