NavList:

Tibor Miseta, you wrote:
"I've tried to estimate by the illuminated disk."

Nice! In practice this is very difficult for any useful level of accuracy, but the concept "illuminates" an important aspect of determining time by lunar distance angles which I sometimes talk about in my Lunars workshops. 

Imagine yourself in any early human society on any continent 5-10 thousand years ago... You could have made plans by the Moon. The phase of the Moon, meaning the apparent fraction of the illumination of the Moon's disk, has been used as a calendar reference "since the dawn of history" (which is to say that just about every society that we know of historically in fact kept a lunar calendar, either by itself or in addition to a solar calendar.

Suppose I lived here in southern New England 5000 years ago... I've got some salted fish to trade. A friend from upriver expects to have a bag of beans in about three weeks. We know the Moon was near full last night, so I suggest to my friend, "When the Moon is exactly half full after the New Moon, I will meet you at the second fork in the river in the afternoon". We agree that the location is unambiguous (even without a "what three words" address). The specific day when the Moon is just half full is easily determined by watching it, and better observers may also have realized that half full occurs when the Moon is at a right angle from the Sun. We would have no problem getting the date right and making our trade as planned. It works. And we have made our appointment using lunar distances for a rough estimate of the time (really just the date in this case but that's only a question of accuracy)!

Wait. Lunar distances? Sure. You noted that there is a simple equation, derived by basic geometry, for the illuminated fraction of the Moon. It can be written as 
  f = (1 - cos LD)/2,
where LD is the "lunar distance" angle from the Sun to the Moon. Estimating the illuminated fraction of the Moon is equivalent to estimating this simple function of the standard Sun-Moon lunar distance. And fans of haversines, take note! That right there is the definition of the haversine so we could re-write this as
  f = hav LD.
There's no practical benefit re-writing it this way, but it may entertain you, and it provides a visual graph of the haversine function (the phase of the Moon over a month ranging from f=0 at LD=0° to 0.5 at LD=90° and then 1.0 at LD=180° and so on exactly demonstrates the values of hav x over its full range).

You also wrote:
"A better estimation if we try to calculate the illuminated disk from GHA differences. By one of Meeus' approximate formula the illuminated fraction is k = (1-cos i)/2, where cos i is approx - cos(delta Lat) * cos (delta Lon), so an estimate is k = (1 - cos(delta GHA))/2. (I neglected the latitude differences, because they are so close, that the cosine of small angles are nearly 1.) "

There may be an issue here with the meaning of Lat and Lon. Suppose we want to calculate the Sun-Moon lunar distance angle from the positions of the Sun and Moon. This is just a great circle distance problem. In any coordinates with a latitude-like coordinate measured from an "equator" and a longitude-like coordinate measured around a central axis perpendicular to that equator, we know that
  cos(LD) = sin(L1) sin(L2) + cos(L1) cos(L2) cos(dLon),
where L1 and L2 are the latitudes in the selected coordinates and dLon is the difference in longitude in those coordinates. This always works, and we can just drop in appropriate values, for example the Dec and GHA of the Sun and Moon taken directly from an almanac or other database. Historically (as recently as twenty years ago), it was not unusual to generate Moon and Sun coordinates from algorithms that started in ecliptic coordinates and then converted to equatorial coordinates (like Dec and GHA). If ecliptic coordinates are available, then there's a nice simplification that occurs since the ecliptic latitude of the Sun is 0°. In the equation for cos(LD) here, we can call object "1" the Sun. Then L1 is 0°. The first term is necessarily zero and cos(L1) is exactly one. What we're left with then is
  cos(LD) = cos(L2) cos(dLon),
and then the fraction illuminated is 
  f = [1 - cos(L2)cos(dLon)]/2.
That's a nice simple result (which, I emphasize, is specific to ecliptic coordinates), and if we have an ephemeris giving the Moon's ecliptic lat and lon and the Sun's ecliptic lon, we can then generate a table of the Moon's fractional illumination as a function of time.

This is all very nice in principle, but in practice there are big problems. It certainly works to determine the time to the nearest day or so, and my alter-ego 5000 years ago would have traded salted fish for beans on the afternoon of the Half Moon --no problem. But the percent illumination changes too slowly to function as any practical check on the hour of UT/GMT. At its fastest rate of change, near First and Last Quarter when the percent illumination is close to 50%, the rate of change is only about half a percent per hour. In a photo like the one we've been discussing I can't imagine estimating the percent illumination to better than the nearest percent which means we could only estimate time to the nearest two hours so so. On the date in question, when the Moon was a slim crescent, the rate of change in percent illumination was even lower, and determining the illumination to the nearest single percent would only yield the absolute time to the nearest four to six hours. 

Of course, what we require is a better way of measuring the fractional illumination of the Moon, or anything equivalent to it. That lunar distance angle is a simple equivalent to the fractional illumination so we need an angle-measuring device. We invent the reflecting quadrant (or octant), improve it as a sextant, and the rest is history. And yet all we're really doing is estimating the phase of the Moon, like our ancestors did thousands of years ago, when we shoot lunars. 

Frank Reed