function [sys,x0,str,ts] = pendeom(t,x,u,flag,g,l,HalfRhoSCd,mass,theta0,thetadot0)
%   S-file for pendulum (with corrected damping term, 1/12/01)
%   This is an S-file subsystem that models a simple pendulum.
%   g is the acceleration of gravity. 
%   l is the length of the pendulum.
%   mass is the mass of the ball that forms the massive part of the pendulum.
%   theta0, thetadot0 are the initial condition on theta and thetadot.
switch flag,      
  case 0,         % Initialization
    [sys,x0,str,ts]=mdlInitializeSizes(theta0,thetadot0);
  case 1,         % Compute derivatives of continuous states
    sys=mdlDerivatives(t,x,u,g,l,HalfRhoSCd,mass) ;
  case 2,
    sys=mdlUpdate(t,x,u);
  case 3,
    sys=mdlOutputs(t,x,u);  % Compute output vector
  case 4,                   % Compute time of next sample
    sys=mdlGetTimeOfNextVarHit(t,x,u);
  case 9,                   % Finished. Do any needed 
    sys=mdlTerminate(t,x,u);
  otherwise                 % Invalid input
    error(['Unhandled flag = ',num2str(flag)]);
end

%***********************************************************
%*                    mdlInitializeSizes                   *
%***********************************************************
function [sys,x0,str,ts]=mdlInitializeSizes(theta0,thetadot0)
% Return the sizes of the system vectors, initial 
% conditions, and the sample times and offets.
sizes = simsizes;   % Create the sizes structure
sizes.NumContStates  = 2;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 2;
sizes.NumInputs      = 0;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;   % at least one sample time is needed
sys = simsizes(sizes);      % load sys with the sizes structure
x0  = [theta0,thetadot0];    % Specify initial conditions for all states 
str = [];    % str is always an empty matrix
ts  = [0 0]; %initialize the array of sample times


%***********************************************************
%*                    mdlDerivatives                       *
%***********************************************************
function sys=mdlDerivatives(t,x,u,g,l,HalfRhoSCd,mass)
% Compute derivatives of continuous states
x_dot = [x(2); -(g/l)*sin(x(1))-HalfRhoSCd*l*x(2)*x(2)*sign(x(2))/mass] ;
sys = [x_dot(1),x_dot(2)] ;

%***********************************************************
%*                    mdlUpdate                            *
%***********************************************************
function sys=mdlUpdate(t,x,u)
% Compute update for discrete states. If necessary, check for
% sample time hits.
sys = [];    % Empty since this model has no discrete states.

%***********************************************************
%*                    mdlOutputs                           *
%***********************************************************
function sys=mdlOutputs(t,x,u)
% Compute output vector given current state, time, and input
sys = x ;

%***********************************************************
%*                    mdlGetTimeOfNextVarHit               *
%***********************************************************
function sys=mdlGetTimeOfNextVarHit(t,x,u)
% Return the time of the next hit for this block.  Note that 
% the result is absolute time.  Note that this function is 
% only used when you specify a variable discrete-time sample
% time [-2 0] in the sample time array in sampleTime = 1;    
sys = [] ;

%************************************************************
%*                    mdlTerminate                          *
%************************************************************
function sys=mdlTerminate(t,x,u)
% Perform any necessary tasks at the end of the simulation
sys = [];

