% 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.
% https://savannah.gnu.org/bugs/?55864
%
% Author: Eddy Xiao.

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
ftol = 10*eps;

b_flip = false;
if lower_bound > upper_bound
  t = upper_bound;
  upper_bound = lower_bound;
  lower_bound = t;
  b_flip = true;
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;  % can be 0.2 or 0.5
  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) <= sqrt(3) * tick_density
  % Major:  1,      2,   3, ..., 9,      10,     20,    30, ...
  % Minor:  1, 1.5, 2, 2.5, ..., 9, 9.5, 10, 15, 20, 25, ...
  ticks_major = fill_log_ticks(lower_bound, upper_bound, 2);
  ticks_minor = fill_log_ticks(lower_bound, upper_bound, 3);
elseif 3 * log10(upper_bound / lower_bound) <= sqrt(3) * tick_density
  % Major:  1, 2,       5,      10, 20,         50, ...
  % Minor:  1, 2, 3, 4, ..., 9, 10, 20, 30, 40, ...
  ticks_major = fill_log_ticks(lower_bound, upper_bound, 1);
  ticks_minor = fill_log_ticks(lower_bound, upper_bound, 2);
elseif log10(upper_bound / lower_bound) <= sqrt(2) * tick_density
  % Major:  1, 10, 100, ...
  % Minor:  1, 2, 3, ... , 9, 10, 20, 30, ...
  ticks_major = fill_log_ticks(lower_bound, upper_bound, 0);
  ticks_minor = fill_log_ticks(lower_bound, upper_bound, 2);
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

if b_flip
  ticks_major = fliplr(ticks_major);
  ticks_minor = fliplr(ticks_minor);
end

end % function

function ticks = fill_log_ticks(lower_bound, upper_bound, density_mode)
%  assert(upper_bound >= lower_bound);

  % density_mode == 0: 1, 10, 100, ...
  % density_mode == 1: 1, 2, 5, 10, 20, 50, ...
  % density_mode == 2: 1, 2, 3, ..., 9, 10, 20, 30, ...
  % density_mode == 3: 1, 1.5, 2, 2.5, ..., 9, 9.5, 10, 15, 20, 25, ...

  switch density_mode
    case 0
      stencil = [1, 10];
    case 1
      stencil = [1, 2, 5, 10];
    case 2
      stencil = 1:10;
    case 3
      stencil = 1:0.5:10;
    otherwise
      error('unsupported density_mode');
  end
  n_stencil = length(stencil) - 1;
  
  lo_lg10 = log10(lower_bound);  % Assume error <= 0.5 ulp,
  hi_lg10 = log10(upper_bound);  % so that floor-ed to correct integer.
  ftol = 4*eps;
  
  magnitude = 10 ^ floor(lo_lg10);
  
  % Find first tick before first 10^n;
  tick_val = lower_bound / magnitude;
  j_stencil = 0;
  while (stencil(j_stencil+1) < tick_val*(1-ftol))
    j_stencil++;
  end
  tick_val = stencil(j_stencil+1) * magnitude;
  
  % locate last tick.
  frac_top = upper_bound / 10 ^ floor(hi_lg10);
  n_top_tick = 0;
  while (stencil(n_top_tick+1) <= frac_top*(1+ftol))
    n_top_tick++;
  end

  n_giant = floor(hi_lg10) - floor(lo_lg10);
  n_ticks = (n_stencil - j_stencil) ...
            + n_stencil * (n_giant - 1) ...
            + n_top_tick;
  ticks = zeros(1, n_ticks);

  idx = 0;
  % fill-in ticks for each order of magnitude
  while magnitude <= upper_bound
    while tick_val < 10*magnitude*(1-ftol) ...
          && tick_val <= upper_bound*(1+ftol)
      ticks(++idx) = tick_val;
      tick_val = stencil(++j_stencil+1) * magnitude;
    end
    magnitude *= 10;
    j_stencil = 0;
  end
  assert(n_ticks == length(ticks));
end
