function m = initModel(nX, nU, nY, nP)
%
% Initializes the data structure
%
% 
% This work is licensed under the Creative Commons Attribution-Noncommercial-No Derivative 
% Works 3.0 United States License. To view a copy of this license, visit 
% http://creativecommons.org/licenses/by-nc-nd/3.0/us/ or send a letter to Creative Commons,
% 171 Second Street, Suite 300, San Francisco, California, 94105, USA.

    % Intialize the sizes of different matrices
    m.nX   = nX;
    m.nU   = nU;
    m.nY   = nY;
    m.nP   = nP;
    m.nR   = 0;
    
    % Initialize the stoichiometry matrix which is size number of states 
    % by number of reactions (zero to start).
    m.spinds.i_S = zeros(3*nX, 1);
    m.spinds.j_S = zeros(3*nX, 1);
    m.spinds.x_S = zeros(3*nX, 1);
    m.spinds.n_S = 0;
    
    % Initializes a new kinetic model of the type:
    %   dx/dt = a1 x + a2 kron(x,x) + b1 u + b2 kron(x,x) + k
    %   y     = c x
    % The arguments are the lengths of the X, U, Y, and P  vectors repectively
    
    m.a1  = sparse(nX,       nX);
    m.a2  = sparse(nX,    nX*nX);
    m.b1  = sparse(nX,       nU);
    m.b2  = sparse(nX,    nU*nX);
    m.k   = sparse(nX,        1);
    m.c   = zeros (nY,       nX);
    m.ic  = zeros (nX,        1);
    m.p   = zeros( nP,        1);
    m.a1P = sparse(nP*nX,    nX);
    m.a2P = sparse(nP*nX, nX*nX);
    m.b1P = sparse(nP*nX,    nU);
    m.b2P = sparse(nP*nX, nU*nX);
    m.kP  = sparse(nP*nX,     1);
    
    m.spinds.i_a1 = zeros(floor(nX^1.2), 1);
    m.spinds.j_a1 = zeros(floor(nX^1.2), 1);
    m.spinds.x_a1 = zeros(floor(nX^1.2), 1);
    m.spinds.n_a1 = 0;
    
    m.spinds.i_a1P = zeros(floor(nX^1.2), 1);
    m.spinds.j_a1P = zeros(floor(nX^1.2), 1);
    m.spinds.x_a1P = zeros(floor(nX^1.2), 1);
    m.spinds.n_a1P = 0;

    m.spinds.i_a2 = zeros(floor(nX^1.2), 1);
    m.spinds.j_a2 = zeros(floor(nX^1.2), 1);
    m.spinds.x_a2 = zeros(floor(nX^1.2), 1);
    m.spinds.n_a2 = 0;
    
    m.spinds.i_a2P = zeros(floor(nX^1.2), 1);
    m.spinds.j_a2P = zeros(floor(nX^1.2), 1);
    m.spinds.x_a2P = zeros(floor(nX^1.2), 1);
    m.spinds.n_a2P = 0;

    m.spinds.i_b1 = zeros(floor(nX^1.2), 1);
    m.spinds.j_b1 = zeros(floor(nX^1.2), 1);
    m.spinds.x_b1 = zeros(floor(nX^1.2), 1);
    m.spinds.n_b1 = 0;
    
    m.spinds.i_b1P = zeros(floor(nX^1.2), 1);
    m.spinds.j_b1P = zeros(floor(nX^1.2), 1);
    m.spinds.x_b1P = zeros(floor(nX^1.2), 1);
    m.spinds.n_b1P = 0;

    m.spinds.i_b2 = zeros(floor(nX^1.2), 1);
    m.spinds.j_b2 = zeros(floor(nX^1.2), 1);
    m.spinds.x_b2 = zeros(floor(nX^1.2), 1);
    m.spinds.n_b2 = 0;
    
    m.spinds.i_b2P = zeros(floor(nX^1.2), 1);
    m.spinds.j_b2P = zeros(floor(nX^1.2), 1);
    m.spinds.x_b2P = zeros(floor(nX^1.2), 1);
    m.spinds.n_b2P = 0;
    
    m.xNames = {};
    m.yNames = {};
    m.uNames = {};
    m.pNames = {};
end