Solution (source code)

= Solution

In a frame with unit <angular velocity>, transforming the <acceleration> introduces the <Coriolis acceleration> and <centrifugal acceleration>. Radiation reduces the attraction of $M_1$ by the factor $1-\beta_1$, so
$$
\ddot{\mathbf r}+2\hat{\mathbf z}\times\dot{\mathbf r}
+\hat{\mathbf z}\times(\hat{\mathbf z}\times\mathbf r)
=-\frac{\mu_1(1-\beta_1)}{r_1^3}(\mathbf r-\mathbf r_1)
-\frac{\mu_2}{r_2^3}(\mathbf r-\mathbf r_2).
$$
With the effective potential
$$
U=\frac12(x^2+y^2)+\frac{\mu_1(1-\beta_1)}{r_1}+\frac{\mu_2}{r_2},
$$
its Cartesian components are exactly
$$
\boxed{\ddot x-2\dot y=U_x,
\qquad
\ddot y+2\dot x=U_y,
\qquad
\ddot z=U_z}.
$$
This is the <photogravitational restricted three-body problem>.