% Make nicer ticks and labels for axies.
% s_ax: Name of axis you want to fix.
%       Can by one of 'x', 'y', 'z', 'c', or combination of them.
%       Repeated letter means that axis boundary need to be mend also.
% tick_density: expected number of major ticks.
%
% Usage example:
% semilogy (logspace (-50, 1, 100))
% fix_axes_ticks ('xxyy');
%
% Author: Eddy Xiao.

function fix_axes_ticks (hax, s_ax, tick_density)
  if ~has_nonempty_var ('hax')
    hax = gca;
  end
  % Called as fix_axes_ticks (s_ax, tick_density).
  if ischar(hax)
    if nargin >= 2
      tick_density = s_ax;
    end
    s_ax = hax;
    hax = gca;
  end
  if ~has_nonempty_var ('s_ax')
    if length (axis ()) == 4
      s_ax = 'xy';   % plot in 2D
    else
      s_ax = 'xyz';  % plot in 3D
    end
  end
  if ~has_nonempty_var ('tick_density')
    tick_density = 4;  % Octave seems using 5.
  end
  if ~isempty(setdiff(s_ax, 'xyzc'))
    error ('invalied character(s) in s_ax');
  end
  
  % Repeated axis name means that axis bound need to be recomputed.
  id_ax_map = struct('x', 1, 'y', 2, 'z', 3, 'c', 4);
  limmodes = 'oooo';  % o: orig, a: auto.
  for ax = s_ax(find (s_ax(2:end) == s_ax(1:end-1)))
    limmodes(id_ax_map.(ax)) = 'a';
  end
  s_ax = unique (s_ax);
  
  % Get data range of all axes (xyzc).
  if any (limmodes ~= 'o')
    orig_lims = axis ();
    axis (hax, 'tight');
    lims = axis ();
    lims
    orig_lims
    axis (hax, orig_lims);
    axis_after = axis ()
    % use tighter bound if possible (e.g. user zoom in)
    if lims(1) < orig_lims(1) || orig_lims(2) < lims(2)
      lims(1:2) = orig_lims(1:2);
    end
    if lims(3) < orig_lims(3) || orig_lims(4) < lims(4)
      lims(3:4) = orig_lims(3:4);
    end
    if length(lims) == 4
      lims = [lims, get(hax,'zlim')];
    end
    if length(lims) == 6
      lims = [lims, get(hax,'clim')];
    end
  end
  
  % Fix the boundary per request and fix the axis ticks.
  for ax = s_ax
    if limmodes(id_ax_map.(ax)) == 'a'
      axlim = lims(2*id_ax_map.(ax) + (-1:0));
      fix_axis_lim (hax, ax, tick_density, axlim);
    end
    fix_axis_ticks (hax, ax, tick_density);
  end
  %get(gca, 'yminortick')  % 'on', 'off'
end

% Nicer boundary by widen the interval axlim up to major ticks.
function axlim = fix_axis_lim (hax, ax, tick_density, axlim)
  scaling_mode = get_scaling_mode (hax, ax);

  % Fix data range, by brute force.
  if scaling_mode == 0
    diff_range = axlim(2) - axlim(1);
    ticksl = calc_axis_ticks (axlim(1) - diff_range, axlim(1),...
                              tick_density, scaling_mode);
    ticksh = calc_axis_ticks (axlim(2), axlim(2) + diff_range,...
                              tick_density, scaling_mode);
    axlim(1) = ticksl(end);
    axlim(2) = ticksh(1);
  else
    diff_range = axlim(2)/axlim(1);
    ticksl = calc_axis_ticks (axlim(1)/diff_range, axlim(1),...
                              tick_density, scaling_mode);
    ticksh = calc_axis_ticks (axlim(2), axlim(2)*diff_range,...
                              tick_density, scaling_mode);
    axlim(1) = ticksl(end);
    axlim(2) = ticksh(1);
  end

  switch ax
    case 'x'
      xlim (hax, axlim);
    case 'y'
      ylim (hax, axlim);
    case 'z'
      zlim (hax, axlim);
    case 'c'
      caxis (hax, axlim);  % not tested
  end
end

% Compute and set ticks/labels in range.
function fix_axis_ticks (hax, ax, tick_density)
  scaling_mode = get_scaling_mode (hax, ax);
  axlim = get (hax, [ax 'lim']);
  fprintf('fix_axis_tics: %s, %g, %g\n', ax, axlim);

  [ticks_major, ticks_minor] =...
    calc_axis_ticks (axlim(1), axlim(2), tick_density, scaling_mode);

  set (hax, [ax 'tick'], ticks_major);
  set (hax, [ax 'minortickvalues'], ticks_minor);
  
  axlim_after = get (hax, [ax 'lim'])
  
  %FIXME fix also labels when suitable (e.g. y is O(1)).
  if all (0.01 < axlim && axlim <= 10000) && scaling_mode == 1
    st_ytk = arrayfun (@(x) sprintf ("%g",x), ticks_major, ...
                       "UniformOutput", false);
    set (hax, [ax 'ticklabel'], st_ytk);
  end
end

function b = has_nonempty_var(var_name)
  b = evalin('caller', ['exist("' var_name '", "var")']) ...
     && ~isempty(evalin('caller', var_name));
end

function scaling_mode = get_scaling_mode (hax, ax)
  if ax == 'c'
    scaling_mode = 0;
  else
    scaling_mode = strcmp (get (hax, [ax 'scale']), 'log');
  end
end
