CM: Need to recover the Hamiltonian, knowing conserved quantities and information about the EOM, possibly via action-angle coordinates

Statement of the problem:

I have a system with 2 degrees of freedom and I have found two independent conserved quantities, without knowledge of the Hamiltonian. I'm looking for a method to recover a Hamiltonian that describes the system.

Action-angle variables seems like a promising approach, but most of the discussions that I've found on this assume knowledge of the Hamiltonian and at least two of the coordinates $(p,q) \rightarrow (I, \phi)$ that represent the canonical transformation, or the generating function, etc.

I also recognize that some guess work may be necessary if no methodical procedure exists, but advice on "guided" guesses would also be helpful. Maybe simply treating the conserved quantities as the action variables and attempting to guess a generating function for the canonical transformation could work...

Further info

I have a coupled set of 4 first order nonlinear differential equations which describe the dynamics in 4D phase space. To be more specific, these equations take the form

\begin{align} \dot{X} = f_1 (Y,S)\\ \dot{Y} = f_2 (X,Z)\\ \dot{Z} = f_3 (X,Y)\\ \end{align}

where $X, Y, Z$ are functions of time and I omit the equation for $\dot{S}$ because there is a hidden relationship $S^2 = X^2 + Y^2 + Z^2$ (these are the Stokes variables). I have no knowledge or "feel" for which coordinates might describe canonical positions or momenta at this point. As stated, I have independent conserved quantities

\begin{align} \dot{J_1}(Y,Z) = 0\\ \dot{J_2}(S,Y,Z) = 0 \end{align}

A direct method for finding a Hamiltonian is of course desirable, but any ideas or direction towards relevant reference material would also be greatly appreciated.

-
You may be interested in Qmechanic's answer on constructing the Hamiltonian from such given PDEs, especially the second part. Unfortunately, you will still have to find the integrating factor for their inexact differential by hand, and this will, I think, heavily depend on the type PDE you are dealing with. –  ACuriousMind Jul 6 at 22:54
Some variation of this might be the way that I need to go. However, it isn't really clear how to do this, because as mentioned, I don't have knowledge of the canonical coordinates $(p,q)$ in this coordinate system. The only assured quantities I have are the conserved $J_i$, which is why I mention action-angle variables. That coordinate system at least directly depends on the $J_i$, but I still don't seem to have enough information to correctly transform to action-angle variables, although I could be wrong. –  Jon Rayner Jul 7 at 10:50
These kind of questions are known as the inverse problem of variational calculus (with symmetries in your case), see for example en.wikipedia.org/wiki/Inverse_problem_for_Lagrangian_mechanics for the Lagrangian viewpoint. Depending on your concrete problem the Helmholtz theory can give you an explicit way to integrate your ODE into a Lagrangian (and hence Hamiltonian) language. –  Tobias Diez Jul 7 at 10:57