Rigid body dynamics: modelling a polygon bounce off ground I'm currently making a physics simulator, but I'm having some trouble making a polygon bounce off the floor. 
I know that collisions are normally modeled as described below by ja72, but I did it somewhat differently.  
What I'm trying to do: 
When the polygon touches (penetrates) the ground, I want to apply a linear impulse $\overrightarrow{P}$ at the contact point such that the energy $E'$ of the polygon after the bounce is some elasticity constant $\epsilon\in[0,1]$ times the energy $E$ of the polygon before the bounce.
$$E'=\epsilon E$$ 
My calculations: 
Consider a cartesian coordinate system where the floor is $y=0$. Then $\overrightarrow{P}(0,p)$ where we expect $p$ to be positive.
Define:


*

*$m$: the mass of the polygon

*$I$: the rotational inertia

*$\overrightarrow{v}$: the initial velocity

*$\overrightarrow{v'}$: the velocity after the collision

*$\overrightarrow{\omega}$: the initial rotational velocity (a vector parallel to the $z$-axis)

*$\overrightarrow{\omega'}$: the rotational velocity after the collision

*$\overrightarrow{r}$: a vector going from the polygon's center of mass to the contact point


Furthermore, for a vector $\overrightarrow{a}$, we denote its length by $|\overrightarrow{a}|$ and we can index it as $a[0]$ or $a[1]$ to get the $x$- and $y$-component (scalars, not vectors), respectively.  
$m$, $I$, $\overrightarrow{v}$, $\overrightarrow{\omega}$ and $r$ are given, and we want to calculate the magnitude $p$ of the impulse $\overrightarrow{P}$.
Because $\overrightarrow{v'}=\overrightarrow{v}+\frac{\overrightarrow{P}}{m}$ and $\overrightarrow{\omega'} = \overrightarrow{\omega} + \frac{\overrightarrow{r}\times\overrightarrow{P}}{I}$, we can calculate:
$$\frac{1}{2}m\cdot\overrightarrow{v'}^2 + \frac{1}{2}I\cdot\overrightarrow{\omega'}^2 = E' = \epsilon E = \epsilon\cdot\left(\frac{1}{2}m\cdot\overrightarrow{v}^2 + \frac{1}{2}I\cdot\overrightarrow{\omega}^2\right)$$
$$\Leftrightarrow m\cdot\left(\overrightarrow{v}+\frac{\overrightarrow{P}}{m}\right)^2 + I\cdot\left(\overrightarrow{\omega}+\frac{\overrightarrow{r}\times\overrightarrow{P}}{I}\right)^2 = \epsilon\cdot\left(m\cdot\overrightarrow{v}^2 + I\cdot\overrightarrow{\omega}^2\right)$$
$$\Leftrightarrow m\cdot\overrightarrow{v}^2 + 2\overrightarrow{v}\cdot\overrightarrow{P} + \frac{\overrightarrow{P}^2}{m} + I\cdot\overrightarrow{\omega}^2 + 2\overrightarrow{\omega}\cdot(\overrightarrow{r}\times\overrightarrow{P}) + \frac{(\overrightarrow{r}\times\overrightarrow{P})^2}{I} = \epsilon\cdot\left(m\cdot\overrightarrow{v}^2 + I\cdot\overrightarrow{\omega}^2\right)$$
And now, we can use the fact that $\overrightarrow{P}(0,p)$ to calculate the dot and cross products:
$$0 = (1-\epsilon)\cdot\left(m\cdot\overrightarrow{v}^2 + I\cdot\overrightarrow{\omega}^2\right) + 2p\cdot v[1] + \frac{p^2}{m} + 2p\cdot r[0]\cdot|\overrightarrow{\omega}| + \frac{(p\cdot r[0])^2}{I}$$
$$\Leftrightarrow \left[\frac{1}{m}+\frac{r[0]^2}{I}\right]\cdot p^2 + \left[2v[1]+2r[0]\cdot|\overrightarrow{\omega}|\right]\cdot p + (1-\epsilon)\cdot\left[m\cdot\overrightarrow{v}^2 + I\cdot\overrightarrow{\omega}^2\right]$$ 
This is a quadratic equation in $p$.
If we set:


*$A = \frac{1}{m}+\frac{r[0]^2}{I}$

*$B = 2v[1]+2r[0]\cdot|\overrightarrow{\omega}|$

*$C = (1-\epsilon)\cdot\left[m\cdot\overrightarrow{v}^2 + I\cdot\overrightarrow{\omega}^2\right]$

We can calculate the determinant $D = B^2-4AC$.
Now, it is clear that $p=\frac{-B\pm\sqrt{D}}{2A}$.
Actually, $p=\frac{-B+\sqrt{D}}{2A}$, as $B$ must be negative to enforce a collision, and the sign of $v[1]$ must be reversed.   
When I run this in a simulation, it looks very natural.
However, it always glitches when $D<0$.
I've tried to resolve this problem in countlessly many ways $(*)$, but it never worked out well.
Can anyone check whether my method and calculations are correct?
And if not, how can I model this collision instead?
Here is a snippet of my code (in Python):
class POLYGON:
    def bounce(self):
        #Checks whether the polygon touches the ground and makes it rebound if necessary

        low_nodes = []    #makes a list of the nodes with y<0
        for node in self.nodes:
            if node[1] <= 0:
                low_nodes.append(node)

        if len(low_nodes) >= 1:
            speed = self.speed
            self.translate([0, -2*lowest_node[1])    #lifts the lowest node above the ground

            for node in low_nodes:
                r = node - self.pos    #calculates r

                A = 1/self.mass + r[0]**2/self.rot_inertia
                B = 2*self.speed[1] + 2*self.rot_speed*r[0]
                C = (1-restituence_constant) * (self.mass*np.linalg.norm(self.speed)**2 + self.rot_inertia*self.rot_speed**2)
                D = B**2 - 4*A*C

                if D<0:
                    impulse = [0, -B/(2*A)]
                else:
                    impulse = [0, (-B-np.sqrt(D)) /(2*A) /len(low_nodes)]
                self.applyLinearImpulse(impulse, node)

$(*)$ 
For example: 
using the absolute value of $D$ 
regarding $p$ as a complex number and applying the impulse $\overrightarrow{P}(\Im(p), \Re(p))$ 
setting $p=\frac{-B}{2A}$ if $D<0$ 
...
 A: The physics of contact can be found at many online resources, so I am just going to summarize below. 
The impulse $\vec{P}$ needs to decomposed into the magnitude $p$ and the known contact normal direction $\vec{n}$, such that $\vec{P} = p\, \vec{n}$. The impulse magnitude is found from the law of collision that states that the relative velocity at the contact point and along the contact normal after the collision is a fraction of the relative impact velocity
$$ v_{\rm relative}' = -\epsilon \, v_{\rm impact} $$
mathematically the above is expressed as follows when an object impacts the ground (single body impact).
$$ \vec{n} \cdot (\vec{v}'+\vec{\omega}' \times \vec{r}) = -\epsilon \;\; \vec{n} \cdot ( \vec{v} + \vec{\omega} \times \vec{r}) $$
where $\vec{v}$ and $\vec{v}'$ are expressed at the center of mass before and after the impact, and the rest according to the author's post.
State the equations of motion & kinematics to get equation 8-18 from the linked paper that gives the impact magnitude

$$ p = \frac{ -(1+\epsilon)\; \vec{n}\cdot \vec{v} }{
 \frac{1}{m} + \vec{n} \cdot \left( I^{-1} (\vec{r} \times \vec{n}) \right) \times \vec{r} }$$

The final velocity will be adjusted then accordingly
$$ \begin{aligned}
  \vec{P} & = p \,\vec{n} \\
  \vec{v}' & = \vec{v} + \tfrac{1}{m} \vec{P} \\
  \vec{\omega}' & = \vec{\omega} + I^{-1} \left( \vec{r} \times \vec{P} \right) \end{aligned} $$

I have also linked an answer to a similar question in 2D.
