Automatic inverse for Lambert’s Azimuthal Equal-Area Projection

Question

I’m interested in finding closed-form inverses for map projection equations in general (of course, this is not always possible). I’m hoping Mathematica can help me with some of the trickier ones.

As a practice run, I thought I would try and use Mathematica to invert Lambert’s Azimuthal Equal-Area Projection. There are of course closed-form inverse equations on MathWorld, but I’m wondering how I would go about finding these using Mathematica.

I’ve tried using Solve and Reduce, for example:

k = Sqrt[2/(1 + s1 Sin[y] + c1 Cos[y] Cos[x])]
Solve[{u, v} == {
  k Cos[y] Sin[x], 
  k (c1 Sin[y] - s1 Cos[y] Cos[x])
}, {x, y}]

But no luck so far. Perhaps I need to specify more constraints?

For example, I know that the longitude and latitude (x and y above) have:

-π ≤ x ≤ π
-π/2 ≤ y ≤ π/2

I wonder if eliminating variables might also help the solving routines.