function [nodes, elements] = mfs_line(nodes, id1, id2, idnew, idelt, ...
                                      type, geom, mat, bias)

# usage: [nodes, elements] = mfs_line(nodes, id1, id2, idnew, idelt, ...
#                                     type, geom, mat, bias)
#
# Input  nodes      Structure array with nodal point data
#        id1        Identifier of first node
#        id2        Identifier of last node
#        idnew      Identifiers of new nodes
#        idelt      Identifiers of elements
#        type       Type of elements to create
#        geom       Structure with geometry data (optional)
#        mat        Structure with material data (optional)
#        bias       Bias (optional)
#
# Output nodes      Structure array with nodal point data
#        elements   Structure array with element data
#
# The function adds nodes along a line between the existing nodes with
# identifiers id1 and id2 to the structure array with nodal point data.
# Subsequently, elements between these nodes are created.
#
# The bias is the ratio between the length of the last interval and
# the first interval. If bias is not given, the nodes are equally spaced.
#
# ------------------------------------------------------------------------

# Copyright (c) 2015 by Johannes Wandinger

# Check arguments

  if (nargin < 6 || nargin > 9 || nargout != 2)
     print_usage();
  endif

  if (! ischar(type))
     error("mfs_line: Element type must be a string\n");
  end

  if (nargin > 6)
     if (! isempty(geom) && ! isstruct(geom))
        error("mfs_line: geom must be a structure\n");
     end
     if (nargin > 7)
        if (! isempty(mat) && ! isstruct(mat))
           error("mfs_line: mat must be a structure\n");
        end
     end
  end

# Determine number of nodes to be created

  nofgen = length(idnew);

# Check number of elements to be created

  nofelt = length(idelt);
  if (nofelt != nofgen + 1)
     error("mfs_line: Number of elements inconsistent with number of nodes\n");
  end

# Create parametric positions

  if (nargin < 9)
     x = linspace(-1, 1, nofgen + 2);
     r = x(2 : nofgen + 1);
  else
     a = 1;
     b = a + (nofgen + 1) * log10(bias) / nofgen;
     x = logspace(a, b, nofgen + 2);
     r = -1 + 2 * (x(2 : nofgen + 1) - x(1)) / (x(nofgen + 2) - x(1));
  end

# Find positions of given nodes

  nofnod = length(nodes);

  ix1 = 1;
  ix2 = 1;

  while(ix1 <= nofnod && nodes(ix1).id != id1)
     ix1++;
  end
  if(ix1 > nofnod)
    error("mfs_linenodes: node %5d does not exist\n", id1);
  end

  while(ix2 <= nofnod && nodes(ix2).id != id2)
     ix2++;
  end
  if(ix2 > nofnod)
    error("mfs_linenodes: node %5d does not exist\n", id2);
  end

# Create nodes

  ix = nofnod + 1;
  for n = 1 : nofgen
      nodes(ix).id = idnew(n);
      nodes(ix).coor  = (1 - r(n)) * nodes(ix1).coor;
      nodes(ix).coor  = 0.5 * (nodes(ix).coor + (1 + r(n)) * nodes(ix2).coor);
      ix++;
  end

# Create elements

  nodid = [id1, idnew, id2];
  for n = 1 : nofelt
      elements(n).id   = idelt(n);
      elements(n).type = type;
      elements(n).nodes = [nodid(n), nodid(n + 1)];
  end

  if (nargin > 6 && ! isempty(geom))
     for n = 1 : nofelt
         elements(n).geom = geom;
     end
  end

  if (nargin > 7 && ! isempty(mat))
     for n = 1 : nofelt
         elements(n).mat = mat;
     end
  end

end
