I am trying to calculate the Spatial Delaunay Elements with actions $(L, G, H)$ and angles $(l, g, h)$ from polar coordinate $\left( \rho, \theta, \phi\right)$ and momentum $\left( P, \Theta, \Phi\right)$, where $x=\rho\sin\phi\cos\theta$, $y=\rho\sin\phi\sin\theta$, $ z=\rho\cos\phi$. The Hamiltonian $$\mathcal{H} = \frac{1}{2}\left( P^2+\frac{\Phi^2}{\rho^2} + \frac{\Theta^2}{\rho^2\sin^2\phi}\right)-\frac{1}{\rho}$$
I have $L = \frac{1}{\sqrt{-2\mathcal{H}}}$, $G = \mathcal{G} = L\sqrt{1-e^2}$, and $H = \Theta=G\cos i$.
In order to get the angles $(l, g, h)$, I set the generating function $$S(I, q)=\int_{\rho_0}^\rho \sqrt{-\frac{1}{L^2}-\frac{G^2}{\rho^2}+\frac{2}{\rho}}\mathrm \rho + \int_{\theta_0}^\theta H\mathrm d \theta + \int_{\phi_0}^\phi \sqrt{G^2-\frac{H^2}{\sin^2\phi}}\mathrm d \phi$$ where $\rho_0 = \rho_{\text{min}}, \theta_0 = 0, \phi_0 = \frac{\pi}{2}$
so $$h = \frac{\partial S}{\partial H}$$
However, the problem is that the first integral gives $$h_\theta = \int_{\theta_0}^\theta \frac{\partial H}{\partial H}\mathrm d \theta=\theta$$ while, using the relation $\sin i\sin(\theta-\Omega) =\cot\phi\cos i$ ($i, \Omega$ are Keplerian Elements), the second gives $$h_\phi=\int_{\phi_0}^\phi \frac{\partial}{\partial H}\sqrt{G^2-\frac{H^2}{\sin^2\phi}}\mathrm d \phi=\int_{\Omega}^{\theta}\mathrm d\theta=\theta-\Omega$$ which could be verified by Mathematica
FullSimplify[D[Sqrt[G^2 - H^2/Sin[\[Phi]]^2], H] D[ArcCot[Tan[i] Sin[\[Theta] - \[CapitalOmega]]], \[Theta]] /. {\[Phi] -> ArcCot[Tan[i] Sin[\[Theta] - \[CapitalOmega]]]} /. {H -> G Cos[i]}, Assumptions -> {0 < i < Pi, G > 0}]
(It yields Sqrt[Cos[\[Theta] - \[CapitalOmega]]^2]Sec[\[Theta] - \[CapitalOmega]] $=1$)
And the two integrals add to $2\theta-\Omega$, which should be $\Omega$
I wish anyone could kindly point out my mistake. thanks