The equations governing the velocity $\vec v_G$ of the center of mass and the angular momentum $\vec L_G$ about the center of mass are \[\left\{\begin{array}{rl}\displaystyle M{d\vec v_G\over dt} & = -Mg \vec e_z + \vec N , \\ \displaystyle {d\vec L_G\over dt} & = \big(-h \vec e_Z\big) \times\vec N . \end{array} \right. \] Here, $M$ is the mass of the top, $\vec N$ is the force from the floor at the axis tip, $g$ is the gravitational acceleration, $\vec e_z$ is the unit vector in the $z$ direction (vertically upward), $h$ is the distance between the axis tip $P$ and the center of mass $G$, and $\vec e_Z$ is the unit vector along the $Z$ axis of the top.
The reaction force $\vec N$ consists of a tangential friction force $\vec F$ and a vertical normal force $\vec R$ \[ \vec N = \vec F + \vec R\, ; \qquad \vec F \,\perp\, \vec e_z, \quad \vec R = R\, \vec e_z\,. \] The friction force $\vec F$ is assumed to be proportional to the velocity $\vec v_P$ of the axis tip $P$ \[ \vec F = -k\, \vec v_P\,. \] Here, $k$ is the friction (resistance) coefficient.
On the other hand, the vertical normal force $R$ is determined by the condition that the tip $P$ is on the floor. This is equivalent to requiring that $\vec v_P$ lies in the $xy$ plane: \[ \vec v_P \cdot \vec e_z = 0\,;\qquad \vec v_P = {d\over dt}\big( \vec r_G -h\vec e_Z\big) = \vec v_G - h\,\vec\omega\times\vec e_Z \] From this, the velocity $\vec v_G$ of the center of mass is expressed using the angular velocity $\vec\omega$, and substituting into the $z$-component of the equation of motion for $\vec v_G$, $R$ can be expressed in terms of $\vec\omega$ and its time derivative.
For the detailed derivation of the equations used in the simulation, see the computation notes(Japanese).