%
% add_reaction
%
function add_reaction(reactants, products, rate)

% Variables
global reaction_reactants;
global reaction_reactants_index;
global reaction_reactants_num;
global reaction_products;
global reaction_rates;
global num_reactions;
global num_external_input_reactions;
global external_input_reactions;
global reaction_has_external_inputs;

num_reactions = num_reactions + 1;
reaction_has_external_inputs(num_reactions) = false;  % may change

disp(['Reaction ' num2str(num_reactions) ':' reactants '->' products]);

% populate reactant names and their indices
reaction_reactants(num_reactions,1) = {reactants};
reactant_indices = [];
for i=1:size(reactants,2)
    ri = lookup_species(reactants(i));
    if (i==1)
        reactant_indices = [ri];
    else
        reactant_indices = [reactant_indices ri];
    end
    % store information about all reactions where any reactant is external
    if (ri < 0)
        num_external_input_reactions = num_external_input_reactions + 1;
        external_input_reactions(num_external_input_reactions) = num_reactions;
        reaction_has_external_inputs(num_reactions) = true;
    end
end
reaction_reactants(num_reactions,2) = {reactant_indices};
reaction_reactants_num(num_reactions) = size(reactant_indices,2);
for k=1:reaction_reactants_num(num_reactions)
    reaction_reactants_index(num_reactions,k) = reactant_indices(k);
end

% populate product names and their indices
reaction_products(num_reactions,1)  = {products};
product_indices = [];
for i=1:size(products,2)
    if (i==1)
        product_indices = [lookup_species(products(i))];
    else
        product_indices = [product_indices lookup_species(products(i))];
    end
end
reaction_products(num_reactions,2) = {product_indices};

% store reaction rate
reaction_rates(num_reactions,1)     = {rate};

end