Express the involved polynomials in cos(T) and sin(T) in terms of cos(2T) and sin(2T); you'll get a linear equation in those, which can be easily solved e.g. by converting it further using a "helper" angle.

As to "the angle that comes closest", what exactly is your measure of closeness? The answer depends on this a lot, but in any case, you can calculate the locus of reachable point pretty easily and proceed from there.

as Darkstrike wrote, you need to transform the sin & cos terms using
sin(2x)=2*sin(x)*cos(x) and
cos(x/2)^2=(1+cos(x))/2

When I calculated it, I used a different direction of the acceleration, but it's just to insert (-a) instead of a.

s^2*h*cos(T)^2 = s^2*sin(T)*cos(T)*d - a*d^2/2

(1+cos(2T))*h/2 = sin(2T)*d/2 - a*d^2/(2*s^2)

(a*d^2+h*s^2)/s^2 + h*cos(2T) = d*sin(2T)

Set (a*d^2+h*s^2)/s^2 = C.

h*cos(2T) + C = d*sin(2T)

cos(2T)*h/d + C/d = sin(2T)

For the rest I made a very nice Paintbrush picture ^^



T = (asin(C/sqrt(d^2+h^2)) + atan(h/d))/2

I could have forgot some special case for asin or atan, but the basics should be correct. Probably a couple of typos too... =)

Sadly, what you described neglects air resistance/wind -- which are rather important, if you want realistic projectile motion.

There is an alternative way to solve such problems. It isn't always best.

Basically, you search for the right firing solution, instead of calculating it.

Start out with sampling the "space" of possible firing options either randomly, or with some best guesses.

Now use those samples to figure out a better set of samples. Repeat until you are close enough.

If you have a simple enough problem, you can often get solutions that 'braket' your target -- one shot is to the left, the other to the right. In that case, a simple binary or dictionary search often rapidly produces an accurate solution.

..

In your particular example, let's first solve for t in terms of h and Vy.

a/2 t^2 + Vy t - h = 0

t = -Vy +/- sqrt(Vy^2 + 2 a h)
----------------------------
a

There are two solutions to the above. So for a given Vy, we have t0 and t1.

Next, see how far away horizontally the shot is at t0 and t1:

d = Vx * t
d0 = t0 Vx
d1 = t1 Vx

That is pretty easy!

Now, we start out by testing a bunch of different random firing solutions. We stop after testing a large number, or after we have 'bracketed' the target (got one to land on the far side, and one to land on the near side, of the target).

Once we have 'bracketed' the target, we simply engage in a binary search through the firing space between the two options.

If we haven't bracketed the target, it is possible the target cannot be hit. We take the closest shot to the target in t0, search the solutions near it, and give up if we stop getting closer. We then do the same with t1 solutions.

All of the above does is solve the problem in a way that doesn't rely on the equations described being 'easy to solve' -- if you added tweaks to the equations (wind resistance), the above solution would be more robust (it wouldn't be guaranteed to work!) compared to the analytical solution.