function label_month(ax_low, ax_up, on_x_axis, reference_day)
% LABEL_MONTH labels an axis of the current figure according to the month
%
% function label_month(ax_low, ax_up, on_x_axis, reference_day)
%
% If there are 3 or less input arguments then it is assumed that the time
% units are years; if there are 4 arguments then the time units are assumed
% to be days and the argument REFERENCE_DAY (see below) is used to specify a
% reference day.
%
%   INPUT ARGUMENTS:
% AX_LOW: is the lower time limit for the plot. If it is a NaN then a 'good'
%         value is chosen. If it is empty the minimum time value is used. If
%         there is no argument passed then the default is 0. 
% AX_UP: is the upper time limit for the plot. If it is a NaN then a 'good'
%         value is chosen. If it is empty the maximum time value is used. If
%         there is no argument passed then the default is 1.
% ON_X_AXIS: if on_x_axis is non-zero then it will be assumed that the x
%            axis is the time axis (this is the default). If on_x_axis is
%            zero then it will be assumed that the y axis is the time axis. 
% REFERENCE_DAY: if this argument is passed then the time units are assumed
%                to be days since a given reference day. The argument may be
%                in a number of forms as listed below:
%             1) REFERENCE_DAY may be a string in one of the date formats
%                0,1,2,6,13,14, 15,16 (as defined by DATESTR). An example is
%                '1-January-1988'.
%             2) REFERENCE_DAY may be the string 'julian' in which case the
%                time units will be assumed to be those used by the CSIRO
%                routines JULIAN and GREGORIAN. (The reference day is that
%                used by astronomers; noon, 24 November, -4713.)
%             3) REFERENCE_DAY may be a number as returned by a matlab call
%                to datenum. e.g. datenum('1-January-1988').
%             4) REFERENCE_DAY may be the "numbers" 0, NaN or []. In this
%                case the time units will be assumed to be those used by the
%                matlab functions DATENUM, DATESTR, DATEVEC, etc. That is the
%                reference day will be 0.

% $Id: label_month.m,v 1.14 2001/07/04 01:53:49 man133 Exp $
% Copyright J. V. Mansbridge, CSIRO, Tue Sep 16 18:01:08 EST 1997

if nargin <= 3
  axis_is_year = 1;
  reference_day = NaN;
  if nargin == 0
    ax_low = 0;
    ax_up = 1;
    on_x_axis = 1;
  elseif nargin == 1
    ax_up = 1;
    on_x_axis = 1;
  elseif nargin == 2
    on_x_axis = 1;
  end
else
  axis_is_year = 0;
  if ischar(reference_day)
    if strcmp(lower(reference_day), 'julian')
      reference_num = datenum(-4713, 11, 24, 12, 0, 0);
    else
      reference_num = datenum(reference_day);
    end
  else
    reference_num = reference_day;
  end
end

% Handle the case of empty arguments

if isempty(ax_low)
  full_month_low = 0;
  ax_low = NaN;
else
  full_month_low = 1;
end
if isempty(ax_up)
  full_month_up = 0;
  ax_up = NaN;
else
  full_month_up = 1;
end
if isempty(on_x_axis) | isnan(on_x_axis)
  on_x_axis = 1;
end
if isempty(reference_day) | isnan(reference_day)
  reference_num = 0;
end

nextpl_gcf = get(gcf, 'NextPlot');
ax = gca;
nextpl_gca = get(ax, 'NextPlot');
set(gcf,'nextplot','add');
set(ax,'nextplot','add');

axis_vals = axis;
mons = ['J' 'F' 'M' 'A' 'M' 'J' 'J' 'A' 'S' 'O' 'N' 'D'];

if on_x_axis ~= 0 % the x axis is labelled
  
  if (isnan(ax_low) | isnan(ax_up))
    hh = get(gca, 'Children');
    xdata = [];
    for ii = 1:length(hh);
      xx = get(hh(ii), 'Type');
      if ((strcmp(xx, 'image')) | (strcmp(xx, 'line')) | ...
	  (strcmp(xx, 'surface')) | (strcmp(xx, 'patch')))
	xx = get(hh(ii), 'XData');
	xdata = [xdata; xx(:)];
      end
    end
  end
  
  if axis_is_year

    % Set the x axes limits
    
    if isnan(ax_low)
      ax_low = axis_vals(1);
    end
    if isnan(ax_up)
      ax_up = axis_vals(2);
    end
    ax_low = round(ax_low*12)/12;
    ax_up = round(ax_up*12)/12;
    if full_month_low
      ax_low_use = ax_low;
    else
      ax_low_use = min(xdata);
    end
    if full_month_up
      ax_up_use = ax_up;
    else
      ax_up_use = max(xdata);
    end
    axis([ax_low_use ax_up_use axis_vals(3) axis_vals(4)])
    xt = ax_low_use:1/12:ax_up_use;
    x_month_mid = (xt(1:end-1) + xt(2:end))/2;
    x_year_mark = ceil(ax_low_use):floor(ax_up_use);
    x_year_mid = (x_year_mark(1:end-1) + x_year_mark(2:end))/2;
    %x_month_mid = (ax_low + 1/24):1/12:ax_up;
    index_month = ceil(mod(x_month_mid, 1)*12);
    %xlabel('month', 'VerticalAlignment', 'top')
  else
    
    % Find the positions of all of the tick marks (i.e., 0 hour on the first
    % day of each month).
    
    if isnan(ax_low)
      x_min = min(xdata);
    else
      x_min = ax_low;
    end
    if isnan(ax_up)
      x_max = max(xdata);
    else
      x_max = ax_up;
    end
    
    [y_st, m_st] = datevec(x_min + reference_num);
    [y_fin, m_fin] = datevec(x_max + reference_num);
    % Begin kluge to work around matlab 6 bug - datevec(0) = [0 0 0 0 0 0]
    if m_st == 0
      m_st = 12;
      y_st = y_st - 1;
    end
    if m_fin == 0
      m_fin = 12;
      y_fin = y_fin - 1;
    end
    % End kluge to work around matlab 6 bug - datevec(0) = [0 0 0 0 0 0]
    if m_st < 12
      year = y_st;
      month = m_st + 1;
    else
      year = y_st + 1;
      month = 1;
    end
    m_fin_p = m_fin + 1;
    if m_fin_p == 13
      y_fin_p = y_fin + 1;
      m_fin_p = 1;
    else
      y_fin_p = y_fin;
    end

    x_begin_month = datenum(y_st, m_st, 1) - reference_num;
    if full_month_low
      ax_low_use = x_begin_month;
    else
      ax_low_use = x_min;
    end
    x_end_month = datenum(y_fin_p, m_fin_p, 1) - reference_num;
    if full_month_up
      ax_up_use = x_end_month;
    else
      ax_up_use = x_max;
    end
    
    axis([ax_low_use ax_up_use axis_vals(3) axis_vals(4)])
    xt = [x_begin_month];
    if month <= 7
      x_year_mark = [datenum(y_st, 1, 1) - reference_num];
    else
      x_year_mark = [];
    end
    index_month = [m_st];
    num_steps = 12*(y_fin - y_st) + m_fin - m_st;
    for ii = 1:num_steps
      xx = datenum(year, month, 1) - reference_num;
      xt = [xt xx];
      index_month = [index_month month];
      month = month + 1;
      if month == 13
	month = 1;
	year = year + 1;
      elseif month == 2
	x_year_mark = [x_year_mark xx];
      end
    end
    xt = [xt x_end_month];
    if month == 1
      x_year_mark = [x_year_mark datenum(year, 1, 1) - reference_num];
    elseif month >= 7
      x_year_mark = [x_year_mark datenum(year+1, 1, 1) - reference_num];
    end
    x_year_mid = (x_year_mark(1:end-1) + x_year_mark(2:end))/2;
    x_month_mid = (xt(1:end-1) + xt(2:end))/2;
  end

  if length(xt) <= 61 % only marks months for 5 years or less
    mark_months = 1;
  else
    mark_months = 0;
  end

  if mark_months
    set(gca, 'XTick', xt)
  else
    set(gca, 'XTick', x_year_mark)
  end
  
  if mark_months | ~axis_is_year
    set(gca, 'XTickLabel', '')
  end
  
  %%%%% Don't write labels outside acceptable range %%%%%%%%%%%%%%%%
  ff = find((x_month_mid < ax_low_use) | (x_month_mid > ax_up_use));
  if ~isempty(ff)
    x_month_mid(ff) = [];
    index_month(ff) = [];
  end 
  ff = find((x_year_mid < ax_low_use) | (x_year_mid > ax_up_use));
  if ~isempty(ff)
    x_year_mid(ff) = [];
  end 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  ydir = get(gca, 'Ydir');
  switch ydir
    case 'normal'
      y_low_1 = axis_vals(3);
      y_low_2 = y_low_1 - 0.05*(axis_vals(4) - axis_vals(3));
    case 'reverse'
      y_low_1 = axis_vals(4);
      y_low_2 = y_low_1 - 0.05*(axis_vals(3) - axis_vals(4));
    otherwise 
      error('weird direction of the y axis')
  end
  if mark_months
    for ii = 1:length(x_month_mid)
      text(x_month_mid(ii), y_low_1, mons(index_month(ii)), ...
	  'HorizontalAlignment', 'Center', ...
	  'VerticalAlignment', 'Top')
    end
  end
  if axis_is_year
    if mark_months
      xlabel('months', 'VerticalAlignment', 'top')
    else
      xlabel('years', 'VerticalAlignment', 'top')
    end
  else
    num_years_print = length(x_year_mid);
    for ii = 1:num_years_print
      timevec = datevec(x_year_mid(ii) + reference_num);
      tt = text(x_year_mid(ii), y_low_2, num2str(timevec(1)), ...
	  'HorizontalAlignment', 'Center', ...
	  'VerticalAlignment', 'Top');
      if num_years_print < 8
	font_size = 15;
      else
	font_size = 105/num_years_print;
      end
      set(tt, 'FontSize', font_size)
    end
  end
  
else % the y axis is labelled
  
  if (isnan(ax_low) | isnan(ax_up))
    hh = get(gca, 'Children');
    ydata = [];
    for ii = 1:length(hh);
      yy = get(hh(ii), 'Type');
      if ((strcmp(yy, 'image')) | (strcmp(yy, 'line')) | ...
	  (strcmp(yy, 'surface')) | (strcmp(yy, 'patch')))
	yy = get(hh(ii), 'YData');
	ydata = [ydata; yy(:)];
      end
    end
  end
  
  if axis_is_year

    % Set the y axes limits
    
    if isnan(ax_low)
      ax_low = axis_vals(3);
    end
    if isnan(ax_up)
      ax_up = axis_vals(4);
    end
    ax_low = round(ax_low*12)/12;
    ax_up = round(ax_up*12)/12;
    if full_month_low
      ax_low_use = ax_low;
    else
      ax_low_use = min(ydata);
    end
    if full_month_up
      ax_up_use = ax_up;
    else
      ax_up_use = max(ydata);
    end
    axis([axis_vals(1) axis_vals(2) ax_low_use ax_up_use])
    xt = ax_low_use:1/12:ax_up_use;
    x_month_mid = (xt(1:end-1) + xt(2:end))/2;
    x_year_mark = ceil(ax_low_use):floor(ax_up_use);
    x_year_mid = (x_year_mark(1:end-1) + x_year_mark(2:end))/2;
    %x_month_mid = (ax_low + 1/24):1/12:ax_up;
    index_month = ceil(mod(x_month_mid, 1)*12);
    %xlabel('month', 'VerticalAlignment', 'top')
  else
    
    % Find the positions of all of the tick marks (i.e., 0 hour on the first
    % day of each month).
    
    if isnan(ax_low)
      y_min = min(ydata);
    else
      y_min = ax_low;
    end
    if isnan(ax_up)
      y_max = max(ydata);
    else
      y_max = ax_up;
    end
    
    [y_st, m_st] = datevec(y_min + reference_num);
    [y_fin, m_fin] = datevec(y_max + reference_num);
    % Begin kluge to work around matlab 6 bug - datevec(0) = [0 0 0 0 0 0]
    if m_st == 0
      m_st = 12;
      y_st = y_st - 1;
    end
    if m_fin == 0
      m_fin = 12;
      y_fin = y_fin - 1;
    end
    % End kluge to work around matlab 6 bug - datevec(0) = [0 0 0 0 0 0]
    if m_st < 12
      year = y_st;
      month = m_st + 1;
    else
      year = y_st + 1;
      month = 1;
    end
    m_fin_p = m_fin + 1;
    if m_fin_p == 13
      y_fin_p = y_fin + 1;
      m_fin_p = 1;
    else
      y_fin_p = y_fin;
    end
    
    x_begin_month = datenum(y_st, m_st, 1) - reference_num;
    if full_month_low
      ax_low_use = x_begin_month;
    else
      ax_low_use = y_min;
    end
    x_end_month = datenum(y_fin_p, m_fin_p, 1) - reference_num;
    if full_month_up
      ax_up_use = x_end_month;
    else
      ax_up_use = y_max;
    end
    
    axis([axis_vals(1) axis_vals(2) ax_low_use ax_up_use])
    xt = [x_begin_month];
    if month <= 7
      x_year_mark = [datenum(y_st, 1, 1) - reference_num];
    else
      x_year_mark = [];
    end
    index_month = [m_st];
    num_steps = 12*(y_fin - y_st) + m_fin - m_st;
    for ii = 1:num_steps
      xx = datenum(year, month, 1) - reference_num;
      xt = [xt xx];
      index_month = [index_month month];
      month = month + 1;
      if month == 13
	month = 1;
	year = year + 1;
      elseif month == 2
	x_year_mark = [x_year_mark xx];
      end
    end
    xt = [xt x_end_month];
    if month == 1
      x_year_mark = [x_year_mark datenum(year, 1, 1) - reference_num];
    elseif month >= 7
      x_year_mark = [x_year_mark datenum(year+1, 1, 1) - reference_num];
    end
    x_year_mid = (x_year_mark(1:end-1) + x_year_mark(2:end))/2;
    x_month_mid = (xt(1:end-1) + xt(2:end))/2;
  end
  
  if length(xt) <= 61 % only marks months for 5 years or less
    mark_months = 1;
  else
    mark_months = 0;
  end

  if mark_months
    set(gca, 'YTick', xt)
  else
    set(gca, 'YTick', x_year_mark)
  end

  if mark_months | ~axis_is_year
    set(gca, 'YTickLabel', '')
  end
  %%%%% Don't write labels outside acceptable range %%%%%%%%%%%%%%%%
  ff = find((x_month_mid < ax_low_use) | (x_month_mid > ax_up_use));
  if ~isempty(ff)
    x_month_mid(ff) = [];
    index_month(ff) = [];
  end 
  ff = find((x_year_mid < ax_low_use) | (x_year_mid > ax_up_use));
  if ~isempty(ff)
    x_year_mid(ff) = [];
  end 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  x_low_1 = axis_vals(1);
  x_low_2 = x_low_1 - 0.05*(axis_vals(2) - axis_vals(1));
  if mark_months
    for ii = 1:length(x_month_mid)
      % text(x_low_1, x_month_mid(ii), mons(index_month(ii)), ...
      % 'HorizontalAlignment', 'Center', ...
      % 'VerticalAlignment', 'Bottom', ...
      % 'Rotation', 90);
      text(x_low_1, x_month_mid(ii), mons(index_month(ii)), ...
	  'HorizontalAlignment', 'Right', ...
	  'VerticalAlignment', 'Middle', ...
	  'Rotation', 0);
    end
  end
  if axis_is_year
    if mark_months
      ylabel('months')
    else
      ylabel('years')
    end
  else
    for ii = 1:length(x_year_mid)
      % tt = text(x_low_2, x_year_mid(ii), num2str(y_st + ii - 1), ...
      %'HorizontalAlignment', 'Center', ...
      % 'VerticalAlignment', 'Bottom', ...
      % 'Rotation', 90);
      timevec = datevec(x_year_mid(ii) + reference_num);
      tt = text(x_low_2, x_year_mid(ii), num2str(timevec(1)), ...
	  'HorizontalAlignment', 'Right', ...
	  'VerticalAlignment', 'Middle', ...
	  'Rotation', 0);
      set(tt, 'FontSize', 15)
    end
  end

end
set(gcf, 'NextPlot', nextpl_gcf);
set(ax, 'NextPlot', nextpl_gca);
