$
\cos ^{-1}\left(-\frac{\tan (\lambda ) (\cos (\lambda ) \cos (\phi ) \cos
(Z)+\sin (\lambda ) \sin (Z))}{\sqrt{(\cos (\lambda ) \sin (Z)-\sin (\lambda)
\cos (\phi ) \cos (Z))^2+\sin ^2(\phi ) \cos ^2(Z)}}\right)-\tan^{-1}
(\cos (\lambda ) \sin (Z)-\sin (\lambda ) \cos (\phi ) \cos (Z),\sin
(\phi ) (-\cos (Z)))
$
is the amount of time from now (see important notes below) a celestial object will set, where:
- $\phi$ is the azimuth of the object
- $Z$ is the altitude of the object above the horizon
- $\lambda$ is the latitude of the observer
Note the two-argument form of arctangent is required for accuracy: https://en.wikipedia.org/wiki/Inverse_trigonometric_functions#Two-argument_variant_of_arctangent
Important notes (MUST READ!):
The result is in radians. To convert to sidereal hours, multiple by $\frac{12}{\pi }$
If the object has a fixed right ascension and declination, divide sidereal hours by 1.002737909350795 to get clock hours.
For the Sun (which doesn't have fixed right ascension/declination), do NOT divide as above. By not dividing, you compensate for the Sun's change in right ascension.
For the Moon, see "The Moon" section below.
The time computed above is for geometric midpoint setting. In reality, refraction near the horizon means the object will set later. Additionally, for objects that have an angular diameter (eg, the Sun), setting usually means when the top edge disappears over the horizon, which will make the set time even later. Accounting for both effects should be possible, but might require numerical methods instead of a closed form formula as above.
If the quantity inside the arccosine is greater than 1, the object is always in the sky and never sets or rises.
If the quantity inside the arccosine is less than -1, the object is never in the sky and thus also never sets or rises.
I did only minimal testing. As always, do not rely on my answers for anything important.
The Moon:
The Moon's right ascension and declination change rapidly, so this calculation does not work well for the moon.
You could compensate for the change in right ascension (and thus hour angle) by approximating the increase as 24 hours every 27.32158 days (its sidereal period) and do an iterative calculation.
An even better compensation for the change in right ascension would be to approximate the moon's movement in ecliptic longitude (which is more constant than its movement its right ascension) as 360 degrees per 27.32158 days and then project the ecliptic longitude back to right ascension, and then iterate.
Compensating for the moon's change in declination is more difficult. The moon's ecliptic latitude (which can be converted to declination) varies sinusoidally, but the equation $\sin (x)=a$ normally has two solutions. Unless you know whether the moon's ecliptic latitude is increasing or decreasing (ie, whether it's between ascending and descending nodes or vice versa), you won't know the direction of declination change.
Less important notes (optional):
Solution notes:
I learned quite a bit answering this question, and thought it was only solvable numerically until I figured out the shortcut:
To convert from the azimuth/elevation sphere to the hour angle/declination sphere, you just rotate around the y axis by $\frac{\pi }{2}-\lambda$ (90 degrees minus the latitude) and then rotate $pi$ (180 degrees) around the z axis.
Once you have the declination, computing the hour angle when an object sets is easy.
You then subtract the setting time from the current hour angle to get the answer.
TODO: I was hoping to add some graphs/visualizations to this answer but couldn't quite get things working.