function [reactant, product] = ReactionVector(Reaction, A, varargin)

% returns the 'reactant' and 'product' vectors related to 'Reaction'.
% Some notes:
% 1- in time-dependent problems A_old should be used for 'Rankine-Hugoniot'
% BC, and 'A' for Dirichlet BCs.
% 
% 2- each reaction with constant 'reaction rate' has 3 elements:
%     Reaction = {reactant_id, product_id, reaction_rate}
%     example: assume we have 2 species 'a', and 'b', then
%     a ---k_1---> b       is defined as:           Reaction = {1,2,k_1} 
%
% 3- each reaction with a function as 'reaction rate' has 4 elements:
%     Reaction = {reactant_id, product_id, reaction_rate_func, independent_var_ids}
%     example: assume we have 2 species 'a', and 'b', then    
%     a ---1+0.2*b*x----> b  is defined as:           Reaction = {1,2,@(b,x) 1+0.2.*b.*x,[2;3]}

% only product will be changed
if isempty(Reaction{1})
    ReactionRate = Reaction{3}(varargin{2:end});
    product = ReactionRate.*A;
    reactant = 0;
    return 
end

% only reactant will be changed
if isempty(Reaction{2})
    ReactionRate = Reaction{3}(varargin{2:end});
    product = 0;
    reactant = -ReactionRate.*A;    
    return
end

% both reactant and product will be changed
ReactionRate = varargin{1}.*Reaction{3}(varargin{2:end});
product = ReactionRate.*A;
reactant = -ReactionRate.*A;

end  % end of function 'ReactionVector'