function plot_pbo_ts(inp_pos,out_ps,ebars,date_range);

% plot_pbo_ts(inp_pos,out_ps,ebars,date_range);
%
% Given the name of a PBO-standard station position time series file and an
% output filename, plot the North, East, Up time series and put the output in
% a file with the given name.  Optionally add errorbars and select time range.
%
% INPUT
%    inp_pos      Name of PBO-standard station position file to be plotted
%    ebars        Integer.  If 0, no errorbars on plot.  If 1, plot 1-sigma
%                 errorbars.  If 2, plot 2-sigma errorbars
%    date_range   Date range for plot.  Format is [YYYY MM DD yyyy mm dd],
%                 where YYYY is the 4-digit start year, MM is the 2-digit
%                 start month, DD is the 2-digit start day of month, 
%                 where yyyy is the 4-digit ending year, mm is the 2-digit
%                 ending month, and dd is the 2-digit ending day of month.
%
% OUTPUT
%    out_ps       Desired name for output PostScript file.

% error checks

if ( nargin ~= 2 & nargin ~= 3 & nargin ~= 4 )
   error('Error: 2, 3, or 4 input arguments required');
end
if ( nargin == 2 )
   ebars = 0;
   date_range = 0;
end
if ( nargin == 3 )
   date_range = 0;
end
if ( exist(inp_pos) ~= 2 ) 
   str = sprintf('Error: %s not found\n',inp_pos);
   error(str);
end
if ( ebars ~= 0 & ebars ~= 1 & ebars ~= 2 )
   error('Error: ebars must be 0, 1, or 2');
end

% read input

pbo_dot = textread(inp_pos,'%*s %*s %s',1,'headerlines',2);
sta_name = textread(inp_pos,'%*s %*s %*s %s',1,'headerlines',3);
first_epoch = textread(inp_pos,'%*s %*s %*s %d',1,'headerlines',4);
[jday,n,e,u,err_n,err_e,err_u] = textread(inp_pos,'%*d %*d %f %*f %*f %*f %*f %*f %*f %*f %*f %*f %*f %*f %*f %f %f %f %f %f %f %*f %*f %*f %*s','headerlines',9);

% set time scale 

if ( date_range == 0 ) % scale is based on what is in the file
   jday_start = 0;
   jday_end = max(jday)-min(jday)+1;
   jday = jday - min(jday);
else
   start_year = date_range(1);
   start_mon = date_range(2);
   start_day = date_range(3);
   end_year = date_range(4);
   end_mon = date_range(5);
   end_day = date_range(6);
   [trash,jday_start] = ymd2julian(start_year,start_mon,start_day);
   [trash,jday_end]   = ymd2julian(end_year,end_mon,end_day);
   if ( jday_end < jday_start )
      error('Error: end date must be after start date for plotting');
   end
   jday = jday-jday_start;
   jday_end = jday_end-jday_start;
   jday_start = 0;
   first_epoch = date_range(3) + date_range(2)*100 + date_range(1)*1e4;
end

% set values from meters to millimeters

n = n*1000;
e = e*1000;
u = u*1000;
err_n = err_n*1000;
err_e = err_e*1000;
err_u = err_u*1000;

% plot data

subplot(3,1,1)
if ( ebars == 0 )
   plot(jday,n,'bo','MarkerFaceColor','b','MarkerSize',2);
   ax = axis;
   axis([jday_start jday_end ax(3) ax(4)]);
else
   h = errorbar(jday,n,err_n*ebars,'bo');
   set(h(2),'MarkerFaceColor','b','MarkerSize',2);
   ax = axis;
   axis([jday_start jday_end ax(3) ax(4)]);
end

subplot(3,1,2)
if ( ebars == 0)
   plot(jday,e,'bo','MarkerFaceColor','b','MarkerSize',2);
   ax = axis;
   axis([jday_start jday_end ax(3) ax(4)]);
else 
   h = errorbar(jday,e,err_e*ebars,'bo');
   set(h(2),'MarkerFaceColor','b','MarkerSize',2);
   ax = axis;
   axis([jday_start jday_end ax(3) ax(4)]);
end

subplot(3,1,3)
if ( ebars == 0 )
   plot(jday,u,'bo','MarkerFaceColor','b','MarkerSize',2);
   ax = axis;
   axis([jday_start jday_end ax(3) ax(4)]);
else
   h = errorbar(jday,u,err_u*ebars,'bo');
   set(h(2),'MarkerFaceColor','b','MarkerSize',2);
   ax = axis;
   axis([jday_start jday_end ax(3) ax(4)]);
end

subplot(3,1,1)
ylabel 'North (mm)'
title_str = [ char(sta_name) ' (' char(pbo_dot) ')' ];
title(title_str,'Interpreter','None');
subplot(3,1,2)
ylabel 'East (mm)'
subplot(3,1,3)
ylabel 'Up (mm)'
x_str = sprintf('Days since %d',first_epoch);
xlabel(x_str);

% print out results

orient landscape
print(gcf,'-dpsc2',out_ps);



function [jday,mjday] = ymd2julian(year,mon,day);

% mjday = ymd2julian(year,mon,day)
%
% Given the year, mon, and day in the Gregorian calendar, and assuming time
% is 00:00:00, return the Julian day number and modified Julian day number.
%
% INPUT
%   year   4-digit Gregorian year
%   mon    2-digit month
%   day    2-digit day of month
%
% OUTPUT
%
%  jday    Julian day
%  mjday   Modified Julian day
%
% SOURCES
%
% The Calendar FAQ,
% http://www.tondering.dk/claus/cal/node3.html#SECTION003162000000000000000
% 2.16.1 Is there a formula for calculating the Julian day number?

% error checks

if ( nargin ~= 3 )
   error('Error: 3 arguments required');
end
if ( mon < 1 | mon > 12 )
   error('Error: month out of range 1-12');
end
if ( day < 1 ) 
   error('Error: day out of range');
end
if ( (mon==1|mon==3|mon==5|mon==7|mon==8|mon==10|mon==12) & day>31 )
   error('Error: day out of range');
end
if ( (mon==4|mon==6|mon==9|mon==11) & day>30 )
   error('Error: day out of range');
end
if ( isleap(year)==1 & mon==2 & day>29 )
   error('Error: day out of range');
end
if ( isleap(year)==0 & mon==2 & day>28 )
   error('Error: day out of range');
end

% do computation

a = floor((14-mon)/12);
y = year + 4800 - a;
m = mon + 12*a - 3;

t2 = floor((153*m + 2)/5);
t4 = floor(y/4);
t5 = -floor(y/100);
t6 = floor(y/400);

jday = day + t2 + 365*y + t4 + t5 + t6 - 32045;
mjday = jday - 2400000.5;





function [result] = isleap(year);

% result = isleap(year)
%
% Given a Gregorian calendar year from 1 AD forward, return 1 if it is a leap
% year or zero otherwise
%
% INPUT
%   year      4-digit year
%
% OUTPUT
%
%   result    0 if not leap year, 1 if leap year

if ( year < 1 )
  error('Error: year must be 1 or later')
end

if ( (mod(year,4)==0 & mod(year,100)~=0 ) | mod(year,400)==0 )
   result = 1;
else 
   result = 0;
end