function sourceval = gm_source(proc, userdata, points)
%  sourceval = gm_source(proc, userdata, points)
% This routine computes source term contributions.  proc is a string
% that defines the function (see the documentation).  points
% is the points to evaluate the function at (zba)  This routine
% is usually invoked by gmfem.
% Returns a zba.

proc = gm_strim(proc);
if proc(1) ~= '('
  error('source-term procedure must start with a parenthesis')
end
if proc(length(proc)) ~= ')'
  error('Unbalanced parentheses in source-term procedure');
end

proc(1) = ' ';
proc(length(proc)) = ' ';
proc = gm_strim(proc);
[first, rest] = strtok(proc);
first = gm_strim(first);
rest = gm_strim(rest);
[n,di] = size(points);
if strcmp(first, 'const')
  sourceval = zba(str2num(rest) * ones(n,1));
elseif strcmp(first, 'vecformula') 
  formula = gm_percentsub(rest, 'double(points(:,', '))');
  sourceval = zba(eval(formula));
elseif strcmp(first, 'formula')
  formula = gm_percentsub(rest, 'double(points(i,', '))');
  sourceval = zba(-ones(n,1));
  for i = 0 : n - 1
    sourceval(i) = eval(formula);
  end
else
  error('Unknown source-term keyword');
end

% ------------------------------------------------------------------
% Copyright (c) 1998 by Cornell University.  All rights reserved.
% See the accompanying file 'Copyright' for authorship information,
% the terms of the license governing this software, and disclaimers
% concerning this software.
% ------------------------------------------------------------------
% This file is part of the QMG software.  
% Version 2.0 of QMG, release date RELDATE
% ------------------------------------------------------------------