reaction(names,stoichiometry,rates,equilibria)

function react=reaction(names,stoichiometry,rates,equilibria)
%
%    react=reaction(names,stoichiometry,[rates],[equilibria])
% 
% reaction class constructor
% % Creates a reaction object containing fields describing a chemical reaction
% fields
%
% reagents are numbered sequentially 1-n but an be referred to by their:
%   reagentnames{} -- cell array
% 
%   stoichiometry -- numeric array 
%
% each row indicates the
% stoichiometry for transforming one set of components into another
% ie A + B --> C would be encoded as
% 
%       A   B   C
%       -1  -1  1
% 
% Note that reactions are *directional* - this means a reaction system could be encoded in 
% more than one way depending upon which were though of as 'reagents' and 'products'
% Also some encodings may be acceptable thermodynamically, but give incorrect kinetics
% 
%   rates -- column vector (optional)
%
% specifies rate constants for each row of the stoichiometry matrix.
% Rate constants default to -1 to indicate the rate is unknown
%
%   equilibria -- column vector (optiona)
%
% specifies the position of the equilibrium for each row of the stoichiometry 
% matrix - note that this also specifes the reverse rate constant so that rows for
% reverse reactions are not required. Equilibrium constants default to Inf. to indicate that
% the reverse reaction is not significant.
% 
%

if (nargin==1) & (isa(names,'reaction'))
    react=names;
else
    % defaults
    if nargin<4 equilibria=[]; end
    if nargin<3 rates=[]; end
    if nargin<2 stoichiometry=[]; end
    if nargin<1 names=cell(0); end

    % organize rows and columns
    names=names(:).';
    rates=rates(:);
    equilibria=equilibria(:);
    
    % calculate sizes
    species=max([size(names,2), size(stoichiometry,2)]);
    steps=max([size(stoichiometry,1), size(rates,1), size(equilibria,1)]);

    % make template and fill it
    react.reagent=cell(1,species);
    react.reagent(1:size(names,2))=names;

    react.stoichiometry=zeros(steps,species);
    if size(stoichiometry,1)~=0
        react.stoichiometry(1:size(stoichiometry,1),1:size(stoichiometry,2))=stoichiometry;
    end
    
    % default reaction rates are -1 - indicating they are not known
    react.rate=-1*ones(steps,1);
    if size(rates,1)~=0
        react.rate(1:size(rates,1),1)=rates;
    end
    % default equilibrium constants are infinite - there is no reverse reaction
    react.equilibrium=inf.*ones(steps,1);
    if size(equilibria,1)~=0
        react.equilibrium(1:size(equilibria,1))=equilibria;
    end

    react=class(react,'reaction');
end