Hi Valérian,
The design of primary/secondary equations in ODE is based on the idea that the primary equations computes derivatives of primary state which are based only on the state itself. Secondary (aka additional) equations can depend on both primary and secondary states, they usually compute the derivatives of secondary state.
As your equations do depend on secondary state, they belong to secondary equations, BUT they also change the primary state derivatives. In a perfect world, I would say this mean primary state is not rich enough and should be extended so your equations really belong to the mrimary state (after all, they do change primary state derivatives). However, this would clearly not suit your needs because the differential equations of the motion implemented in Orekit only depend on orbit and attitude.
However, this case despite being very rare can be handled. If you look at the javadoc for SecondaryODE.computeDerivatives, you will see that the primaryDot array can be changed by the secondary equation. The javadoc for this parameter reads:
primaryDot array containing the derivative of the primary state vector
(the method is allowed to change the derivatives here, when the additional
equations do have an effect on the primary equations)
This feature is available at Orekit level too (see lines 729-733 in AbstractIntegratedPropagator). If AdditionalEquations.computeDerivatives return an array instead of returning null, then this array will be considered to be an additional contribution to the primary state derivatives (i.e. the orbit).
So instead of implementing your equations as a ForceModel, you can implement it fully as an AdditionalEquations, andhave this additional equation have a side effect on orbit because it is an acceleration. You will have to deal with the orbit type by yourself, though. If you are propagating in Cartesian parameters, it is simple: as the state is (pos, vel), the derivative of the state is (vel, acc) so your additional equation should return an array where the 3 first elements are0 and the 3 last elements are the contribution of this model to the acceleration. This vector will be added to the other known contributions (typically Keplerian attraction, gravity field, …). If you want your equation to be usable also when propagating with other orbit types, it will be more complex, you will have to use compute the Jacobian of the orbit with respect to Cartesian parameters and multiply by the (vel, acc) vector so you get the derivatives or the orbital parameters (basically this is how Gaussian equations work). You can look at the NumericalPropagator.Main internal class and its addNonKeplerianAcceleration method for a way to convert an accelration to a primary state derivative.