function dotq = central_force (t,q)
% 中心力が作用する質点の運動方程式
    global a;

    r = q(1);    theta = q(2);
    v = q(3);    omega = q(4);

    dotr = v;
    dottheta = omega;
    dotv = r*omega^2 - a/(r^2);
    dotomega = - 2*v*omega/r;

    dotq = [dotr; dottheta; dotv; dotomega];
end
