function [freq, nfreq] = mfs_addfreq(f, fr, D, nb)

# usage  [freq, nfreq] = mfs_addfreq(f, fr, D, nb)
#
# Input  f(:)    List of excitation frequencies
#        fr(:)   List of resonance frequencies
#        D(:)    List of modal damping ratios
#        nb      Number of frequencies per halfpower bandwidth
#
# Output freq(:) Augmented list of excitation frequencies
#        nfreq   Number of excitation frquencies
#
# The function adds nb equally spaces excitations frequencies in the
# halfpower bandwidth to the list of excitation frequencies.
#
# ------------------------------------------------------------------------

# Resonance frequencies in frequency range

  fmin = min(f);
  fmax = max(f);

  k1 = lookup(fr, fmin) + 1;
  k2 = lookup(fr, fmax);

# Append new excitation frequencies

  for k = k1 : k2
      fu = (1 - D(k)) * fr(k);
      fo = (1 + D(k)) * fr(k);
      f  = [f, linspace(fu, fo, nb)];
  end

# Sort frequencies

  f = sort(f);

# Remove duplicate frequencies and frequencies out of range

  ix1 = lookup(f, fmin, "m");
  ix2 = lookup(f, fmax, "m");

  nfreq   = 1;
  freq(1) = fmin;

  s = 1 + 1e6 * eps;

  for k = ix1 + 1 : ix2
      if (f(k) > s * freq(nfreq))
         freq(++nfreq) = f(k);
      end
  end

end
