function [h2]=timeplt(jd,u,istack,ylabels,ylims);
%
% TIMEPLT  Gregorian time series panel plotting
%
%   USAGE: timeplt(jd,u,[istack],[ylabels],[ylims]);
%      where
%         jd = Julian Day time vector  (e.g produced by JULIAN.M)
%         u = column vector or matrix of column vectors containing time
%             series data.  If the column is complex, it will be plotted
%             as a stick plot.  
%         istack = vector of indices indicating which panel you want
%             to plot the time series data.  istack=[1 2] would make
%             two panels on on top of the other and plot the first 
%             column of u in the first panel and the second column of
%             u in the second panel.  If any column in u is complex,  
%             istack must be specified.  If not specified, all the 
%             columns will be plotted in the first panel.
%         ylabels = text matrix containing ylabels.  There should be
%              1 row for each panel.  If there are three panels, then
%              ylabels might look like ['label1';'label2';'label3'].
%              NOTE: all rows must have the same number of characters!
%         ylims  = [npanels x 2] matrix containing the ylimits of 
%              the panel plots.  If you are plotting two panels and 
%              you want the limits of both plots to be from -10 to 15,
%              then set ylims=[-10 15; -10 15].
%

% requires GREGAXD, GREGAXM, GREGAXH, GREGAXY, GREGORIAN, JULIAN
%
% Rich Signell rsignell@usgs.gov
%
% 8-29-94 fixed bug in stick orientation 
% 8-30-94 fixed problem with upper panel too close to top of figure
[m,n]=size(u);
if(m==1),u=u(:);m=n;n=1;end; %if row vector, make column vec
clg
if(nargin==2),
  istack=ones(n,1);
end
njd=length(jd);
jd0=jd(1);
jd1=jd(njd);
ndays=jd1-jd0;
%
% determine how many stack plots to make
%
nstack=max(istack);
% determine height of each stack plot
hstack=.70/nstack;
xmin=.13;
ymin=.15;
width=0.77;
h(1)=axes('position',[xmin ymin width hstack]);
up=u(:,find(istack==1));
if(isreal(up)),
  h2=plot(jd,up);... 
  if(nargin>4);...
    set(gca,'ylim',ylims(1,:));...
  end;...
else
  if(nargin>4),
    ylim=ylims(1,:);
  else
    y0=gmin(imag(up(:)));
    y1=gmax(imag(up(:)));
    ylim=[y0 y1];
  end
  xlim=[jd0 jd1];
  phys=get(gcf,'position');
  uscale=hstack/diff(ylim)*diff(xlim)/width*phys(4)/phys(3);
  vp=imag(up);
  up=uscale.*real(up);
  nsub=1;
  j=[1:nsub:length(up)];
  x=jd;
  xp=x;
  yp=zeros(size(xp));
  xplot=ones(length(xp),2);
  yplot=xplot;
  xplot(:,1)=x(:);
  xplot(:,2)=xp(:)+up(:);
  yplot(:,1)=yp(:);
  yplot(:,2)=yp(:)+vp(:);
  plot([jd0 jd1],[0 0],'w-',xplot(j,:)',yplot(j,:)','g-')
  if(nargin>4),...
      set(gca,'ylim',ylim);...
  end
end
if(nargin>3),ylabel(ylabels(1,:)),end;...
set(gca,'units','pixels');...
pos=get(gca,'pos');xlen=pos(3);ylen=pos(4);...
font=get(gca,'fontsize');...
label_width=5*font;
nlabel=floor(xlen/label_width);
set(gca,'units','normalized');
fac=ndays/nlabel;  %number of days per label
if( fac > 250),	% greater than one year
 % disp('gregaxy')
  gregaxy(jd,floor(fac/250));
elseif(300 >fac & fac > 20),	% one month to one year
  %disp('gregaxm')
  gregaxm(jd,floor(fac/20));
elseif(20 > fac & fac > .2),	% one day to one month
  %disp('gregaxd')
  gregaxd(jd,ceil(fac));		% force one day as min 
elseif(.2 > fac & fac > 0.05); 	% less than one day, .2 = 4.8 h/label
  %disp('gregaxh')
  gregaxh(jd,max(1,floor(fac*48)));
elseif(.05 > fac); 	% less than one day, .05 = 1.2 h/label = 72 min/label
  %disp('gregaxh')
  gregaxmi(jd,max(1,floor(fac*120*24)));
end
if(nstack==1),return,end
pos_norm=get(gca,'position');
vsep=1.3*font*pos_norm(4)/pos(4);
for iplot=2:nstack,
  h(iplot)=axes('position',[xmin ymin+(hstack+vsep)*(iplot-1) width hstack]);
  up=u(:,find(istack==iplot));
if(isreal(up)),
  plot(jd,up);
  if(nargin>4),
    set(gca,'ylim',ylims(iplot,:))
  end
else
  if(nargin>4),
    ylim=ylims(iplot,:);
  else
    y0=gmin(imag(up(:)));
    y1=gmax(imag(up(:)));
    ylim=[y0 y1];
  end
  xlim=[jd0 jd1];
  phys=get(gcf,'position');
  uscale=hstack/diff(ylim)*diff(xlim)/width*phys(4)/phys(3);
  vp=imag(up);
  up=uscale.*real(up);
  nsub=1;
  j=[1:nsub:length(up)];
  x=jd;
  xp=x;
  yp=zeros(size(xp));
  xplot=ones(length(xp),2);
  yplot=xplot;
  xplot(:,1)=x(:);
  xplot(:,2)=xp(:)+up(:);
  yplot(:,1)=yp(:);
  yplot(:,2)=yp(:)+vp(:);
  plot([jd0 jd1],[0 0],'w-',xplot(j,:)',yplot(j,:)','g-')
  if(nargin>4),...
      set(gca,'ylim',ylim);...
  end
end
  set(gca,'xlim',get(h(1),'xlim'),'xtick',get(h(1),'xtick'),'xticklabels',[]);
  if(nargin>3),ylabel(ylabels(iplot,:)),end;...
end
