function crfp2(year, psfile, ax, RefMatName, SmMatName, min_speed, ... 
		max_speed, max_gap_ratio, cruise, note)
%  1.  NAME
%
%  crfp2 - Plot raw and smoothed velocities
%
%  2.  SYNOPSIS
%
%  crfp2(year, psfile, ax, RefMatName, SmMatName, min_speed, max_speed,
%  max_gap_ratio, cruise, note )
%
%  3.  DESCRIPTION
%
%  This is a modified version of Matlab function callrefp which uses more
%  parameters. It is used to plot raw and smoothed position fixes. It can
%  be called from Matlab script crfp(1m). Parameters min_speed, max_speed
%  and max_gap_ratio should have the same value as those used in smoothr.
%
%  This implementation will also calculates the velocity uncertainty
%  between raw and smoothed velocities, giving an idea of smoothing
%  effects, and therefore data quality.
%
%  4.  PARAMETERS
%
%     year
%        Base year used in input files.
%
%     psfile
%        Name for output postscript file.
%
%     ax Vector with start time,time intervall,start velocity and
%        velocity intervall for plot axes.
%
%     RefMatName
%        Name of output mat file from refabs (parsed through ref2mat).
%
%     SmMatName
%        Name of output mat file from smoothr (parsed through sm2mat).
%
%     min_speed
%        Minimal allowed ship speed for processing good fixes.
%
%     max_speed
%        Maximal allowed ship speed for processing good fixes.
%
%     max_gap_ratio
%        The maximal gap ratio allowed for processing good fixes.
%
%     cruise
%        String to insert in the title.
%
%     note
%        Additionnal string to print on plot.
%
%  5.  AUTHORS
%
%  Pierre Jaccard, Geophysical Institute, University of Bergen, 1999
%

%--------------------------------------------------------------------------
% 
% crfp2.m -- 
% 
% Author          : Jaccard Pierre
% Created On      : 2000/06/29 11:07:42
% Last Modified By: Jaccard Pierre
% Last Modified On: 2001/10/31 09:09:50
% Directory       : /home5/jaccard/pcd/codas/2000-11-21/codas3/gfi/matlab/matlab50/
% ----------------------------------------------------------------------- +
% VERSION CONTROL: 
% 
%   $Author$
%   $Date$
%   $Id$
%   $Name$
%   $Revision$
%   $State$
%   $Log$
% 
%--------------------------------------------------------------------------


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

if (isnan(ref(1,1)) & (isnan(sm(1,1)) | (size(find(sm(:,2) ~=NaN)) == 0)))
  fprintf('No data available.\n')
  return
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;

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);
  end
end

% Here I want to get an idea of the error velocity:
sm_t  = sm(:, j_smtime)';
t_min = max(min(sm(:,j_smtime)), min(tt));
t_max = min(max(sm(:,j_smtime)), max(tt));
t_eps = 0.5 * min( abs(sm_t(2:(length(sm_t))) ...
		- sm_t(1:(length(sm_t)-1))) );
i_sm  = find( (sm_t >= t_min) & (sm_t <= t_max));
err_u = 0;
err_v = 0;
err_n = 0;
for i=i_sm,	
	i_ref = find( (tt > sm_t(i)-t_eps) & (tt < sm_t(i)+t_eps));
	if ~isempty(i_ref) & ~isnan(sm(i,j_sm_u)) & ~isnan(sm(i,j_sm_v)),
		i_ref = i_ref(1);
		err_u = err_u + (sm(i, j_sm_u) - uu(i_ref))^2;
		err_v = err_v + (sm(i, j_sm_v) - vv(i_ref))^2;
		err_n = err_n + 1;
	end
end
if err_n > 0,
	err_u = sqrt(err_u/err_n);
	err_v = sqrt(err_v/err_n);
else
	err_u = -1.0;
	err_v = -1.0;
end

figure
set(gcf, 'PaperType', 'a4');

% plot absolute reference layer velocity vs. time
%subplot(411),
  axes('position',[.2 .73 .65 .18])
  h = plot(tt, uu, '-y', ...
	   notime, null, '+r', ...
	   sm(:,j_smtime), sm(:,j_sm_u), '.g');
  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);
		text(ax(2), ax(4), cruise, 'horizontalalignment', 'right', 'verticalalignment', 'bottom')
	end
  grid
	% Print the error:
	str = 'err = 1.0e38';
	if err_u > 0,
		str = [ 'err = ' sprintf('%-8.3f', err_u) ];
	end
	text(ax(2), ax(3), str, 'rotation', 90, 'FontSize' , 8, ...
			'horizontalalignment', 'left', 'verticalalignment', 'top');
	
%subplot(412)
  axes('position',[.2 .53 .65 .18])
  h = plot(tt, vv, notime, null, '+', sm(:,j_smtime), sm(:,j_sm_v), '.');
  set(h, 'markersize', 5)
  axis(ax);
  set(gca, 'xtick', [ax(1):0.2:ax(2)])
  ylabel('V (m/s)'),
  xlabel([int2str(year), ' Days']),
  grid
	% Print the error:
	str = 'err = 1.0e38';
	if err_v > 0,
		str = [ 'err = ' sprintf('%-8.3f', err_v) ];
	end
	text(ax(2), ax(3), str, 'rotation', 90, 'FontSize' , 8, ...
			'horizontalalignment', 'left', 'verticalalignment', 'top');

% 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)))
  i_goodsm = find(~isnan(sm(:,j_smlon)));
  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
%subplot(413)
  axes('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), '.');   %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), '+'); %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
%subplot(414)
  axes('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), '.');  %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), '+'); %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 -dps -append ' psfile])


%--------------------------------------------------------------------------
% 2001/10/20 PJ
% Add plot of velcity amplitudes

if 0,
figure
orient landscape

subplot(211);
h = plot(tt, curramp(uu,vv), '-b', ...
				 notime, null, '+r', ...
				 sm(:,j_smtime), curramp(sm(:,j_sm_u), sm(:,j_sm_v)), '.r');
set(h, 'markersize', 5)
axis([ax(1), ax(2), 0, 1]);
set(gca, 'xtick', [ax(1):0.2:ax(2)])
ylabel('Speed (m/s)'),
grid on;
title([ 'Absolute Reference Layer Speed - ' cruise ])

subplot(212);
ii = 2:length(tt);
ii = find(tt(ii) ~= tt(ii-1));
ispd = interp1(tt(ii), curramp(uu(ii),vv(ii)), sm(:,j_smtime));
h = plot(sm(:,j_smtime), curramp(sm(:,j_sm_u), sm(:,j_sm_v))-ispd, '-b');
set(h, 'markersize', 5)
axis([ax(1), ax(2), -0.25, +0.25]);
set(gca, 'xtick', [ax(1):0.2:ax(2)])
ylabel('Speed Difference (sm-ref) (m/s)'),
grid on

eval(['print -dps -append ' psfile])


function spd=curramp(u,v)
spd = sqrt(u.*u + v.*v);
return;

end


