function vprime = damppend(t,vin)

global r;
global a;
global b;
global k;


x = vin(1);
y = vin(2);
z = vin(3);

% compute derivatives for damped pendulum, 
% from K Tung p. 227, b=omega^2, k= gamma, r=omega_D, a= D 

xprime = y;
yprime = -b*sin(x) - k*y + a*cos(z);
zprime = r;

% map back to vector format
vprime=[xprime;yprime;zprime];

end