Let's say we have a 1D harmonic oscillator, its Hamiltonian is given by
$$H=\frac{p^2}{2m}+\frac{1}{2}m\omega^2x^2$$
we wish to solve it via the Hamilton-Jacobi equation so we have
$$\frac{1}{2m}\Big(\frac{\partial S}{\partial x}\Big)^2+\frac{1}{2}m\omega^2x^2+\frac{\partial S}{\partial t}=0$$
since the Hamiltonian doesn't explicitly depend on time we can write $$S(x,p_0,t)=W(x,p_0)-Et$$ so that
$$\frac{1}{2m}\Big(\frac{\partial S}{\partial x}\Big)^2+\frac{1}{2}m\omega^2x^2=E$$
and we can solve for $W$ a well just by simple integration and get S
$$S=\sqrt {mk}\int\sqrt{2E/k-x^2}dx-Et$$
if we solve the integral we get some function of x
what I don't get is the following: In order to obtain the solution for $x(t)$ what we do is derive $S$ by the energy?
$$\beta=\frac{\partial S(x,t,E)}{\partial E}=f(x,t,E) $$ and then we invert $f(x,t,E)$ to get $x(t)$
But what I don't get is how we can differentiate wrt. the energy because it is not like we are looking at how does our solution change when the energy of our system changes.
Also if we had a 2D system and wanted to find $y(x)$ or $x(y)$ that is how one variable changes with respect to the other it is common to derive by $p_0$
$$\beta_2=\frac{\partial S(x,t,E)}{\partial p_{0x}}=g(x,y,E)$$
(because in the 2-D case we would have $p_{0x}$ and $p_{0y}$this is even more puzzling to me, how can we derive by $p_0$ if we know it to be a constant for that system, in deriving the Hamilton-Jacobi equation we wanted to find a new coordinate system for which the new Hamiltonian is zero so that the equations of motion become trivial $\dot{q}=0$ and $\dot{p}=0$ ($q=q_0$ , $p=p_0$) so if we know p_0 is a constant, but we none the less derive by it? how and why?