function timeplt_draw()
% TIMEPLT_DRAW:  Takes care of drawing routines.
%

%disp ( sprintf ( 'here in timeplt_draw' ) );

global timeplt_obj;


N = timeplt_obj_index;


eval ( 'delete ( timeplt_obj{N}.plot_axes );', '' );


%
% Redo all the axes.
num_axes = max(timeplt_obj{N}.istack(:) );

for i = 1:num_axes
    h(i) = subplot(num_axes,1,i);

    if ( ~isempty(timeplt_obj{N}.ylims) )
        set ( h, 'YLim', timeplt_obj{N}.ylims(i,:) );
    end

end

set ( h,     'NextPlot', 'Add' );

            
h = flipud(h(:));





jd = timeplt_obj{N}.jd;
data = timeplt_obj{N}.data;
istack = timeplt_obj{N}.istack;


%
% colordef white color scheme line colors
defaultcmap = get ( 0, 'defaultaxescolororder' );
if ( size(istack,1) > size(defaultcmap,1) )
    cmap = jet(size(istack,1));
else
    cmap = defaultcmap;
end
axis_color_index = zeros(length(istack),1);


	


jd0 = gmin(jd(:));
jd1 = gmax(jd(:));

%
% Make the time axis always extend 5 percent past the
% data.
timediff = jd1 - jd0;
xlim = [(jd0 - timediff*0.05) (jd1 + timediff*0.05)];
set ( h, 'XLim', xlim );


%
% plot the data
for i = 1:length(istack)
    axes ( h(istack(i)) );

	%
	% set the index into the colormap correctly
	% each pane is to be the same as the order for regular matlab plots.
	% While I'm at it, I think I will cure world hunger. 
	axis_color_index(istack(i)) = axis_color_index(istack(i)) + 1;
	
    if ( isreal ( data(:,i) ) )
        graph = plot ( jd, data(:,i) );
        if ( ~isempty(timeplt_obj{N}.ylims) )
            set ( gca, 'YLim', timeplt_obj{N}.ylims(istack(i),:) );
        end
    else

        %
        % Retrieve the ylims.  If ylims were specified
        % in the timplt call, then it has already been
        % set in the gca.  Otherwise construct it from the
        % data.
        if ( isempty(timeplt_obj{N}.ylims) )
            y0 = min ( gmin(imag(data(:,i))), 0);
            y1 = max ( gmax(imag(data(:,i))), 0);
            ylim = [y0 y1];
            set ( gca, 'YLim', ylim );
        else
            ylim = get ( gca, 'ylim' );
        end

        %
        % Construct the scaling factor to accurately do
        % stick plots.
        old_units = get ( gca, 'Units' );
        set ( gca, 'Units', 'Pixels' );
        pos = get ( gca, 'Position' );
        set ( gca, 'Units', old_units );

        scale_factor = (diff(xlim)/diff(ylim)) * (pos(4)/pos(3));

        u_data = scale_factor * real(data(:,i));
        v_data = imag(data(:,i));

        x=jd;
        xp=x;
        yp=zeros(size(xp));
        xplot=ones(length(xp),2);
        yplot=xplot;
        xplot(:,1)=x(:);
        xplot(:,2)=xp(:)+u_data(:);
        xplot(:,3)=x(:);
        yplot(:,1)=yp(:);
        yplot(:,2)=yp(:)+v_data(:);
        yplot(:,3)=yp(:)*nan;
        xplot=xplot';
        yplot=yplot';
        if(~isempty(find(finite(u_data(:)))))
            graph(1) = plot( [jd0 jd1],[0 0] );
            graph(2) = plot ( xplot(:), yplot(:) );...
        end
            
        
    end

%	index = mod(i-1, length(cmap))+1;
%    set ( graph, 'color', cmap(index,:) );
	set ( graph, 'color', cmap(axis_color_index(istack(i)),:) );

end



%
% We want to do our own labels, thank you very much.
set ( h, 'XTickLabel', [] );

%
% What's the separation of the ticks?
xt = get(gca,'XTick' );
xt_sep = max(diff(xt(:)));


axes ( h(1) );
if (timeplt_obj{N}.year_cut_specified)
    timeplt_gregax_year ( timeplt_obj{N}.jd );
elseif (timeplt_obj{N}.month_cut_specified)
    timeplt_gregax_month ( timeplt_obj{N}.jd );
elseif (timeplt_obj{N}.day_cut_specified)
    timeplt_gregax_day ( timeplt_obj{N}.jd );
elseif (timeplt_obj{N}.hour_cut_specified)
    timeplt_gregax_hour( timeplt_obj{N}.jd );
elseif (timeplt_obj{N}.minute_cut_specified)
    timeplt_gregax_minute( timeplt_obj{N}.jd );
elseif ( xt_sep > timeplt_obj{N}.year_cut )
    timeplt_gregax_year ( timeplt_obj{N}.jd );
    timeplt_obj{N}.year_cut_specified = 1;
elseif (xt_sep > timeplt_obj{N}.month_cut) 
    timeplt_gregax_month ( timeplt_obj{N}.jd );
    timeplt_obj{N}.month_cut_specified = 1;
elseif (xt_sep > timeplt_obj{N}.day_cut)
    timeplt_gregax_day ( timeplt_obj{N}.jd );
    timeplt_obj{N}.day_cut_specified = 1;
elseif (xt_sep > timeplt_obj{N}.hour_cut)
    timeplt_gregax_hour( timeplt_obj{N}.jd );
    timeplt_obj{N}.hour_cut_specified = 1;
else
    timeplt_gregax_minute ( timeplt_obj{N}.jd );
    timeplt_obj{N}.minute_cut_specified = 1;
end


%
% Now make sure that all axes have the same tick marks as the first.
set ( h, 'XTick', get(h(1), 'XTick' ) );

    


    
timeplt_obj{N}.plot_axes = h;


%
% The last thing we do is to calculate the printing properties.
% We want the paperposition to come out as close as possible to
% what's on screen.
old_units = get( timeplt_obj{N}.figure, 'Units' );
set ( timeplt_obj{N}.figure, 'Units', 'Points' );
old_paper_units = get ( timeplt_obj{N}.figure, 'PaperUnits' );
set ( timeplt_obj{N}.figure, 'PaperUnits', 'Points' );
set ( timeplt_obj{N}.figure, 'PaperType', 'usletter' );


onscreen_position = get ( timeplt_obj{N}.figure, 'Position' );
onpaper_position = get ( timeplt_obj{N}.figure, 'PaperPosition' );

%
% The maximum width that the plot can be is 612 points.  Anything
% larger must be scaled down.
% The maximum height is 792 points.
if ( onscreen_position(3) > 612 )
    hreduction_factor = 612 / onscreen_position(3);
else
    hreduction_factor = 1;
end

onpaper_position(3) = onscreen_position(3) * hreduction_factor;
onpaper_position(4) = onscreen_position(4) * hreduction_factor;

if ( onscreen_position(4) > 792 )
    vreduction_factor = 792 / onscreen_position(4);
else
    vreduction_factor = 1;
end

onpaper_position(3) = onpaper_position(3) * vreduction_factor;
onpaper_position(4) = onpaper_position(4) * vreduction_factor;

onpaper_position(1) = (612 - onpaper_position(3))/2;
onpaper_position(2) = (792 - onpaper_position(3))/2;

set ( timeplt_obj{N}.figure, 'PaperPosition', onpaper_position );
set ( timeplt_obj{N}.figure, 'Units', old_units );
set ( timeplt_obj{N}.figure, 'PaperUnits', old_paper_units );


return;






function timeplt_gregax_day( jd )
%TIMEPLT_GREGAX_DAY:  Writes gregorian axes tick labels

%disp ( 'present in timeplt_gregax_day' );

month_label = [    'Jan'; ...
                'Feb';...
                'Mar';...
                'Apr';...
                'May';...
                'Jun';...
                'Jul';...
                'Aug';...
                'Sep';...
                'Oct';...
                'Nov';...
                'Dec'];

%
% First figure out the height of the bottom plot in points.
% Then allow for 3 lines to fit here.
old_units = get ( gca, 'Units' );
set ( gca, 'Units', 'Points' );
pos = get(gca, 'Position' );
set ( gca, 'Units', old_units )
height = pos(2);

%
% Fontsize should be the same as the ylabel.
fsize = get ( get(gca,'YLabel'), 'Fontsize' );


%
% Now figure the height of the bottom plot in user coords.
ylim = get ( gca, 'YLim' );
cfactor = diff(ylim) / pos(4);


xt = get ( gca, 'XTick' );
greg = gregorian ( xt(:) );

xlim = get(gca,'xlim');


%
% Make sure that the minute ticks start right on the day,
% and not in the middle.
% Stop the ticks somewhere around the end of xlim.  Don't really
% care so much about this.
jd0 = gmin(jd(:)); jd1 = xlim(2);
start = gregorian(jd0);
start = [start(1:3) 0 0 0];
jd0 = julian(start);
stop = gregorian(jd1);


%
% We need to know how many labels to write.
% Assume that we want a maximum of 50 pts between labels.
max_width = 40;
point_width = pos(3);  
dx = diff(xlim);
scale_factor = dx/point_width;  % number of days per point

day_tic = max ( max_width*scale_factor, 1 );

%
% round this up to the nearest day.
day_tic = ceil(day_tic);

%
% do some more refinement to get good day splits
if ( (day_tic >= 4) & (day_tic <6) )
    day_tic = 5;
elseif ( (day_tic >= 6) & (day_tic <8.5) )
    day_tic = 7;
elseif ( (day_tic >= 8.5) & (day_tic <12) )
    day_tic = 10;
elseif ( (day_tic >= 12) & (day_tic <16) )
    day_tic = 14;
end

xt = [jd0:day_tic:jd1];

ifind = find((xt>=xlim(1)) & (xt<=xlim(2)));
xt = xt(ifind);
set ( gca, ...
        'XTick', xt, ...
        'XTickLabel', [] );

greg = gregorian(xt);





%
% Now label the xticks in days.
gt_string = {     num2str(greg(1,3)); ...
                month_label(greg(1,2),:); ...
                num2str(greg(1,1)) };

ypoint = ylim(1) - 0.5*fsize*cfactor;
xt_label(1) = text ( xt(1), ypoint, 0, gt_string );


%
% Keep track of where to put a label for a new month or year
% by differencing the gregorian dates.  Any entry that is not zero
% means that the entry changed from the previous.
tdiff = diff(greg(:,1:3));

for i = 2:length(xt)

    if ( tdiff(i-1,1)  )
        year_str = num2str(greg(i,1));
    else
        year_str = '';
    end

    if ( tdiff(i-1,2)  )
        month_str = month_label(greg(i,2),:);
    else
        month_str = '';
    end

    if ( tdiff(i-1,3) )
        day_str = num2str(greg(i,3));
    else
        day_str = '';
    end

    gt_string = { day_str; month_str; year_str };

    xt_label(i) = text ( xt(i), ypoint, 0, gt_string );

        
end


set ( xt_label, ...
        'Fontsize', fsize, ...
        'HorizontalAlignment', 'Center', ...
        'VerticalAlignment', 'cap' );





return;










function timeplt_gregax_hour( jd )
%TIMEPLT_GREGAX_HOUR:  Writes gregorian axes tick labels

month_label = [	'Jan'; ...
				'Feb';...
				'Mar';...
				'Apr';...
				'May';...
				'Jun';...
				'Jul';...
				'Aug';...
				'Sep';...
				'Oct';...
				'Nov';...
				'Dec'];

%
% First figure out the height of the bottom plot in points.
% Then allow for 3 lines to fit here.
old_units = get ( gca, 'Units' );
set ( gca, 'Units', 'Points' );
pos = get(gca, 'Position' );
set ( gca, 'Units', old_units )
height = pos(2);

%
% Fontsize should be the same as the ylabel.
fsize = get ( get(gca,'YLabel'), 'Fontsize' );

%
% Fontsize shouldn't be too large.
%fsize = min ( height / 4, 12 );


%
% Now figure the height of the bottom plot in user coords.
ylim = get ( gca, 'YLim' );
cfactor = diff(ylim) / pos(4);


xt = get ( gca, 'XTick' );
greg = gregorian ( xt(:) );

xlim = get(gca,'xlim');


%
% Make sure that the hour ticks start right on the hour,
% and not in the middle.
% Stop the ticks somewhere around the end of xlim.  Don't really
% care so much about this.
% Also, the first tick starts where the data starts.
jd0 = gmin(jd(:)); jd1 = xlim(2);
start = gregorian(jd0);
start = [start(1:4) 0 0];
jd0 = julian(start);
stop = gregorian(jd1);


%
% We need to know how many labels to write.
% Assume that we want a maximum of 50 pts between labels.
max_width = 40;
point_width = pos(3);  
dx = diff(xlim);
scale_factor = dx/point_width;  % points per number of days

hour_tic = max_width * scale_factor;

%
% Now round this up to the nearest minute.
hour_tic = ceil ( hour_tic * 24 ) / 24;

%hour_tic = ceil(max(diff(xt(:)))*24)/24;


xt = [jd0:hour_tic:jd1];

ifind = find((xt>xlim(1)) & (xt<xlim(2)));
xt = xt(ifind);
set ( gca, ...
		'XTick', xt, ...
		'XTickLabel', [] );

greg = gregorian(xt);



%
% Now label the xticks in days.
hour_str = sprintf ( '%02i:00', greg(1,4) );
day_str = sprintf ( '%s %2.0f', month_label(greg(1,2),:), greg(1,3) );
year_str = num2str(greg(1,1));
gt_string = { hour_str; day_str; year_str };

ypoint = ylim(1) - 0.5*fsize*cfactor;
xt_label(1) = text ( xt(1), ypoint, 0, gt_string );


%
% Keep track of where to put a label for a new month or year
% by differencing the gregorian dates.  Any entry that is not zero
% means that the entry changed from the previous.
tdiff = diff(greg(:,1:4));

for i = 2:length(xt)

	if ( tdiff(i-1,4)  )
		hour_str = sprintf ( '%02i:00', greg(i,4) );
	else
		hour_str = '';
	end

	if ( tdiff(i-1,3)  )
		day_str = sprintf ( '%s %2.0f', month_label(greg(i,2),:), greg(i,3) );
	else
		day_str = '';
	end

	if ( tdiff(i-1,3) )
		year_str = num2str(greg(i,1));
	else
		year_str = '';
	end

	gt_string = { hour_str; day_str; year_str };

	xt_label(i) = text ( xt(i), ypoint, 0, gt_string );

		
end


set ( xt_label, ...
		'Fontsize', fsize, ...
		'HorizontalAlignment', 'Center', ...
		'VerticalAlignment', 'cap' );





return;











function timeplt_gregax_minute( jd )
%TIMEPLT_GREGAX_MINUTE:  Writes gregorian axes tick labels

month_label = [	'Jan'; ...
				'Feb';...
				'Mar';...
				'Apr';...
				'May';...
				'Jun';...
				'Jul';...
				'Aug';...
				'Sep';...
				'Oct';...
				'Nov';...
				'Dec'];

%
% First figure out the height of the bottom plot in points.
% Then allow for 3 lines to fit here.
old_units = get ( gca, 'Units' );
set ( gca, 'Units', 'Points' );
pos = get(gca, 'Position' );
set ( gca, 'Units', old_units )
height = pos(2);

%
% Fontsize shouldn't be too large.
fsize = get ( get(gca,'YLabel'), 'Fontsize' );
%fsize = min ( height / 4, 12 );


%
% Now figure the height of the bottom plot in user coords.
ylim = get ( gca, 'YLim' );
cfactor = diff(ylim) / pos(4);


xt = get ( gca, 'XTick' );
greg = gregorian ( xt(:) );

xlim = get(gca,'xlim');


%
% Make sure that the minute ticks start right on the minute,
% and not in the middle.
% Stop the ticks somewhere around the end of xlim.  Don't really
% care so much about this.
jd0 = gmin(jd(:)); jd1 = xlim(2);
start = gregorian(jd0);
start = [start(1:5) 0];
jd0 = julian(start);
stop = gregorian(jd1);


%
% We need to know how many labels to write.
% Assume that we want a maximum of 50 pts between labels.
max_width = 40;
point_width = pos(3);  
dx = diff(xlim);
scale_factor = dx/point_width;  % points per number of days

min_tic = max_width * scale_factor;

%
% Now round this up to the nearest minute.
min_tic = ceil ( min_tic * 24 * 60 ) / 24 / 60;

%min_tic = ceil(max(diff(xt(:)))*24*60)/(24*60);




xt = [jd0:min_tic:jd1];

ifind = find((xt>=xlim(1)) & (xt<=xlim(2)));
xt = xt(ifind);
set ( gca, ...
		'XTick', xt, ...
		'XTickLabel', [] );

greg = gregorian(xt);



%
% Now label the xticks in days.
year_str = num2str(greg(1,1));
day_str = sprintf ( '%s %2.0i', month_label(greg(1,2),:), greg(1,3) ); 
min_str = sprintf ( '%.0f:%02i', greg(1,4), greg(1,5) );
gt_string = { min_str; day_str; year_str };

ypoint = ylim(1) - 0.5*fsize*cfactor;
xt_label(1) = text ( xt(1), ypoint, 0, gt_string );


%
% Keep track of where to put a label for a new month or year
% by differencing the gregorian dates.  Any entry that is not zero
% means that the entry changed from the previous.
tdiff = diff(greg(:,1:5));

for i = 2:length(xt)

	if ( tdiff(i-1,1)  )
		year_str = num2str(greg(i,1));
	else
		year_str = '';
	end

	if ( tdiff(i-1,3)  )
		day_str = sprintf ( '%s %2.0i', month_label(greg(i,2),:), greg(i,3) ); 
	else
		day_str = '';
	end

	if ( tdiff(i-1,5) )
		min_str = sprintf ( '%.0f:%02i', greg(i,4), greg(i,5) );
	else
		min_str = '';
	end

	gt_string = { min_str; day_str; year_str };

	xt_label(i) = text ( xt(i), ypoint, 0, gt_string );

		
end


set ( xt_label, ...
		'Fontsize', round(fsize), ...
		'HorizontalAlignment', 'Center', ...
		'VerticalAlignment', 'cap', ...
		'Tag', 'xtick labels' );




return;












function timeplt_gregax_month( jd )
%TIMEPLT_GREGAX_MONTH:  Writes gregorian axes tick labels

month_label = [	'Jan'; ...
				'Feb';...
				'Mar';...
				'Apr';...
				'May';...
				'Jun';...
				'Jul';...
				'Aug';...
				'Sep';...
				'Oct';...
				'Nov';...
				'Dec'];


%
% First figure out the height of the bottom plot in points.
% Then allow for 3 lines to fit here.
old_units = get ( gca, 'Units' );
set ( gca, 'Units', 'Points' );
pos = get(gca, 'Position' );
set ( gca, 'Units', old_units )
height = pos(2);

%
% Fontsize should be the same as the ylabel.
fsize = get ( get(gca,'YLabel'), 'Fontsize' );



%
% Now figure the height of the bottom plot in user coords.
ylim = get ( gca, 'YLim' );
cfactor = diff(ylim) / pos(4);


xt = get ( gca, 'XTick' );
greg = gregorian ( xt(:) );

xlim = get(gca,'xlim');


%
% Make sure that the minute ticks start right on the day,
% and not in the middle.
% Stop the ticks somewhere around the end of xlim.  Don't really
% care so much about this.
jd0 = gmin(jd(:)); jd1 = xlim(2);
start = gregorian(jd0);
start = [start(1:2) 1 0 0 0];
jd0 = julian(start);
stop = gregorian(jd1);


%
% We need to know how many labels to write.
% Assume that we want a maximum of 50 pts between labels.
max_width = 40;
point_width = pos(3);  
dx = diff(xlim);
scale_factor = dx/point_width;  % points per number of days

mon_tic = ceil ( max_width * scale_factor / 31 );


%mon_tic = ceil(max(diff(xt(:)/31)));

gticks = start;
last_gtick = gticks;

last_jtick = julian(last_gtick);

while ( last_jtick <= xlim(2) )
	this_month = last_gtick(2);
	this_year = last_gtick(1);
	if ( (this_month - 1 + mon_tic) >= 12 )
		this_year = this_year+1;
	end
		
	next_month = mod ( this_month -1 + mon_tic, 12) + 1;

	
	last_gtick = [this_year next_month 1 0 0 0];
	last_jtick = julian(last_gtick);
	gticks = [gticks; last_gtick ];
end




xt = julian(gticks);

ifind = find((xt>=xlim(1)) & (xt<=xlim(2)));
xt = xt(ifind);
set ( gca, ...
		'XTick', xt, ...
		'XTickLabel', [] );

greg = gregorian(xt);



%
% Now figure the height of the bottom plot in user coords.
ylim = get ( gca, 'YLim' );
cfactor = diff(ylim) / pos(4);

%
% Now label the xticks in days.
year_str = num2str(greg(1,1));
month_str = month_label(greg(1,2),:);
gt_string = { month_str; year_str };

ypoint = ylim(1) - 0.5*fsize*cfactor;
xt_label(1) = text ( xt(1), ypoint, 0, gt_string );


%
% Keep track of where to put a label for a new month or year
% by differencing the gregorian dates.  Any entry that is not zero
% means that the entry changed from the previous.
tdiff = diff(greg(:,1:2));

for i = 2:length(xt)

	if ( tdiff(i-1,1)  )
		year_str = num2str(greg(i,1));
	else
		year_str = '';
	end

	if ( tdiff(i-1,2)  )
		month_str = month_label(greg(i,2),:);
	else
		month_str = '';
	end

	gt_string = { month_str; year_str };

	xt_label(i) = text ( xt(i), ypoint, 0, gt_string );

		
end


set ( xt_label, ...
		'Fontsize', fsize, ...
		'HorizontalAlignment', 'Center', ...
		'VerticalAlignment', 'cap' );




return;









function timeplt_gregax_year( jd )
%TIMEPLT_GREGAX_DAY:  Writes gregorian axes tick labels


%
% First figure out the height of the bottom plot in points.
% Then allow for 3 lines to fit here.
old_units = get ( gca, 'Units' );
set ( gca, 'Units', 'Points' );
pos = get(gca, 'Position' );
set ( gca, 'Units', old_units )
height = pos(2);

%
% Fontsize should be the same as the ylabel.
fsize = get ( get(gca,'YLabel'), 'Fontsize' );



%
% Now figure the height of the bottom plot in user coords.
ylim = get ( gca, 'YLim' );
cfactor = diff(ylim) / pos(4);


xt = get ( gca, 'XTick' );
greg = gregorian ( xt(:) );

xlim = get(gca,'xlim');


%
% Make sure that the minute ticks start right on the day,
% and not in the middle.
% Stop the ticks somewhere around the end of xlim.  Don't really
% care so much about this.
jd0 = gmin(jd(:)); jd1 = xlim(2);
start = gregorian(jd0);
start = [start(1) 1  1 0 0 0];
jd0 = julian(start);
stop = gregorian(jd1);


%
% We need to know how many labels to write.
% Assume that we want a maximum of 50 pts between labels.
max_width = 40;
point_width = pos(3);  
dx = diff(xlim);
scale_factor = dx/point_width;  % points per number of days

year_tic = ceil ( max_width * scale_factor / 365.24 );

%year_tic = ceil(max(diff(xt(:)/365.24)));


gyears = start;
last_gyear = start; 
last_jyear = julian(start);
while ( last_jyear <= jd1 )
	last_gyear(1) = last_gyear(1) + year_tic;
	last_jyear = julian(last_gyear);
	gyears = [gyears; last_gyear];
end

xt = julian(gyears);

ifind = find((xt>=xlim(1)) & (xt<=xlim(2)));
xt = xt(ifind);
set ( gca, ...
		'XTick', xt, ...
		'XTickLabel', [] );

greg = gregorian(xt);

year_labels = num2str(greg(:,1));
set ( gca, 'XTickLabel', year_labels );





return;










function xnew=gmin(x)
%function xnew=gmin(x)
% just like min, except that it skips over bad points
[imax,jmax]=size(x);

for j=1:jmax
   good=find(finite(x(:,j)));
   if length(good)>0
       xnew(j)=min(x(good,j));
   else
       xnew(j)=NaN;
   end
end








function xnew=gmax(x)
%function xnew=gmax(x)
% just like max, except that it skips over bad points
[imax,jmax]=size(x);

for j=1:jmax
       good=find(finite(x(:,j)));
       if length(good)>0
          xnew(j)=max(x(good,j));
       else
          xnew(j)=NaN;
       end
end