function [ a, f ] = assemble ( node_num, node_xy, nnodes, element_num, ...
  element_node, quad_num, wq, xq, yq, element_area, ib, time )

%% ASSEMBLE assembles the coefficient matrix A and right hand side F.
%
%  Discussion:
%
%    The matrix is known to be banded.  A special matrix storage format
%    is used to reduce the space required.  Details of this format are
%    discussed in the routine DGB_FA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Input, integer NNODES, the number of nodes used to form one element.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Input, integer QUAD_NUM, the number of quadrature points used in assembly.
%
%    Input, real WQ(QUAD_NUM), quadrature weights.
%
%    Input, real XQ(QUAD_NUM,ELEMENT_NUM), YQ(QUAD_NUM,ELEMENT_NUM), the
%    coordinates of the quadrature points in each element.
%
%    Input, real ELEMENT_AREA(ELEMENT_NUM), the area of elements.
%
%    Input, integer IB, the half-bandwidth of the matrix.
%
%    Input, real TIME, the current time.
%
%    Output, real A(3*IB+1,NODE_NUM), the NODE_NUM by NODE_NUM
%    coefficient matrix, stored in a compressed format.
%
%    Output, real F(NODE_NUM), the right hand side.
%
%  Local parameters:
%
%    Local, real BB, BX, BY, the value of some basis function
%    and its first derivatives at a quadrature point.
%
%    Local, real BBB, BBX, BBY, the value of another basis
%    function and its first derivatives at a quadrature point.
%

%
%  Initialize the arrays to zero.
%
  f(1:node_num) = 0.0;
  a(1:3*ib+1,1:node_num) = 0.0;
%
%  The actual values of A and F are determined by summing up
%  contributions from all the elements.
%
  for element = 1 : element_num
%
%  Consider a quadrature point QUAD, with coordinates (X,Y).
%
    for quad = 1 : quad_num

      x = xq(quad,element);
      y = yq(quad,element);
      w = element_area(element) * wq(quad);
%
%  Consider one of the basis functions, which will play the
%  role of test function in the integral.
%
%  We generate an integral for every node associated with an unknown.
%  But if a node is associated with a boundary condition, we do nothing.
%
      for test = 1 : nnodes

        i = element_node(test,element);

        [ bi, dbidx, dbidy ] = qbf ( x, y, element, test, node_xy, ...
          element_node, element_num, nnodes, node_num );

        f(i) = f(i) + w * rhs ( x, y, time ) * bi;
%
%  Consider a basis function used to form the value of the solution function.
%
%  If this basis function is associated with a boundary condition,
%  then subtract the term from the right hand side.
%
%  Otherwise, this term is included in the system matrix.
%
%  Logically, this term goes in entry A(I,J).  Because of the
%  band matrix storage, entry (I,J) is actually stored in
%  A(I-J+2*NHBA+1,J).
%
        for basis = 1 : nnodes

          j = element_node(basis,element);

          [ bj, dbjdx, dbjdy ] = qbf ( x, y, element, basis, node_xy, ...
            element_node, element_num, nnodes, node_num );

          aij = dbidx * dbjdx + dbidy * dbjdy;

          a(i-j+2*ib+1,j) = a(i-j+2*ib+1,j) + w * aij;

        end

      end

    end

  end

  return
end