function refplot(year, psfile, ax, RefMatName, SmMatName, min_speed, max_speed, max_gap_ratio, cruise)
%
% Usage:  refplot( year, psfile, ax, RefMatName, SmMatName [, cruise] )
% where:  year       = base year for decimal days (x-axis label)
%         psfile     = output PostScript filename
%         ax         = axis parameters [ t0 dt y0 dy ]
%         RefMatName = input Mat-file, unsmoothed (output of ref2mat)
%         SmMatName  = input Mat-file, smoothed (output of sm2mat)
%         cruise     = cruise ID, optional
%
% refplot plots:
% 1) absolute reference layer velocity (smoothed and/or unsmoothed) vs. time
% 2) ship position (longitude and latitude) vs. time
%
% Following variables must be initialized within this function:
%      llmarg -- degrees; plot margin for lon. and lat vs. time
%

llmarg = 0.1;  % degrees; plot boundary for lat-long plot

if (nargin < 8)
  error('Invalid no. of arguments.');
end;

if (length(RefMatName) ~= 0)
  eval(['load ' RefMatName]);
  ref = ref';
end
if (length(SmMatName) ~= 0)
  eval(['load ' SmMatName]);
  sm = sm';
end

[ref_nrow ref_ncol] = size(ref);
[sm_nrow sm_ncol] = size(sm);
if (ref_nrow == 0)
  ref = [NaN NaN NaN NaN NaN NaN NaN NaN NaN];
end
if (sm_nrow == 0)
  sm  = [NaN NaN NaN NaN NaN NaN NaN NaN];
end
% refabs and smoothr column indices:
j_reftime  = 1;
j_lon      = 2;
j_lat      = 3;
j_absu     = 4;
j_absv     = 5;
j_dfixtime = 6;
j_gapu     = 7;
j_gapv     = 8;
j_gaptime  = 9;

j_smtime   = 1;
j_sm_u     = 2;
j_sm_v     = 3;
j_smlon    = 6;
j_smlat    = 7;
if (~isnan(sm(1,1)))
  i_goodsm = find(~isnan(sm(:,j_smlon)));
  if isempty(i_goodsm), sm(1,1) = NaN; end
end


if (isnan(ref(1,1)) & (isnan(sm(1,1)) | (size(find(sm(:,2) ~=NaN)) == 0)))
  fprintf('No data available.\n')
  return
end

tt = NaN;                  % initialize
uu = NaN;
vv = NaN;
notime = NaN;
null = NaN;
% edit the refabs input as in smoothr:
if (ref_nrow > 1)
  minv = ones(size(ref(:, j_reftime))) * min_speed;   % generate a row of min_speed
  maxv = ones(size(ref(:, j_reftime))) * max_speed;   % generate a row of max_speed
  sped = sqrt((ref(:,j_gapu).*ref(:,j_gapu))+(ref(:,j_gapv).* ref(:,j_gapv)));
  speed= max(sped, minv);
  ratio=(speed .* ref(:,j_gaptime) )./ (maxv .* ref(:,j_dfixtime) );
  % ratio = (gap velocity * gaptime) / (maximum velocity * interval of fix );

  % EDITING
  range = find(ratio <  max_gap_ratio);   % find data within max_gap ratio range
  rang1 = find(ratio >= max_gap_ratio);   % find data outside max_gap ratio range
  % assign those data that within max_gap ratio range
  if (length(range) ~= 0)
    abs_u = ref(range, j_absu);     % abs_u
    abs_v = ref(range, j_absv);     % abs_v
    dt    = ref(range, j_dfixtime); % (fix1_time - fix2_time)
    t     = ref(range, j_reftime);  % fix time
  end
  % assign those data that outside max_gap ratio range
  notime = ref(rang1, j_reftime);    %fix time
  null   = notime*0.0 + ax(3) + 0.05*(ax(4)-ax(3));
  if (size(notime) == [0 0])
    notime = NaN;
    null = NaN;
  end
  if (length(range) ~= 0)  % make abs_u and abs_v "bar-shaped"
    n = length(t);
    tt(1:2:(2*n-1)) = t;
    tt(2:2:(2*n-2)) = t(1:n-1) + dt(1:n-1)/(60*24);
    uu(1:2:(2*n-1)) = abs_u;
    uu(2:2:(2*n-2)) = abs_u(1:n-1);
    vv(1:2:(2*n-1)) = abs_v;
    vv(2:2:(2*n-2)) = abs_v(1:n-1);
  else
    fprintf('No good reference layer velocities available.')
    return
  end
end

%clf
%figure('units','normal','position',[.5 .0 .5 1.0])
figure
% plot absolute reference layer velocity vs. time
h1=subplot(411),
  set(h1,'position',[.2 .73 .65 .18])
  h = plot(tt, uu, '.r', ...
           notime, null, '+m', ...
           sm(:,j_smtime), sm(:,j_sm_u), '-b');
  set(h, 'markersize', 5)
  axis(ax);
  set(gca, 'xtick', [ax(1):0.2:ax(2)])
  set(gca, xtl_str, [''])
  ylabel('U (m/s)'),
  title('Absolute Reference Layer Velocity')
  if (nargin == 9)
%    xloc = 0.95 - 0.01444 * length(cruise);
cruise
    text(ax(2), ax(4), cruise, 'horizontalalignment', 'right', 'verticalalignment', 'bottom')
  end
  grid
h2=subplot(412)
  set(h2,'position',[.2 .53 .65 .18])
  h = plot(tt, vv, '.r', ...
            notime, null, '+m', ...
            sm(:,j_smtime), sm(:,j_sm_v), '-b');
  set(h, 'markersize', 5)
  axis(ax);
  set(gca, 'xtick', [ax(1):0.2:ax(2)])
  ylabel('V (m/s)'),
  xlabel([int2str(year), ' Days']),
  grid

% plot the positions (vs. time)
if (~isnan(ref(1,1)))
  if (max(abs(diff(ref(:,j_lon)))) > 180) % ship is crossing the dateline
    neg_ind = find(ref(:,j_lon) < 0);
    ref(neg_ind, j_lon) = ref(neg_ind,j_lon) + 360;
  end;
  lon(1) = min(ref(:,j_lon)) - llmarg;       %longitude lower boundary
  lon(2) = max(ref(:,j_lon)) + llmarg;       %longitude upper boundary
  lat(1) = min(ref(:,j_lat)) - llmarg;       %latitude lower boundary
  lat(2) = max(ref(:,j_lat)) + llmarg;       %latitude upper boundary
end

if (~isnan(sm(1,1)))
  if (max(abs(diff(sm(i_goodsm,j_smlon)))) > 180) %ship is crossing the dateline
    neg_ind = find(sm(i_goodsm,j_smlon) < 0 );
    sm(i_goodsm(neg_ind),j_smlon) = sm(i_goodsm(neg_ind),j_smlon) + 360;
  end;
  lon(1) = min(sm(i_goodsm,j_smlon)) - llmarg;    %longitude lower boundary
  lon(2) = max(sm(i_goodsm,j_smlon)) + llmarg;    %longitude upper boundary
  lat(1) = min(sm(i_goodsm,j_smlat)) - llmarg;    %latitude lower boundary
  lat(2) = max(sm(i_goodsm,j_smlat)) + llmarg;    %latitude upper boundary
end
h3=subplot(413)
  set(h3,'position',[.2 .29 .65 .18])
  dll = max([abs(diff(lat)), abs(diff(lon))]);
  axll = ax;
  axll(3:4) = [lon(1), (lon(1) + dll)];  %only change vertical boundary
  if (~isnan(sm(1,1)))
     h = plot(sm(:,j_smtime), sm(:,j_smlon), '-b');   %lon vs. time smoothed)
     %set(h, 'markersize', 2)
     set(gca, xtl_str, [''])
  end
  set(gca, 'xtick', [ax(1):0.2:ax(2)])
  if (~isnan(ref(1,1)) & ~isnan(sm(1,1)))
    hold on;
  end
  if (~isnan(ref(1,1)))
     h = plot(ref(:,j_reftime), ref(:,j_lon), '.r'); %lon vs. time (un-smoothed)
     set(h, 'markersize', 5)
     set(gca, xtl_str, [''])
  end
  if (~isnan(ref(1,1)) & ~isnan(sm(1,1)))
    hold off;
  end
  axis(axll)
  ylabel('Longitude'),
  grid
h4=subplot(414)
  set(h4,'position',[.2 .09 .65 .18])
  axll(3:4) = [lat(1), (lat(1) + dll)];  %only change vertical boundary
  if (~isnan(sm(1,1)))
     h = plot(sm(:,j_smtime), sm(:,j_smlat), '-b');  %lat vs. time (smoothed)
     %set(h, 'markersize', 2)
  end
  set(gca, 'xtick', [ax(1):0.2:ax(2)])
  if (~isnan(ref(1,1)) & ~isnan(sm(1,1)))
    hold on;
  end
  if (~isnan(ref(1,1)))
     h = plot(ref(:,j_reftime), ref(:,j_lat), '.r'); %lat vs. time (un-smoothed)
     set(h, 'markersize', 5)
  end
  if (~isnan(ref(1,1)) & ~isnan(sm(1,1)))
    hold off;
  end
  axis(axll)
  ylabel('Latitude'),
  xlabel([int2str(year), ' Days']),
  grid
  text(-0.25, -0.35, pctime, 'units', 'normalized', 'rotation', 90.0);    %print PC (current) time on the plot
  orient tall
  eval(['print '  psfile ' -dppm'])
  close all
  