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

data = read_asc(file);
if exist('cal.ps')
  delete cal_hdg.ps
end

% start the calculation -----------------------------------------------

ind = find((abs(data(:,1)) == 1));
t = data(ind,2);
phase = data(ind,11);
acc = data(ind,1);
du = data(ind,5);
dv = data(ind,6);
sign_acc = sign(acc);
t_fit = [];
if(timeOrder >= 1) t_fit = t; end
heading = 90-atan2(dv.*sign_acc, du.*sign_acc)*180/pi;
ilt = find(heading < 0);
if ~isempty(ilt), heading(ilt) = heading(ilt) + 180; end

indErr = 1;
label1 = [label ' Water Track Calibration (with all data points)']; 
orient tall
sp = 1;
clf
while length(indErr) > 0
   if (sp == 1) clf; end
   subplot(2,1,sp)
   [est err] = fitphase([heading phase],label1,t_fit,timeOrder);
   indErr = find(abs(err) > trimErr);
   if(~isempty(indErr))
      heading(indErr) = [];
      phase(indErr) = [];
      t(indErr) = []; 
      if(timeOrder >= 1) t_fit = t; end
      label1 = [label ' Water Track Calibration (delete ' num2str(length(indErr)) ' data points)']; 
   end
   if (sp == 2)
     text(-0.25, -0.35, pctime, 'units', 'normalized', 'rotation', 90.0);
     print -append cal.ps
     paused
   end
   sp = sp + 1;
   if (sp > 2) 
     sp = 1;
   end
end
if (sp == 2)
   text(-0.25, -0.35, pctime, 'units', 'normalized', 'rotation', 90.0);
   print -append cal.ps
   paused
end

if(plot_hdg == 1) 
   clf
   subplot(211)
   plot(t,heading,'*')
   xlabel('Time (days)')
   ylabel('Ship''s Heading');
   title(label);
   grid
end

len = length(at_hdg);
if(len < 1) 
   text(-0.25, -0.35, pctime, 'units', 'normalized', 'rotation', 90.0);
   print -append cal.ps
   paused
   return;
end

at_hdg = sort(at_hdg);
mean_ph = [];
mean_hdg = [];
N = [];
for ii = 1:len
   lower = at_hdg(ii) - del_hdg;
   upper = at_hdg(ii) + del_hdg;
   ind = find(heading >= lower &  heading <= upper);
   if(length(ind) ~= 0)
     ph = phase(ind);
     h = heading(ind);
     mean_ph = [mean_ph; mean(ph)];
     mean_hdg = [mean_hdg; mean(h)];
     N = [N; length(ph)];
   end
end
len = length(N);
if((360+min(at_hdg)-max(at_hdg)) <= del_hdg)
  avg_ph = (N(1)*mean_ph(1)+N(len)*mean_ph(len))/(N(1)+N(len));
  avg_h = (N(1)*(mean_hdg(1)+360)+N(len)*mean_hdg(len))/(N(1)+N(len));
  if(avg_h > 360) avg_h = avg_h - 360; end 
  mean_hdg(1) = avg_h;
  mean_ph(1) = avg_ph;
  mean_hdg(len) = [];
  mean_ph(len) = [];
end
subplot(212)
[est err] = fitphase([mean_hdg mean_ph],label);
text(-0.25, -0.35, pctime, 'units', 'normalized', 'rotation', 90.0);
print -append cal.ps
paused
