Calculating the illumination of the moon I am trying to find a formula that will enable me to calculate the illumination of the moon down to one thousandth of a percent, given that the Gregorian year, month, day, and hour is known. Can anyone help me with this formula?
 A: I'm unsure what is meant with

[...] calculate the illumination of the moon down to one thousandth of a percent [...]

Assuming that you are trying to determine if a certain point on the lunar surface is illuminated or not, I could possibly give an approach for that - I've done that within my Master's Thesis, since one of my tasks is the development of a Moon surface illumination simulation software. (If you are interested in this topic, first results can be seen here; the thesis itself will be published in spring 2012.)
The required astrodynamical calculations for this purpose are quite complex; I don't believe that a single equation can be derived. For high accuracy astrodynamical calculations an extensive set of differential equations needs to be considered and evaluated - according to the level of accuracy which is needed.
Hence, I would advise to use an external library like the NASA NAIF SPICE toolkit, which is available for many environments (C, C++, Matlab, IDL, Fortran). This way, you have to write a short program in one of those languages.
Use SPICE to calculate the Sun's position at a certain point of time (your Gregorian date is the input) w.r.t. the Moon's Mean Earth/Polar Axis reference frame (ME/PA). This reference frame is defined bodycentric for the Moon. Its origin is the lunar center, z-axis is the polar axis, x-axis is the mean Earth direction. This way, the location (0° N, 0° E) on the Moon's surface, expressed by the pair (lat, lon) (latitude and longitude), will be at $(r, 0, 0)$ in the ME/PA reference frame, whereas $r$ is the lunar volumetric mean radius. More examples: $(0^\circ\,\text{N}, 90^\circ\,\text{E}) \mapsto (0,r,0)$, $(90^\circ\,\text{N}, 0^\circ\,\text{E}) \mapsto (0,0,r)$, $(0^\circ\,\text{N}, 90^\circ\,\text{W}) \mapsto (0,-r,0)$.
Knowing the Sun's position in this frame, you are able to determine whether a sunray will hit the given surface point or not by using a simple ray tracing algorithm.
Note: This neglects possible shadow occultations caused by the lunar topography. If such considerations are necessary in your case, there is likely no simple solution, as you must maintain an entire lunar digital elevation model.
A: I let the programmers of my planetarium software do all the heavy lifting for me with questions like this. I use Starry Night myself, and have come to trust the answers it gives me. If you insist on doing it yourself, Jean Meeus' books are the place to go.
