% Funcao destinada ao uso da funcao ODE45 para a determinaçao numerica
% da equaçao diferencial associada a dinamica do pendulo simples
%
% Equacao diferencial:  D2theta+c/(M*L^2)*Dtheta+g/L*sin(theta)=0
%       theta: Orientacao do pendulo a partir da posiçao de repouso
%       Dtheta: Primeira derivada de theta em ordem ao tempo
%       D2theta: Segunda derivada de theta em ordem ao tempo
%       M: Massa da esfera na extremidade do pendulo
%       L: Comprimento da haste do pendulo
%       c: Coeficiente de fricçao da junta
%       g: Aceleracao da gravidade
%
% Como a funçao ODE45 so suporta derivadas de primeira ordem, a tecnica de
% reduçao de ordem da equaçao teve que ser usada, convertendo-a em duas
% equacoes.
% p1=theta
% p2=Dtheta
% Dp2=D2theta=-c/(M*L)*Dtheta-g/L*sin(theta) = -c/(M*L)*p2-g/L*sin(p1)
% Dp=D(p1;p2)= [ p2; -c/(M*L^2)*p2-g/L*sin(p1) ] = [f1(t,p);f2(t,p)]=f(t,p)
%
% Condicoes inicias: p(0)=[p1(0);p2(0)]=[theta(0);Dtheta(0)]
%
% (c) MRS, 2005/09/14
% Referencia: Bucknell University, Using ODE45 (MatLab Help)

function dpdt=theta_odefun(t,p)

pendulum_data

dpdt=zeros(size(p));
dpdt(1)=p(2);
dpdt(2)=-1/(M*L^2)*(c*p(2)+k*p(1)+M*g*L*sin(p(1)));
