M'=M+\sum{\Delta m_i}

M'V'=MV+\sum{\Delta m_i v_i}

V' = V \left ( \frac{1}{1+\epsilon} \right ) +

M'R' = M'R + \sum{\Delta m_i \left ( r_i - R \right )}=MR + \sum{\Delta m_i r_i }

S'=S+\sum{\Delta m_i \left ( r_i -R \right ) \times v_i}

However, a better treatment is to conserve angular momentum.

L=MR \times V

J=L+S

J'=J +\sum{\Delta m_i r_i \times v_i}

S'=S + MR \times V - M'R' \times V' + \sum{\Delta m_i r_i \times v_i}

S'=S + MR \times V - M'R' \times V' + \sum{\Delta m_i r_i \times v_i}

S' = S+MR \times V - \left( MR + \sum{\Delta m_i r_i} \right ) \times \frac{MV + \sum{\Delta m_i v_i}}{M+\sum{m_i}}+\sum{\Delta m_i r_i \times v_i}

which if we approximate to first order in \frac{\sum{\Delta m_i}}{M}

S' \approx S+MR \times V - \left( MR + \sum{\Delta m_i r_i} \right ) \times \left ( V + \frac{\sum{\Delta m_i v_i}}{M} \right ) \left (1 -\frac{\sum{\Delta m_i}}{M} \right )+\sum{\Delta m_i r_i \times v_i}

S' \approx S+MR \times V - MR \times V - \sum{\Delta m_i r_i } \times V - R \times \sum{\Delta m_i v_i} + R \times \sum{\Delta m_i} V+\sum{\Delta m_i r_i \times v_i}

S' \approx S +\sum{\Delta m_i \left (r_i - R \right ) \times \left (v_i - V \right )}

which implies our method is only accurate if V << v_i and \sum{\Delta m_i} << M