function [statedots, cm, grf, mforces, mforcesCE, mtorqueDOF, stick, outputsFEA, Premind] = legs2d(times, states, controls, biofile);

% Ordinary Differential Equation for 2D musculoskeletal model : dx/dt = f(x,u,t);
%
% This model is identical to the one described in Gerritsen et al., Motor Control, 1998, but muscle models are
% implemented as described in McLean et al., J Biomech Eng 2003. 
%
% Usage: [statedots, cm, grf, mforces, stick, mforcesCE, Premind, outputsFEA] = legs2d(times, states, controls, footm);
%
% Source code for the MEX function legs2d.dll is in legs2d.c
%
% Input:
%	time		Times at which the ODE must be evaluated (1 x Nt)
%	states      States (50 x Nt) : q(9), qdot(9), Lce(16), act(16)
%	controls	Stims (16 x Nt): stimulation level for 16 muscles
%   biofile    (optional) Name of .bio file and _i file.  File must have 16 muscles.
%   The default one is legs2d.bio and legs3d_i in the same directory if
%   there is no this input argument. In the foot contact model, the right
%   foot comes before the left foot
%
% Output
%   statedots       State derivatives (NSTATES x Nt)
%   cm          	(optional) Center of mass motion (4 x Nt); the four columns are: x,y,vx,vy
%   grf         	(optional) Ground reaction force ny spring damper elements (6 x Nt); the six columns are:
%                   Fxleft, Fyleft, COPleft, Fxright, Fyright, COPPright. Note: returns zero if footm = 1.
%   mforces         (optional) Muscle forces (16 x Nt); forces in the 16 muscles
%   mforcesCE       (optional) Muscle CE forces (16 x Nt); contractile element forces in the 16 muscles 
%   mtorqueDOF      (optional) the total muscle torque on each DOF (9 x Nt)
%   stick           global XYZ coordinates of all joints (7) and all points (10*2) of each GRF elements: 3*(7+2*10)
%                   (81 x Nt)
%   outputsFEA      (optional) forces and moments at ankle (12 x Nt); returns ground contact forces at rigth and left ankles due to
%                   contact model specified in right.m and left.m, respectively:  Fxright, Fyright, Fzright, Mxright, 
%                   Myright, Mzright, Fxleft, Fyleft, Fzleft, Mxleft, Myleft, Mzleft. Note: returns zero if footm = 0
%   Premind         (optional) returns error messages

% Notes:
%   (1) When used with an ODE solver, legs2d is usually evaluated at a single time (Nt=1).
%       The capability to evaluate at multiple times is needed for solving optimal control
%       problems or boundary value problems using collocation.
%       In the demo program (simulate.m), we use Nt=100 at the end to generate results with
%       constant time interval.  This could also have been done by calling the ODE solver once for
%       each time interval, this would have been somewhat slower though because the ODE solver has
%       to be restarted multiple times.
%   (2) States are defined as follows:
%       x(1)		X coordinate of trunk center of mass (m), increases with forward motion
%       x(2)		Y coordinate of trunk center of mass (m), increases with upward motion
%       x(3)		orientation of trunk (rad), zero when model is lying horizontally facing up, -pi/2 when standing
%       x(4)		left hip angle (rad), zero when standing, positive for counterclockwise rotation of femur (flexion)
%       x(5)		left knee angle (rad), zero when standing, positive for counterclockwise rotation of tibia (extension)
%       x(6)        left ankle angle (rad), +pi/2 when standing, increasing with dorsiflexion
%       x(7-9)      right side joint angles, defines as x(4-6)
%       x(10-18)    generalized speeds, time derivatives of x(1-9)
%       x(19-34)    contractile element (CE) lengths of the muscles, normalized to each muscle's LCEopt.
%                   muscles are ordered as in .bio file
%       x(35-50)    active state of the muscles.  Full activation is 1.0.
%   (3) COP (center of pressure) is zero when vertical force (Fy) is zero, and unreliable when Fy is very small.
%       It is recommended to ignore grf results when Fy is less than 1 N or so.
%   (4) all the output arguments except statedots are optional. For example, in order to
%       obtain stick use
%       [statedots, cm, grf, mforces, stick] = legs2d(times, states, controls, footm)
%   (5) user-defined contact models can be used by setting footm = 1. In this case, legs2d calls the user-defined matlab
%       function right.m and left.m for contact models at the right and left feet, respectively. legs2d provides the position and velocity of
%       the ankle joint, and the orientation and angular velocity of the foot to right.m and left.m, which should return the force and moment at 
%       the respective ankle joint. right.m and left.m should look like 
%       function F = right(q)
%       where:
%       q   ankle and foot kinematics (12 x 1): global ankle position (X, Y and Z components), global foot orientation (consecutive rotations about X, Y
%           and Z), global ankle velocity (X, Y and Z components), global foot angular velocity (X, Y and Z components) 
%       F   force and moments at the ankle (6 x 1): global force at ankle (X, Y and Z components), global moment at ankle (X, Y and Z components)    