function f = pendulum(t,y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% f = pendulum(t,y)                                      %
%  -- Evaluation of 1st Order System of ODEs             %
%     modeling a mechanical pendulum--                   %
%                                                        %
% pendulum.m computes                                    %
% f(t,y) =                                               %
%    [y(2);u - c*y(2) - m*g*y(1))/(m*l)]                 %
%        (mydata.linarflag == 0)                         %
%  or                                                    %
%    [y(2);u - c*y(2) - m*g*sin(y(1)/(m*l)]              %
%        (mydata.linarflag == 1)                         %
% where the ODE y' = f(t,y) describes the motion         %
% of a damped, driven pendulum.                          %
%                                                        %
% input   t        Real nonnegative number denoting      %
%                  the time at which f(t,y) is computed. %
%                                                        %
%         y        Real column vector                    %
%                  y = [theta(t);theta'(t)]              %
%                                                        %
% output  f        Real column vector denoting           %
%                  f = [theta(t);theta'(t)]              %
%                                                        %
% Yalin Sagduyu and Dianne P. O'Leary                    %            
% 12/02 and 01/03                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global mydata 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% mydata.m              is the mass                      %
% mydata.l              is the length of the pendulum    %
% mydata.g = 9.81       is the acceleration of gravity   %
% mydata.c              is the damping constant          %
% mydata.forcedflag = 0 indicates no applied force       %
%                   = 1 indicates applied force          %
%                   = 2 indicates applied force          %
%                        with control                    %
% mydata.linearflag = 0 indicates the nonlinear model    %
% mydata.forceangle =   forceangle.                      %
%                       u(t) = m*g*sin(forceangle)       %
%                        when forcedflag == 1            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (mydata.forcedflag == 0)
    u = 0;
elseif (mydata.forcedflag == 1)
    u = mydata.m * mydata.g * sin(mydata.forceangle);
else
    u = mydata.m * mydata.g * sin(mydata.forceangle)...
        + mydata.m * mydata.l * (mydata.Bc * y(2));
end
    
if (mydata.linearflag == 0)
    tt = sin(y(1));
else
    tt = y(1);
end

f(1,1) = y(2);
f(2,1) = (u - mydata.c*y(2) - mydata.m * mydata.g * tt) ...
                   /(mydata.m * mydata.l);