function [h]=timeplt(jd,u,istack,ylims);
%
% TIMEPLT  time series stack plots with Gregorian Time labels on x axes
%
%   USAGE: [h]=timeplt(jd,u,[istack],[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 one on top of the other and plot the first 
%             column of u in the lower panel and the second column of
%             u in the upper panel.  If any column in u is complex,  
%             istack must be specified.  If istack is not specified, all the 
%             columns will be plotted in the first panel.
%         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].  Autoscales if ylims
%              is not set
%
%       outputs:
%            h = handles for stack plots (axes)

% 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
% 11-1-94 eliminated ylabels from the argument list.  Just use the
%         axes handles returned by timeplt to change axes and then
%         use the standard "ylabel" command.
% 
%-------------------------------------------------------------------
%  Set the cutoff for different types of Gregorian axis types
%  You can adjust these to suit your preferences.  For example, if your
%  plot got labeled with days, but you want hours, increase "daycut" 
%  (until it's larger than the "fac" for your plot).
  
  yearcut=250;
  moncut=20;
  daycut=.2;
%
[m,n]=size(u);
if(m==1),u=u(:);m=n;n=1;end; %if row vector, make column vec
clg
if(~exist('istack')),
  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]);
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
% adjust xfactor for subsequent stretching by Gregorian Date
  if( fac > yearcut),
    xlim=[jd0-180 jd1+180];
  elseif(yearcut > fac & fac > moncut),
    xlim=[jd0-15 jd1+15];
  elseif(moncut > fac & fac > daycut),
    xlim=[floor(jd0)-.5 ceil(jd1)+.5];
  elseif(daycut > fac);
    xlim=[jd0-1/48 jd1+1/48];
  end
up=u(:,find(istack==1));
if(isreal(up)),
  h2=plot(jd,up);... 
  if(exist('ylims')),
    set(gca,'ylim',ylims(1,:));...
  end
else
  if(exist('ylims')),
    ylim=ylims(1,:);
  else
    y0=min(gmin(imag(up(:))),0);
    y1=max(gmax(imag(up(:))),0);
    ylim=[y0 y1];
  end
  set(gca,'units','pixels');
  ppos=get(gca,'position');
  set(gca,'units','norm');
  d=diff(xlim)
  uscale=(diff(xlim)/diff(ylim))*(ppos(4)/ppos(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,:)','w-');...
  set(gca,'ylim',ylim);...
end
set(gca,'xlim',xlim);
if( fac > yearcut),
  gregaxy(jd,floor(fac/yearcut));
elseif(yearcut>fac & fac>moncut),
  gregaxm(jd,floor(fac/moncut));
elseif(moncut>fac & fac>daycut),
  gregaxd(jd,ceil(fac));
elseif(daycut>fac);
  gregaxh(jd,max(1,floor(fac*48)));
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(exist('ylims')),
    set(gca,'ylim',ylims(iplot,:))
  end
else
  if(exist('ylims')),
    ylim=ylims(iplot,:);
  else
    y0=min(0,gmin(imag(up(:))));
    y1=max(0,gmax(imag(up(:))));
    ylim=[y0 y1];
  end
  set(gca,'units','pixels');
  ppos=get(gca,'position');
  set(gca,'units','norm');
  d=diff(xlim)
  uscale=(diff(xlim)/diff(ylim))*(ppos(4)/ppos(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,:)','w-')
  set(gca,'ylim',ylim);...
end
  set(gca,'xlim',xlim,'xtick',get(h(1),'xtick'),'xticklabels',[]);
end
