function dots = rotation_quaternion_ODE_params (t,s, inertiaJ, inertiaJinv, alpha)
    q = s(1:4); p = s(5:8);
    dotq = p;
    H = matrix_H(q);
    r = quaternion_constraint(q,p,alpha);
    Hp = H*p; JHp = inertiaJ*Hp;
    dotp = -r*q - 2*H'*inertiaJinv*(cross(Hp,JHp)-(1/4)*external_torque(t));
    dots = [dotq;dotp];
end
