% Generate major and minor ticks for an axis range [lower_bound, upper_bound]
% tick_density is an expected number of major ticks.
% scaling_mode: 0 <-> linear scale, 1 <-> logarithmic scale.

function [ticks_major, ticks_minor] =...
         calc_axis_ticks(lower_bound, upper_bound, tick_density, scaling_mode)
if isinf(lower_bound+upper_bound) || isnan(lower_bound+upper_bound)
  error("calc_axis_ticks(): value invalid");
end
if lower_bound == upper_bound
  warning("calc_axis_ticks(): lower_bound == upper_bound");
  ticks_major = [lower_bound];
  ticks_minor = [lower_bound];
  return;
end
if scaling_mode == 1 && (lower_bound <=0 || upper_bound <= 0)
  error("calc_axis_ticks(): logscale lower_bound <=0 || upper_bound <= 0 !!");
end

if lower_bound > upper_bound
  t = upper_bound;
  upper_bound = lower_bound;
  lower_bound = t;
end

if scaling_mode == 0 || upper_bound / lower_bound <= 10
  % Linear increment.
%  if tick_density < 5/sqrt(2*5)
%    % Ensure that there are at least 1 major ticks.
%    % The constant here is because: if step_size_desired = sqrt(2*5),
%    % then step_major = 5, this is the most sparse case.
%    % i.e. step size become larger by a factor of 5/sqrt(2*5).
%    % By setting tick_density = 5/sqrt(2*5), we make sure the
%    % resulting minimum n_ticks_major >= 1.
%    tick_density = 5/sqrt(2*5);
%  end
  if scaling_mode == 1
    % reduce tick_density when ratio upper_bound/lower_bound is large
    tick_density /= log(49 * upper_bound/lower_bound + 1) / log(50);
    step_size_desired = (upper_bound/lower_bound) ^ (1/tick_density);
    step_size_desired = (step_size_desired-1) * lower_bound;
  else
    step_size_desired = (upper_bound - lower_bound) / tick_density;
  end
  
  exponent = 10 ^ floor(log10(step_size_desired));
  fraction = step_size_desired / exponent;
  % step size is one of 1, 2, 5, 10, ...
  if fraction < sqrt(1*2)
    step_major = 1;
    step_minor = 0.2;
  elseif fraction < sqrt(2*5)
    step_major = 2;
    step_minor = 1.0;
  elseif fraction < sqrt(5*10)
    step_major = 5;
    step_minor = 1;
  else
    step_major = 10;
    step_minor = 2;
  end
  step_major *= exponent;
  step_minor *= exponent;

  tick_start_major = ceil(lower_bound / step_major) * step_major;
  n_ticks_major = floor((upper_bound - tick_start_major) / step_major) + 1;
  ticks_major = tick_start_major + (0:n_ticks_major-1)*step_major;

  tick_start_minor = ceil(lower_bound / step_minor) * step_minor;
  n_ticks_minor = floor((upper_bound - tick_start_minor) / step_minor) + 1;
  ticks_minor = tick_start_minor + (0:n_ticks_minor-1)*step_minor;
elseif 9 * log10(upper_bound / lower_bound) <= 2 * tick_density
  % Major:  1,      2,   3, ..., 9,      10,     20,    30, ...
  % Minor:  1, 1.5, 2, 2.5, ..., 9, 9.5, 10, 15, 20, 25, ...
  step_major = 10 ^ floor(log10(lower_bound));
  step_minor = step_major / 2;
  tick_major = ceil(lower_bound / step_major) * step_major;
  tick_minor = ceil(lower_bound / step_minor) * step_minor;

  % compute number of tick points
  low_scale  = 10 ^ ceil(log10(lower_bound));
  high_scale = 10 ^ floor(log10(upper_bound));
  n_ticks_giant = floor(log10(upper_bound)) - ceil(log10(lower_bound)) + 1;
  n_ticks_major = floor((low_scale - lower_bound) / step_major) ...
                  + 9 * (n_ticks_giant - 1) ...
                  + floor((upper_bound - high_scale) / high_scale) + 1;
  n_ticks_minor = floor((low_scale - lower_bound) / step_minor) ...
                  + 18 * (n_ticks_giant - 1) ...
                  + floor((upper_bound - high_scale) / (high_scale/2)) + 1;
  ticks_major = zeros(1, n_ticks_major);
  ticks_minor = zeros(1, n_ticks_minor);
  idx_tick = 0;
  idx_val = 0;
  % fill-in ticks for each order of magnitude
  while step_major <= upper_bound
    while tick_major < 10*step_major && tick_major <= upper_bound
      ticks_major(++idx_val) = tick_major;
      tick_major += step_major;
    end
    while tick_minor < 10*step_major && tick_minor <= upper_bound
      ticks_minor(++idx_tick) = tick_minor;
      tick_minor += step_minor;
    end
    step_major *= 10;
    step_minor = step_major / 2;
  end
  %FIXME: one may add an extra tick display level, e.g.
  % Major:  1, 2,       5,      10, 20,         50, ...
  % Minor:  1, 2, 3, 4, ..., 9, 10, 20, 30, 40, ...
elseif log10(upper_bound / lower_bound) <= sqrt(2) * tick_density
  % Major:  1, 10, 100, ...
  % Minor:  1, 2, 3, ... , 9, 10, 20, 30, ...
  step_major = 10 ^ floor(log10(lower_bound));
  tick_minor = ceil(lower_bound / step_major) * step_major;

  % compute number of tick points
  low_scale  = 10 ^ ceil(log10(lower_bound));
  high_scale = 10 ^ floor(log10(upper_bound));
  n_ticks_major = floor(log10(upper_bound)) - ceil(log10(lower_bound)) + 1;
  n_ticks_minor = floor((low_scale - lower_bound) / step_major) ...
                  + 9 * (n_ticks_major - 1) ...
                  + floor((upper_bound - high_scale) / high_scale) + 1;
  ticks_major = zeros(1, n_ticks_major);
  ticks_minor = zeros(1, n_ticks_minor);
  idx_tick = 0;
  idx_val = 0;
  % fill-in ticks for each order of magnitude
  while step_major <= upper_bound
    while tick_minor < 10*step_major && tick_minor <= upper_bound
      ticks_minor(++idx_tick) = tick_minor;
      tick_minor += step_major;
    end
    if step_major >= lower_bound
      ticks_major(++idx_val) = step_major;
    end
    step_major *= 10;
  end
else
  % linear on log scale
  % As long as log10(upper_bound/lower_bound) / tick_density >= sqrt(2),
  % we are sure that ticks_minor are always power of 10.
  [ticks_major, ticks_minor] =...
    calc_axis_ticks(log10(lower_bound), log10(upper_bound), tick_density, 0);
  ticks_major = 10 .^ ticks_major;
  ticks_minor = 10 .^ ticks_minor;
end

