function varargout = vgeo_sta(yb, fn, x, plt, inst, info)
%------------------------------------------------------------------
% NAME
%
%   vgeo_sta - Plot data from a sta vgeo file
%
% SYNOPSIS
%
%   vgeo_sta(year_base, file_name, x, title, institute, info)
%
% PARAMETERS
%
%
%   year_base
%
%   file_name
%      Name of vgeo sta file, i.e. output of programme vgeo when stat_output
%      option is set. This file has the following columns:
%
%         1. time (decimal day)
%         2. longitude
%         3. latitude
%         4. reference depth
%         5. reference velocity
%         6. error on reference velocity
%         7. number of points used to calculate reference velocity
%
%      Only give the root name, since suffix '.sta' is appended inside the
%      function. 
%
%   x
%      Index to x variable, according to the columns in the data file (i.e.
%      most probably 1, 2 or 3). If negative, the axis direction will be
%      reversed. 
%
%   title
%
%   institute
%
%   info
%      Information string.
%
% DESCRIPTION
%
%   Plot reference velocities with error bars, as well as reference levels and
%   number of points.
%
% OUTPUT
%
%
%   <file_name>_sta.ps
%      Output postscript file.
%
%------------------------------------------------------------------

reverse = 0;
if x < 0,
  reverse = 1;
  x = abs(x);
end
ps_file = [ fn '_sta.ps' ];

%--- Read data:
file = [ fn '.sta' ];
d = read_asc(file);
d = tonan(d);
nc = size(d, 2);
if nc ~= 7,
	error([ 'Wrong number of columns in file ' file ]);
end

%--- Indexes to array:
t = 1; z = 4; v = 5; e = 6; n = 7;
d = d';

%--- Plot data:
figure;
orient tall;

subplot(311);
errorbar(d(x,:), d(v,:), d(e,:), 'o');
ylabel('Reference Velocity');
grid on;
if reverse,
  set(gca, 'XDir', 'reverse');
end

subplot(312);
plot(d(x,:), d(z,:), '+', 'MarkerSize', 5);
ylabel('Reference Depth (m)');
grid on;
if reverse,
  set(gca, 'XDir', 'reverse');
end
set(gca, 'YDir', 'reverse');

subplot(313);
plot(d(x,:), d(n,:), '+', 'MarkerSize', 5);
ylabel('Number of Points');
grid on;
if reverse,
  set(gca, 'XDir', 'reverse');
end

switch(x)
	case 1,
		xlabel([ 'Decimal Days (' int2str(yb) ')' ]);
	case 2,
		xlabel('Longitude');
	case 3,
		xlabel('Latitude');
end

%--- Get name of output file:
%ps_file = newext(next_psf([ fn '.ps' ]);
lgstr = [ '  - PS File: ' ps_file ];
info = [ 'Reference Level - ' info ];

%--- Plot residual things:
subtitle = [ to_date(yb, d(t,1)) '  to  ' ...
			to_date(yb, d(t,length(d(t,:)))) ];
plotstuf('-title', plt, '-subtitle', subtitle, '-note', lgstr, ...
		'-institute', inst, '-info', info);

%--- Print output:
eval([ 'print -dps ' ps_file ]);











