% 93/11/19 JR changed to use harmfit = more generalized version of intfit
%        periods have to be expressed in same unit as t (decimal day)

% 00/05/10 EF added min_inertial and min_good
%    "     JF added check of inertial for plot title 
%    "     JH added parentheses

if ~exist('min_inertial', 'var'), min_inertial = 2.7; end
if ~exist('min_good', 'var'), min_good = 5; end


eval(['load ' infile '_uv.mat'])
eval(['load ' infile '_xy.mat'])
if ~isempty(outfile)
  outfile = [outfile '.ps']
end

if ~append & exist(outfile)
  eval(['delete ' outfile])
end

[m,n] = size(uv);
uvc = uv(:,1:2:n-1) + i * uv(:,2:2:n);
clear uv

% delete the columns which is all NaN, if there are a few NaN in the
% bottom of the column, harmfit will not give NaN.

%igood = find(~isnan(uvc(1,:)));

igood = find(sum(~isnan(uvc)) >= min_good);

uvc = uvc(:,igood);
xyt = xyt(:,igood);
d = - zc;

t = xyt(x_index,:);
if x_index == 1                         % longitude on x-axis
  Xlabel = 'Longitude';
  harmonic = 0;
elseif x_index == 2                     % latitude on x-axis
  Xlabel = 'Latitude';
  harmonic = 0;
elseif x_index == 3                     % time on x-axis
  Xlabel = [num2str(year) ' Days'];
end

% give harmonic reasonable value.
% if number of on station day < 1 day, har = 2
% if number of on station day < min_inertial days, har = 3
% if number of on station day >= min_inertial days, har = 4
if harmonic ~= 0
   osday = t(length(t)) - t(2);
   if (osday  < 1),   harmonic = 2;
   elseif ((osday < min_inertial) & (osday >= 1)),  harmonic = 3;
   elseif (osday >= min_inertial),   harmonic = 4;
   end
end
harmonic
periods = [1 1];
if harmonic > 1, periods = [periods 12.42/24]; end
if harmonic > 2, periods = [periods 1   ]; end
if harmonic > 3, periods = [periods inertial_period]; end
if harmonic < 2, semidiurnal = 0; end
if harmonic < 3, diurnal     = 0; end
if harmonic < 4, inertial    = 0; end

if mean_trend | residual | inertial | diurnal | semidiurnal | harmonic
  [b,x,r,f,v] = harmfit(uvc, t, periods);
  eval([' save ' infile ' b x r f v periods t ']);
end

uscale = vscale * (1/72) * (376 / abs(xyrange(2)-xyrange(1)));
nplots = 0;
clf
orient tall
if velocity == 1
  subplot(2,1,1+rem(nplots, 2))
  tdvec(uvc, 1, t, d, xyrange, uscale, key, arrowlength);
  if x_index == 3
    title(['Velocity On Station' plot_title]);
  else
    title(plot_title)
  end
  xlabel(Xlabel);
  ylabel('Depth (m)')
  ax=axis;
  text(ax(2),ax(4),cruise,'HorizontalAlignment','right','VerticalAlignment','bottom');
  nplots = nplots + 1;
end

if mean_trend == 1
  subplot(2,1,1+rem(nplots,2))
  nplots = nplots + 1;
  tdvec(f, 1, t, d, xyrange, uscale, key, arrowlength);
  if inertial == 0
    title(['Mean, Trend, Tides' plot_title])
  else
    title(['Mean, Trend, Tides, Inertial' plot_title])
  end
  xlabel(Xlabel)
  ylabel('Depth (m)')
  if rem(nplots,2) == 0
    if ~isempty(outfile) fprintf(['Writing ' outfile '...']); eval(['print -append ' outfile]); disp('done.'); end
    if pause_flag , disp('paused');pause; end
    clf
  else
    ax=axis;
    text(ax(2),ax(4),cruise,'HorizontalAlignment','right','VerticalAlignment','bottom');
  end
end

if residual == 1
  subplot(2,1,1+rem(nplots,2))
  tdvec(r, 1, t, d, xyrange, uscale, key, arrowlength);
  nplots = nplots + 1;
  title(['Residuals' plot_title])
  xlabel(Xlabel)
  ylabel('Depth (m)')
  if ~isempty(outfile) & rem(nplots,2) == 0
    if ~isempty(outfile) fprintf(['Writing ' outfile '...']); eval(['print -append ' outfile]); disp('done.'); end
    if pause_flag , disp('paused');pause; end
    clf
  else
    ax=axis;
    text(ax(2),ax(4),cruise,'HorizontalAlignment','right','VerticalAlignment','bottom');
  end
end

if inertial == 1
  subplot(2,1,1+rem(nplots,2))
  tdvec(x(:,7:8), b(7:8,:), t, d, xyrange, uscale, key, arrowlength);
  nplots = nplots + 1;
  title(['Inertial cw and ccw' plot_title])
  xlabel(Xlabel)
  ylabel('Depth (m)')
  if rem(nplots,2) == 0
    if ~isempty(outfile) fprintf(['Writing ' outfile '...']); eval(['print -append ' outfile]); disp('done.'); end
    if pause_flag , disp('paused');pause; end
    clf
  else
    ax=axis;
    text(ax(2),ax(4),cruise,'HorizontalAlignment','right','VerticalAlignment','bottom');
  end
end

if diurnal == 1
  subplot(2,1,1+rem(nplots,2))
  tdvec(x(:,5:6), b(5:6,:), t, d, xyrange, uscale, key, arrowlength);
  nplots = nplots + 1;
  title(['Diurnal cw and ccw' plot_title])
  xlabel(Xlabel)
  ylabel('Depth (m)')
  if rem(nplots,2) == 0
    if ~isempty(outfile) fprintf(['Writing ' outfile '...']); eval(['print -append ' outfile]); disp('done.'); end
    if pause_flag , disp('paused');pause; end
    clf
  else
    ax=axis;
    text(ax(2),ax(4),cruise,'HorizontalAlignment','right','VerticalAlignment','bottom');
  end
end

if semidiurnal == 1
  subplot(2,1,1+rem(nplots,2))
  tdvec(x(:,3:4), b(3:4,:), t, d, xyrange, uscale, key, arrowlength);
  nplots = nplots + 1;
  title(['Semidiurnal cw and ccw' plot_title])
  xlabel(Xlabel)
  ylabel('Depth (m)')
  if rem(nplots,2) == 0
    if ~isempty(outfile) fprintf(['Writing ' outfile '...']); eval(['print -append ' outfile]); disp('done.'); end
    if pause_flag , disp('paused');pause; end
    clf
  else
    ax=axis;
    text(ax(2),ax(4),cruise,'HorizontalAlignment','right','VerticalAlignment','bottom');
  end
end

if harmonic > 0
  xyrange2 = [0 4 xyrange(3) xyrange(4)];
  oldxrange = abs(xyrange(2) - xyrange(1));
  key2 = [4*(key(1)-xyrange(1))/oldxrange key(2) key(3)];
  rescale = 1 / (4 / oldxrange);
  uscale2 = uscale * rescale;
  subplot(2,1,1+rem(nplots,2))
  step = xyrange2(1)+[1:7]*(xyrange2(2)-xyrange2(1))/8;
  h = tdvec(x(:,1),b(1,:), t, d, xyrange2, uscale2, key2, arrowlength, [step(1) 1]);
  set(gca, xtl_str, [''])
  nplots = nplots + 1;
  hold on
  tdvec(x(:,1:2),b(1:2,:), t, d, xyrange2, uscale2, [step(1) 1]);
  title(['Harmonic Analysis of Velocity' plot_title])
  ylabel('Depth (m)')
  text(0.5,xyrange2(3)/0.97, 'mean', 'HorizontalAlignment', 'center')        % x coordinate
  text(0.5,xyrange2(3)/0.93, '+ trend', 'HorizontalAlignment', 'center')     % x coordinate
  grid

  if harmonic > 1
    tdvec(x(:,3:4),b(3:4,:), t, d, xyrange2, uscale2, [step(2) 2]);
    tdvec(x(:,3:4),b(3:4,:), t, d, xyrange2, uscale2, arrowlength, [step(3) 3]);
    xx = 1 - xyrange2(2)*(1-0.31)/4;
    text(1.25,xyrange2(3)/0.97, 'semidiurnal', 'HorizontalAlignment', 'center')     % x coordinate
    text(1.25,xyrange2(3)/0.93, '12.42 hours', 'HorizontalAlignment', 'center')     % x coordinate

    if harmonic > 2
      tdvec(x(:,5:6),b(5:6,:), t, d, xyrange2, uscale2, [step(4) 2]);
      tdvec(x(:,5:6),b(5:6,:), t, d, xyrange2, uscale2, arrowlength, [step(5) 3]);
      text(2.25,xyrange(3)/0.97, 'diurnal', 'HorizontalAlignment', 'center')        % x coordinate
      text(2.25,xyrange(3)/0.93, '24 hours', 'HorizontalAlignment', 'center')       % x coordinate

      if harmonic > 3
        tdvec(x(:,7:8),b(7:8,:), t, d, xyrange2, uscale2, [step(6) 2]);
        tdvec(x(:,7:8),b(7:8,:), t, d, xyrange2, uscale2, arrowlength, [step(7) 3]);
        text(3.25,xyrange(3)/0.97, 'inertial', 'HorizontalAlignment', 'center')     % x coordinate
        text(3.25,xyrange(3)/0.93, sprintf('%0.0f hours', inertial_period*24), 'HorizontalAlignment', 'center')
      end
    end
  end
  %set(h, 'markersize', 1)
  if rem(nplots,2) == 0
    if ~isempty(outfile) fprintf(['Writing ' outfile '...']); eval(['print -append ' outfile]); disp('done.'); end
    if pause_flag , disp('paused');pause; end
  else
    ax=axis;
    text(ax(2),ax(4),cruise,'HorizontalAlignment','right','VerticalAlignment','bottom');
  end
  hold off
end

if ~isempty(outfile) & rem(nplots,2) ~= 0
  fprintf(['Writing ' outfile '...']); eval(['print -append ' outfile]); disp('done.');
end
