Moon phase data correct only for waxing crescent/gibbous

[aendie]

I have a problem trying to get Moon Phase information from Skyfield 1.11.
I want to replicate the results I get correctly from PyEphem using this (simplified) code:

def moonphase(date):
    # returns the moon's elongation (angle to the sun)
    obs = ephem.Observer()
    obs.date = date
    m = ephem.Moon()
    m.compute(date+0.5)
    #phase = int(round(m.elong.norm*180/math.pi))	# moon phase (0:new to 180:full to 360:new degrees)
    phase = m.elong.norm+0.0    # moon phase as float (0:new to π:full to 2π:new)
    return phase

PyEphem gives me values from 0 to π in the WAXING phase and from π to 2π in the WANING phase. However Skyfield gives me only values from 0 to π in both WAXING and WANING phases. Attached here are the results (in degrees) I get from Skyfield for selected dates while calculating a Nautical Almanac for 2019:
phase_angle.zip

I am using this kind of code in Skyfield:

eph = load('de421.bsp')	# ephemeris
def moonphase(d):
    # returns the moon's elongation (angle to the sun)
    t12 = ts.utc(d.year, d.month, d.day, 12, 0, 0)    # moon phase is calculated at noon
    phase_angle = almanac.phase_angle(eph, 'moon', t12)
    elong = phase_angle.radians
    if elong >= 0 and elong <= math.pi:
        phase = math.pi - phase_angle.radians
    if elong > math.pi:
        phase = (3 * math.pi) - phase_angle.radians   # ??????
    print d, phase_angle.degrees       # FOR DEBUGGING ONLY
    return phase

For example, the 'phase' results are correct between NEW MOON on 6th Jan 2019 up to FULL MOON on 21 Jan 2019. During this time phase_angle.degrees goes from 180 to zero (assumed correct) but then till the next NEW MOON on 4th. Feb 2019 it goes from zero to 180 again (assumed incorrect)... I expected it to go from 360 to 180 degrees.

In other words, the Skyfield value range is 0 to π instead of 0 to 2π.

Do I understand correctly that phase_angle should be able to return a value in radians between 0 and 2π?

Kind thanks in appreciation of your time & effort to look into this!

[ghost]

Technically, the angle between two vectors, as given by the arc-cosine of their dot product divided by their norms, can only be between 0 and pi. I would argue that pyephem's result, while more helpful, is technically incorrect. However, it would be nice to include a flag re whether the angle is increasing or decreasing.

[aendie]

Yes, definitely (include a flag). As mentioned, I am now trying to use Skyfield to create Nautical Almanacs, which can include a graphic of the moon phase. So the "technically correct" answer doesn't help me. I would otherwise have to employ

t, y = almanac.find_discrete(t0, t1, almanac.moon_phases(e))

to determine at least when Full Moon occurs, and then create the phase angle (0 to 2π) that builds the correct moon graphic.

It sounds a little heavy-handed for what should be a simple operation. PyEphem was good at a lot of things - but we have to migrate away from it. Another point ... compared with PyEphem, you can take a coffee break while Skyfield is calculating. (An almanac typically requires time accuracy to hours and minutes, e.g. for transits - the exception is Equation of Time that is given in minutes and seconds.)

[brandon-rhodes]

but we have to migrate away from it.

Note that there is no imperative to migrate away from PyEphem; it is merely not going to be receiving any big new features. Thanks to new tools that have recently been made available to the Python community for using remote machines to compile and release Mac and Windows versions of packages, I have even been able to make a recent release. I’ve updated the home page to sound a little less dire and to make clear the difference between “receiving no new features” and “abandoned”:

https://rhodesmill.org/pyephem/

Hopefully that will offer more clarity!

I am now trying to use Skyfield to create Nautical Almanacs, which can include a graphic of the moon phase.

Moon phase, happily, is not defined by elongation. The Naval Observatory https://aa.usno.navy.mil/faq/docs/moon_phases.php clarifies that, “the phases New Moon, First Quarter, Full Moon, and Last Quarter are defined to occur when the excess of the apparent ecliptic (celestial) longitude of the Moon over that of the Sun is 0, 90, 180, and 270 degrees, respectively.” So you will want to do something like (I'm adapting this from Skyfield's almanac.py):

        _, mlon, _ = e.observe(moon).apparent().ecliptic_latlon('date')
        _, slon, _ = e.observe(sun).apparent().ecliptic_latlon('date')
        angle = mlon.radians - slon.radians

That angle should be in the range [0,2𝜋] and let you draw the crescent Moon facing either left or right. I had to write similar code long ago to illustrate calendar days in the Lord of the Rings using the correctly shaped crescent moon:

https://shire-reckoning.com/moon.html

[aendie]

Hi Brandon! I’m delighted and honoured to hear from you. I’d like to share with you what I’ve done so far… attached are the PyEphem and Skyfield versions of the code I originally took from the version that Enno Rodegerdts originally wrote: https://github.com/rodegerdts/Pyalmanac Attached are the Python 2 versions. Both produce identically formatted output, so one can easily see the differences visually by flipping between the tabs in a PDF viewer. I have not started on checking the accuracy yet, but I thought that Skyfield must produce better results. By the way, Skyfield is slow because two methods take up around 80% of the processing time: ‘daylength’ and ‘moonday’ in alma_skyfield.py. No wonder, really. Also attached are sample PDFs that demonstrate days that have two moonsets or moonrise on one day (in ‘traditional’ and ‘modern’ format style). I only started learning Python at the beginning of the year … so don’t expect brilliant code. I’ve mastered the tricky bits of LaTeX to get just the optical effects I want. I guess I made a couple of mistakes, however I also cured quite a few mistakes that Enno made, listed here: https://github.com/aendie/Pyalmanac I have all the moon phases by now … probably based incorrectly on elongation. They look right to me on initial inspection though. I don’t want to drag you away from your work … just happy when you may have some time to take a look at this! Also hoping that it kind of “advertises” Skyfield for what it can do. (I only cheated in grabbing the planet magnitudes from PyEphem.) Best regards from Munich, Germany! Andrew [Pyalmanac-v3h-python2.zip](https://github.com/skyfielders/python-skyfield/files/3548300/Pyalmanac-v3h-python2.zip) [SFalmanac-v2b-Skyfield.zip](https://github.com/skyfielders/python-skyfield/files/3548302/SFalmanac-v2b-Skyfield.zip) [sample almanac output.zip](https://github.com/skyfielders/python-skyfield/files/3548309/sample.almanac.output.zip)

[brandon-rhodes]

I’m glad you were able to correct several of the mistakes of the earlier programmer — I also found that it took work to get units to display correctly with the final digit rounded, which took work in Skyfield’s timelib module to do correctly.

With the recent release I did of PyEphem, you might find that the navigational stars are now included in its main star database?

Have you ever looked at ReportLab? I have sometimes found it easier to use than LaTeX for simple tables.

Are you also going to publish the output somewhere online, or is this only to be expected to be used by people who can install the software?

Finally, let me know if this issue can be closed, or if there's still an angle Skyfield is returning that you think is incorrect. Thanks!

[aendie]

Yes, you can close this thread as I have implemented "before or after Full Moon" as indication of the angle. I didn't think it is incorrect - I thought that Skyfield left out important information in providing an angle that's only between 0 and π instead of between 0 and 2π. So the suggestion (from barrycarter) to add a flag 'angle increasing' or 'angle decreasing' I think is very appropriate.

I haven't looked at ReportLab - I continued using LaTeX that the original developer, Enno Rodegerdts, chose. I have one open issue with pdflatex (included in MiKTeX) that I'll pick up with them directly. But thanks for the suggestion. I'll also look into the navigational stars in the latest PyEphem. Thanks!

The PyEphem output PDFs are already published online here: https://www.thenauticalalmanac.com/
I am in touch with the site owner and my name is mentioned very kindly there too. He does some clever post-processing before publishing. I've just finished the Skyfield version, so I expect it will be published when we're sure the data is okay.

By the way, the moon image was fun to create. I overlaid a semi-transparent shadow onto the moon. First I simply deleted the black pixels from a moon image and the result was hilarious: a "hairy" moon (due to moonshine). I then wrote code to place transparent pixels outside a circle, carefully positioning it over the moon, which is almost a circle. Finally I thought maybe it looks better without the moon, so I included a parameter in config.py to replace it with the original data in that table location.

Thanks for your helpful comments!

P.S. My comment about Skyfield performance... clearly PyEphem uses different algorithms, but do you perhaps see a way to make the following methods execute faster?

def daylength(topos, degBelowHorizon):
    # Build a function of time that returns the daylength.
    topos_at = (earth + topos).at

    def is_sun_up_at(t):
        t._nutation_angles = iau2000b(t.tt)
        # Return `True` if the sun has risen by time `t`.
        return topos_at(t).observe(sun).apparent().altaz()[0].degrees > -degBelowHorizon

    is_sun_up_at.rough_period = 0.5  # twice a day
    return is_sun_up_at

and

def moonday(topos, degBelowHorizon):
    # Build a function of time that returns the "moon day" length.
    topos_at = (earth + topos).at

    def is_moon_up_at(t):
        t._nutation_angles = iau2000b(t.tt)
        # Return `True` if the moon has risen by time `t`.
        return topos_at(t).observe(moon).apparent().altaz()[0].degrees > -degBelowHorizon

    is_moon_up_at.rough_period = 0.5  # twice a day
    return is_moon_up_at

... should I enter this as a new thread?

[aendie]

With the recent release I did of PyEphem, you might find that the navigational stars are now included in its main star database?

Thanks - I already had all the stars I needed for the Nautical Almanac in the correct sequence, with names as used officially in almanacs.

In PyEphem 3.7.7 the following stars are duplicates (with just slightly different data):

• Sirrah + Alpheratz
• Adhara + Adara (mispelled)
• Gienah + Gienah Corvi
• Alcaid + Alkaid
• Hagar + Agena
• Etamin + Eltanin
• Fomalhaut + Formalhaut (mispelled)

I would suggest removing duplicates (keeping the better data) and sort all by RA (as in PyEphem 3.7.6).

[brandon-rhodes]

Yes, you can close this thread…

Thanks, I'll do so. Yes, please open the performance question separately.

Thanks for the links to the PDF output, it's interesting to see the result!

I do see now what you mean about PyEphem — I don't remember looking at its strange nonstandard definition of “elongation” before. Skyfield will stick for now with the Wikipedia definition of a simple angle <180°. Here's what PyEphem does, though:

/* given geocentric ecliptic longitude and latitude, lam and bet, of some object
 * and the longitude of the sun, lsn, find the elongation, el. this is the
 * actual angular separation of the object from the sun, not just the difference
 * in the longitude. the sign, however, IS set simply as a test on longitude
 * such that el will be >0 for an evening object <0 for a morning object.
 * to understand the test for el sign, draw a graph with lam going from 0-2*PI
 *   down the vertical axis, lsn going from 0-2*PI across the hor axis. then
 *   define the diagonal regions bounded by the lines lam=lsn+PI, lam=lsn and
 *   lam=lsn-PI. the "morning" regions are any values to the lower left of the
 *   first line and bounded within the second pair of lines.
 * all angles in radians.
 */
static void
elongation (double lam, double bet, double lsn, double *el)
{
	*el = acos(cos(bet)*cos(lam-lsn));
	if (lam>lsn+PI || (lam>lsn-PI && lam<lsn)) *el = - *el;
}

It looks like it uses a difference in ecliptic longitude as the trigger for changing the sign, so my original recommendation — that Skyfield users compare ecliptic longitudes when interested in comparing a planet’s position with the Sun’s — still stands, and should produce compatible results.

I'll go update PyEphem’s documentation to remove the misleading description of that hybrid value as “elongation” and include the sign-flipping behavior in the description.

[tkircher]

Ran into this exact problem today, expecting skyfield to give me values in either 0 to 2pi or -pi to pi. It's definitely wrong for the phase angle of the first and last quarter to be the same.

[brandon-rhodes]

It's definitely wrong for the phase angle of the first and last quarter to be the same.

@tkircher — Skyfield is following both the United States Astronomical Almanac and (for what it’s worth) the Wikipedia in its definition of “phase angle” as always 0°≤p≤180°. While I won’t go so far as to claim that makes Skyfield “definitely right” as against your claim that it’s “definitely wrong”, Skyfield is certainly “in agreement with the Astronomical Almanac” which is the high standard to which I generally hold it.

If we can find a technical term for “a number running from 0° to nearly 360° as the Moon goes from New Moon to the next New Moon” (which definitely isn’t “phase angle”; I just double-checked the Almanac’s Explanatory Supplement to be sure), I will be happy to add a new function with that name to the almanac module. It will look like the moon_phase_at() function returned by moon_phases(), but instead of reducing the result to an integer, will simply subtract the two angles in degrees and then return the result modulo 360°.

(It’s the same number, I suppose, as we use with other bodies for conjunction and opposition: a comparison of an ecliptic longitude with the Sun’s ecliptic longitude. I feel like there’s a name for that?)

[tkircher]

How about signed_phase_angle to go along with phase_angle, analogous to atan() and atan2().

[brandon-rhodes]

I'll think about that proposal. The problem is that “phase angle” is by definition (for those who look it up; I readily admit that other reasonable definitions might come to mind from simply that pair of words) is the real-world angle between two things, whereas if we want to know the phase of the Moon where that means "new - quarter - full" we really want a difference in ecliptic longitudes that completely ignores whatever real angle is between the two bodies ­— since the Sun and Moon are never really 180° apart, not exactly, yet I think for your purposes you want an angle that unambiguously says "full Moon" without having to do maxima detection ("ah, yes, this month the maxima is 179.1°, I'll look for that angle to mean 'full'!").

What about a routine called moon_phase() that returns an angle 0°–360°? Or do you really need the angle that one month might only reach 179.1°, but you need it to swap to -179.1° at exactly the moment of Full?

[tkircher]

I mean, it doesn't matter to me what it's called, so moon_phase() is great. Thanks again for the response. I understand the desire to maintain correctness relative to the USAA and also appreciate your flexibility in offering to implement this utility function.

[brandon-rhodes]

@tkircher — Thanks for the response! One final clarification; do you want:

1. New Moon 0° → First Quarter 90° → Full 180° → Last Quarter 270° → 359.99° right before the next New Moon,
2. OR do you want something more like a phase angle, so 0.7° New moon → around 90° at First Quarter → 179.1° right before Full / –179.1° right after Full → around –90° at Last Quarter → 0.65° or whatever the angular separation is between the Sun and Moon at the next New Moon.

My guess is the first, since you mentioned 0 and 2𝜋 in your first comment, but I wanted to double-check before adding the routine. 🙂

[tkircher]

Yeah, the first one would actually be great!

[brandon-rhodes]

@tkircher — If you run:

pip install https://github.com/skyfielders/python-skyfield/archive/master.zip

— then your Skyfield should now support almanac.moon_phase(eph, t) that will return the Moon phase as an Angle object. If you get the chance to try it out, let me know if it gives you sensible values!

[tkircher]

@brandon-rhodes gave it a test and it works great. Thanks!