Plotting parametric curve with implicit parameter function

Question

I need to plot a Smale diagram as a parametric curve $(p(x), h(x))$ where $h$ is amended potential and $p = \omega^2$ is parameter of bifurcation: $$h:=V_\omega(x,y)=-2\cos(x)-\cos(y)-\frac{1}{2}\omega^2\bigl(\sin^2(x)+(\sin(x)+\sin(y))^2\bigr)$$ There is system of equilibrium equations (partial derivatives of amended potential are zero): $$h_x'=2\sin(x)-\omega^2\cos(x)(2\sin(x)+\sin(y))=0$$ $$h_y'=\sin(y)-\omega^2\cos(y)(\sin(x)+\sin(y))=0$$

I am trying to get $y = y(x)$ and $p = p(x) = p(x, y(x))$ from this system. For example, isolating $\sin(y)$ and $\cos(y)$, and using $\sin^2(y)+\cos^2(y)=1$, it can be turned into: $$4\bigl(1-\omega^2\cos(x)\bigr)^2\Bigl(\cos^2(x)+\bigl(2-\omega^2\cos(x)\bigr)^2\sin^2(x)\Bigr)=\omega^4\cos^2(x)(2-\omega^2\cos(x))^2$$ Then, substituting $\cos(y)$ and $\sin(y)$ in $h$: $$h=-2\cos(x)-\frac{2\bigl(1-\omega^2\cos(x)\bigr)}{\omega^2\bigl(2-\omega^2\cos(x)\bigr)}-\frac{1}{2}\omega^2\sin^2(x)\biggl(1+\Bigl(1+\frac{2\bigl(1-\omega^2\cos(x)\bigr)}{\omega^2\cos(x)}\Bigr)^2\biggr)$$

So I got something like parametrization for $h$ and $p$ with the parameter $x$, but for $p(x)$ it is implicit. I had attempted to plot this in Maple 13, but unsuccessfully.

(it's described here)

What I want to know - is there a way to plot this parametric curve on $(h,\omega^2)$ coordinates with computational methods, even numerically (Wolfram Mathematica or Maple would be the best)?

Particularly, for one curve $h_-(\omega^2)$, separating from point $(x,y,\omega^2)=(0,0,2-\sqrt{2})$ with the following conditions: $$\omega^2>2-\sqrt{2}$$ $$h(2-\sqrt{2})=-3$$ $$\sin(x)\sin(y)>0, (x,y)\in O_\varepsilon (0)$$

Answer 1

With MATHEMATICA the leaves's graphics for $p(x) = \omega^2(x)$

ContourPlot[4 (1 - p Cos[x])^2 (Cos[x]^2 + (2 - p Cos[x])^2) Sin[x]^2 - p^2 Cos[x]^2 (2 - p Cos[x])^2 == 0, {x, -Pi, Pi}, {p, 0, 2 Pi}, ContourStyle -> {Thick, Blue}]

Follows a script to plot the parametric $(p(x),h(x))$ for one of $p(x)$ leaves,

Clear[p, equ, h]
equ[p_, x_] := 4 (1 - p Cos[x])^2 (Cos[x]^2 + (2 - p Cos[x])^2) Sin[x]^2 - p^2 Cos[x]^2 (2 - p Cos[x])^2
h[p_, x_] := -2 Cos[x] - 2 (1 - p Cos[x])/(p (2 - p Cos[x])) - 1/2 p Sin[x]^2 (1 + (1 + 2 (1 - p Cos[x])/(p Cos[x]))^2);
x = -Pi;
Sols = {};
dx = 0.01;

While[x < Pi,
  sol = Solve[{equ[p, x] == 0, p > 0}, p, Reals];
  If[sol != {}, AppendTo[Sols, {x, Flatten[p /. sol][[1]]}]];
  x += dx
]

par = {};
For[k = 1, k <= Length[Sols], k++, 
  {xx, pp} = Sols[[k]]; 
  If[Max[pp, Abs[h[pp, xx]]] < 10, AppendTo[par, {pp, h[pp, xx]}]]
]
ListLinePlot[par, PlotStyle -> {Thick, Blue}]

Answer 2

From $h'_y=0$, you have $$\sin (x)+\sin (y)=\frac{\tan (y)}{\omega ^2}$$ Replacing in $h'_x=0$, then $$y=\tan ^{-1}\left(2 \tan (x)-\omega ^2 \sin (x)\right)$$