Still cleaning this up, but here are the basic ideas...

The end points of the main beam are given by
\begin{align} {x_1} & = L_{1} \sin\left({\theta}\right) + {x} \\ {y_1} & = L_{1} \cos\left({\theta}\right) \\ {x_2} & = -L_{2} \sin\left({\theta}\right) + {x} \\ {y_2} & = -L_{2} \cos\left({\theta}\right) \\ \end{align}
Here we allow for translation of the trebuchet by x(t) from the origin. Now the position of the beam center of gravity is
\begin{align} {x_b} & = \frac{1}{2} \, {\left(L_{1} - L_{2}\right)} \sin\left({\theta}\right) + {x} \\ {y_b} & = \frac{1}{2} \, {\left(L_{1} - L_{2}\right)} \cos\left({\theta}\right) \\ \end{align} The position of the projectile (m)
\begin{align} {x_m} & = L_{m} \sin\left({\psi} + {\theta}\right) - L_{2} \sin\left({\theta}\right) + \, {x} \\ {y_m} & = L_{m} \cos\left({\psi} + {\theta}\right) - L_{2} \cos\left({\theta}\right) \\ \end{align} The position of the counterweight (M)
\begin{align} {x_M} & = L_{M} \sin\left({\phi} + {\theta}\right) + L_{1} \sin\left({\theta}\right) + \, {x} \\ {y_M} & = L_{M} \cos\left({\phi} + {\theta}\right) + L_{1} \cos\left({\theta}\right) \\ \end{align} The Lagrangian is defined \begin{align} \mathcal{L} = T - V \end{align} where T is the kinetic energy and U is the potential energy of e the system. Will need the velocities
\begin{align} {\dot x_M} & = L_{M} {\left({\dot \phi} + {\dot \theta}\right)} \cos\left({\phi} + {\theta}\right) + L_{1} {\dot \theta} \cos\left({\theta}\right) + \, {\dot x} \\ {\dot x_M} & = -L_{M} {\left({\dot \phi} + {\dot \theta}\right)} \sin\left({\phi} + {\theta}\right) - L_{1} {\dot \theta} \sin\left({\theta}\right) \\ {\dot x_m} & = L_{m} {\left({\dot \psi} + {\dot \theta}\right)} \cos\left({\psi} + {\theta}\right) - L_{2} {\dot \theta} \cos\left({\theta}\right) + \, {\dot x} \\ {\dot x_m} & = -L_{m} {\left({\dot \psi} + {\dot \theta}\right)} \sin\left({\psi} + {\theta}\right) + L_{2} {\dot \theta} \sin\left({\theta}\right) \\ {\dot x} & = {\dot x} \\ \end{align}

\begin{align} {I_b} = m_{b} \frac{1}{12} \, {\left({\left(L_{1} + L_{2}\right)}^{2} + 12 \, L_{b}^{2}\right)} \end{align} \begin{align} {T} = \frac{1}{2} {m} {\left({\dot x_m}^2 + {\dot y_m}^2\right)} + \frac{1}{2} {M} {\left({\dot x_M}^2 + {\dot y_M}^2\right)} + \frac{1}{2} {\left({m_b + m_c}\right)} {\dot x^2} + \frac{1}{2} {I_b} {\dot \theta^2} \end{align}

Define KE and TE here... After considerable algebra (SageMath), the kinetic energy

\begin{gather} L_{1} L_{M} M {\dot \phi} {\dot \theta} \cos\left({\phi}\right) - L_{2} L_{m} m {\dot \psi} {\dot \theta} \cos\left({\psi}\right) + \frac{1}{2} \, L_{M}^{2} M {\dot \phi}^{2} + \frac{1}{2} \, L_{m}^{2} m {\dot \psi}^{2} + L_{M}^{2} M {\dot \phi} {\dot \theta} + L_{m}^{2} m {\dot \psi} {\dot \theta} + L_{M} M {\dot \phi} {\dot x} \cos\left({\phi} + {\theta}\right) + L_{M} M {\dot \theta} {\dot x} \cos\left({\phi} + {\theta}\right) + L_{m} m {\dot \psi} {\dot x} \cos\left({\psi} + {\theta}\right) + L_{m} m {\dot \theta} {\dot x} \cos\left({\psi} + {\theta}\right) + L_{1} M {\dot \theta} {\dot x} \cos\left({\theta}\right) - L_{2} m {\dot \theta} {\dot x} \cos\left({\theta}\right) + \frac{1}{2} \, {\left(2 \, L_{1} L_{M} M \cos\left({\phi}\right) - 2 \, L_{2} L_{m} m \cos\left({\psi}\right) + L_{1}^{2} M + L_{M}^{2} M + L_{2}^{2} m + L_{m}^{2} m + I_{b}\right)} {\dot \theta}^{2} + \frac{1}{2} \, {\left(M + M_{C} + m\right)} {\dot x}^{2} \end{gather} Likewise, the potential energy for each mass \begin{gather} {U} = {\left(L_{M} \cos\left({\phi} + {\theta}\right) + L_{1} \cos\left({\theta}\right) + H\right)} M g + {\left(L_{m} \cos\left({\psi} \\ + {\theta}\right) - L_{2} \cos\left({\theta}\right) + H\right)} g m + \frac{1}{2} \, {\left({\left(L_{1} - L_{2}\right)} \cos\left({\theta}\right) + 2 \, H\right)} g m_{b} \end{gather} The equations of motion are found with: \begin{align*} \frac{d}{dt} \left( \frac{\partial \mathcal{L}}{\partial \dot{q}_i} \right) - \frac{\partial \mathcal{L}}{\partial q_i} &= 0 \end{align*} After substituting for each state variable:

θ \begin{gather} -L_{1} L_{M} M {\dot \phi}^{2} \sin\left({\phi}\right) - 2 \, L_{1} L_{M} M {\dot \phi} {\dot \theta} \sin\left({\phi}\right) + L_{2} L_{m} m {\dot \psi}^{2} \sin\left({\psi}\right) + 2 \, L_{2} L_{m} m {\dot \psi} {\dot \theta} \sin\left({\psi}\right) - L_{M} M g \sin\left({\phi} + {\theta}\right) - L_{m} g m \sin\left({\psi} + {\theta}\right) - L_{1} M g \sin\left({\theta}\right) + L_{2} g m \sin\left({\theta}\right) - \frac{1}{2} \, L_{1} g m_{b} \sin\left({\theta}\right) + \frac{1}{2} \, L_{2} g m_{b} \sin\left({\theta}\right) + {\left(L_{1} L_{M} M \cos\left({\phi}\right) + L_{M}^{2} M\right)} {\ddot \phi} - {\left(L_{2} L_{m} m \cos\left({\psi}\right) - L_{m}^{2} m\right)} {\ddot \psi} + \frac{1}{12} \, {\left(24 \, L_{1} L_{M} M \cos\left({\phi}\right) - 24 \, L_{2} L_{m} m \cos\left({\psi}\right) + 12 \, L_{1}^{2} M + 12 \, L_{M}^{2} M + 12 \, L_{2}^{2} m + 12 \, L_{m}^{2} m + L_{1}^{2} m_{b} + 2 \, L_{1} L_{2} m_{b} + L_{2}^{2} m_{b} + 12 \, L_{b}^{2} m_{b}\right)} {\ddot \theta} + {\left(L_{M} M \cos\left({\phi} + {\theta}\right) + L_{m} m \cos\left({\psi} + {\theta}\right) + L_{1} M \cos\left({\theta}\right) - L_{2} m \cos\left({\theta}\right)\right)} {\ddot x} =0 \end{gather}

φ \begin{gather} L_{1} L_{M} M {\dot \theta}^{2} \sin\left({\phi}\right) + L_{M}^{2} M {\ddot \phi} + L_{M} M {\ddot x} \cos\left({\phi} + {\theta}\right) - L_{M} M g \sin\left({\phi} + {\theta}\right) + {\left(L_{1} L_{M} M \cos\left({\phi}\right) + L_{M}^{2} M\right)} {\ddot \theta} \end{gather}

ψ \begin{gather} -L_{2} L_{m} m {\dot \theta}^{2} \sin\left({\psi}\right) + L_{m}^{2} m {\ddot \psi} + L_{m} m {\ddot x} \cos\left({\psi} + {\theta}\right) - L_{m} g m \sin\left({\psi} + {\theta}\right) - {\left(L_{2} L_{m} m \cos\left({\psi}\right) - L_{m}^{2} m\right)} {\ddot \theta} =0 \end{gather}

x \begin{gather} -L_{M} M {\dot \phi}^{2} \sin\left({\phi} + {\theta}\right) - 2 \, L_{M} M {\dot \phi} {\dot \theta} \sin\left({\phi} + {\theta}\right) - L_{M} M {\dot \theta}^{2} \sin\left({\phi} + {\theta}\right) - L_{m} m {\dot \psi}^{2} \sin\left({\psi} + {\theta}\right) - 2 \, L_{m} m {\dot \psi} {\dot \theta} \sin\left({\psi} + {\theta}\right) - L_{m} m {\dot \theta}^{2} \sin\left({\psi} + {\theta}\right) - L_{1} M {\dot \theta}^{2} \sin\left({\theta}\right) + L_{2} m {\dot \theta}^{2} \sin\left({\theta}\right) + L_{M} M {\ddot \phi} \cos\left({\phi} + {\theta}\right) + L_{m} m {\ddot \psi} \cos\left({\psi} + {\theta}\right) + {\left(L_{M} M \cos\left({\phi} + {\theta}\right) + L_{m} m \cos\left({\psi} + {\theta}\right) + L_{1} M \cos\left({\theta}\right) - L_{2} m \cos\left({\theta}\right)\right)} {\ddot \theta} + {\left(M + M_{C} + m\right)} {\ddot x} \end{gather}
Have to deal with constraint that projectile might initially be sliding on the base of the trebuchet (no friction). Use Lagrange multipliers. \begin{align*} \frac{d}{dt} \left( \frac{\partial \mathcal{L}}{\partial \dot{q}_i} \right) - \frac{\partial \mathcal{L}}{\partial q_i} - {\lambda}\frac{\partial {f}}{\partial q_i} = 0 \end{align*}
The constraint is: \begin{gather} {f} = L_{m} \cos\left(\theta\left(t\right) + \psi\left(t\right)\right) - L_{2} \cos\left(\theta\left(t\right)\right) + H = 0 \\ \end{gather} Since this is equal to zero, we can differentiate twice with respect to time - this puts the constraint in the same form as the other equations. \begin{gather} \frac{\partial^2 {f}}{\partial {t}^2} = {\ddot {f}} = -L_{m} {\left({\dot \psi} + {\dot \theta}\right)}^{2} \cos\left({\psi} + {\theta}\right) + L_{2} {\dot \theta}^{2} \cos\left({\theta}\right) - L_{m} {\ddot \psi} \sin\left({\psi} + {\theta}\right) - {\left(L_{m} \sin\left({\psi} + {\theta}\right) - L_{2} \sin\left({\theta}\right)\right)} {\ddot \theta} \\ \end{gather} \begin{gather} \frac{\partial {f}}{\partial \theta} = -L_{m} {\left({\dot \psi} + {\dot \theta}\right)}^{2} \sin\left({\psi} + {\theta}\right) + L_{2} {\dot \theta}^{2} \sin\left({\theta}\right) + L_{m} {\ddot \psi} \cos\left({\psi} + {\theta}\right) + {\left(L_{m} \cos\left({\psi} + {\theta}\right) - L_{2} \cos\left({\theta}\right)\right)} {\ddot \theta} \\ \frac{\partial {f}}{\partial \psi} = -L_{m} {\left({\dot \psi} + {\dot \theta}\right)}^{2} \sin\left({\psi} + {\theta}\right) + L_{m} {\ddot \psi} \cos\left({\psi} + {\theta}\right) + L_{m} {\ddot \theta} \cos\left({\psi} + {\theta}\right) \\ \end{gather} For the simulation we need to solve simultaneous differential equations described above. There are several distinct phases of the operation:

Phase 1: The projectile is sliding along the deck

\begin{gather} \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} & a_{15}\\ a_{21} & a_{22} & a_{23} & a_{24} & a_{25}\\ a_{31} & a_{32} & a_{33} & a_{34} & a_{35}\\ a_{41} & a_{42} & a_{43} & a_{44} & a_{45}\\ a_{51} & a_{52} & a_{53} & a_{54} & a_{55} \end{bmatrix} \begin{bmatrix} \ddot \theta\\ \ddot \phi\\ \ddot x\\ \ddot \psi\\ \lambda \end{bmatrix} = \begin{bmatrix} b_{1}\\ b_{2}\\ b_{3}\\ b_{4}\\ b_{5} \end{bmatrix} \end{gather}

Phase 2: The projectile is airborne (triple pendulum)

\begin{gather} \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14}\\ a_{21} & a_{22} & a_{23} & a_{24}\\ a_{31} & a_{32} & a_{33} & a_{34}\\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix} \begin{bmatrix} \ddot \theta\\ \ddot \phi\\ \ddot x\\ \ddot \psi\\ \end{bmatrix} = \begin{bmatrix} b_{1}\\ b_{2}\\ b_{3}\\ b_{4}\\ \end{bmatrix} \end{gather}

Phase 3: The projectile is released (double pendulum)

\begin{gather} \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} \begin{bmatrix} \ddot \theta\\ \ddot \phi\\ \ddot x\\ \end{bmatrix} = \begin{bmatrix} b_{1}\\ b_{2}\\ b_{3}\\ \end{bmatrix} \end{gather} After release + \begin{align*} {\ddot x_m} = 0\\ {\ddot y_m} = -g \end{align*}

...all programmed in Javascript.