function dotq = damped_pendulum (t,q)
% equation of motion of damped pendulum
    global mass; global length; global grav;
    global viscous;
    theta = q(1);
    omega = q(2);
    dottheta = omega;
    dotomega = 1/(mass*length*length)* ...
      (external_torque(t) - mass*grav*length*sin(theta) - viscous*omega);
    dotq = [dottheta; dotomega];
end
