Lunar Magnitude

[mkbrewer]

Hello Brandon, It's been awhile since I contacted you. I have a few issues with PyEphem that I'd like you to address when you have the time. First of all, and this has been a peeve of mine for a long time, but I never brought it up with you before, there are issues with the diameter of the Moon and its elongation from the Sun. The diameter is off a bit, so I always use:

AU = 149597870.7 #km
MEAN_MOON_RADIUS = 1737.4 #km
radius = 3600.0 * np.degrees(MEAN_MOON_RADIUS / (moon.earth_distance * AU)) #arc seconds

which yields the correct value. Perhaps you are using a different value for the mean radius. Mine comes from the NASA fact sheet. https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html I use the volumetric value.

Second, for some inexplicable reason, the Moon's elongation from the Sun is always geocentric. I need to compute the correct topocentric value on my own. This also affects the phase angle and illumination fraction.

#Moon's` topocentric elongation from the Sun.
elong = np.arccos(np.sin(sun.dec) * np.sin(moon.dec) + np.cos(sun.dec) * np.cos(moon.dec) * np.cos(moon.ra - sun.ra))
#Moon's topocentric phase angle (Sun-Moon-Observer).
phase = np.arctan2(sun.earth_distance * np.sin(elong), moon.earth_distance - sun.earth_distance * np.cos(elong))
illum_frac = 0.5 * (1.0 + np.cos(phase))

Now that that's over with, I recently responded to a question on the Astropy mailing list about computing the lunar magnitude and surface brightness by recommending either PyEphem or Skyfield. However, when I looked into this further, I found that Skyfield doesn't actually provide lunar magnitude, only that of the planets. I also found that the lunar magnitude provided by PyEphem and JPL Horizons differ by up to 3 magnitudes!

In order to get to the bottom of this, I used the widely cited Lane, Adair P.; Irvine, William M. 1973. Monochromatic phase curves and albedos for the lunar disk. Astronomical Journal, Vol. 78, pp. 267-277 to model the magnitude of the Moon in the visual band. Following their prescription, I fit the V band column in their Table V to a line out to 40 degrees phase angle and to a cubic from 40 to 120 degrees. This yielded:

if np.degrees(phase) < 40:
    p = np.degrees(phase) - 20.0
    magnitude = -12.72 + 0.0267 * p + 0.534
else:
    p = np.degrees(phase) - 80.0
    magnitude = -12.72 + p * (0.03188 + p * (1.9621e-4 + p * 1.7556e-6)) + 2.14

Since their numbers are given for mean Moon-Earth distance and a Moon-Sun distance of 1 AU, I also added these corrections:

MEAN_MOON_ANGULAR_DIAMETER = 1896.0 #arc seconds
magnitude += 2.5 * np.log10((MEAN_MOON_ANGULAR_DIAMETER / (2.0 * radius)) ** 2)
moon_sun_distance = sun.earth_distance - moon.earth_distance #AU
magnitude -= 2.5 * np.log10((1.0 / moon_sun_distance) ** 2)

and finally:

surface_brightness = magnitude + 2.5 * np.log10(np.pi * illum_frac * radius ** 2)

I compared the results of this model for two two month periods sampled at 12 hour intervals centered on the January 2nd, 2024 Earth perihelion and the July 5th, 2024 Earth aphelion with those of Horizons. The average difference for both periods was zero and the standard deviation was 0.07 magnitude. The difference is generally less than 0.1 magnitude except very near the new moon where it approaches 0.2 magnitude. The surface brightness equation appears to be exact in that it produces the same difference as that of the magnitude to 3 decimal places.

My conclusion is that there is something seriously wrong with PyEphem's calculation of lunar magnitude.

[mkbrewer]

I noticed a couple of slight errors in the above. The cubic equation should have read:

magnitude = -12.72 + p * (0.03188 + p * (1.9621e-4 + p * 1.7256e-6)) + 2.14

This was a transcription error on my part. Also, I missed one term in the equation for the Moon-Sun distance. That should have been:

moon_sun_distance = sun.earth_distance - moon.earth_distance * np.cos(elong) #AU

These corrections make little practical difference, but they do result in a slight improvement in the agreement with Horizons (0.01 - 0.03 magnitude depending on phase).

[mkbrewer]

By the way, the person who asked about this on the Astropy mailing list stated that their intended use of the lunar magnitude and surface brightness was to study the effect of the Moon on the behavior of marine mammals over large geographic extents and decadal time scales. They said that querying JPL Horizons was much too slow, which is why I recommended PyEphem. It is interesting to me to see how widely varied the applications of an astronomical ephemeris program can be.

[drbitboy]

Using JPL/NAIF SPICE is the "right" was to get geometry programmatically from JPL.

Cf. https://naif.jpl.nasa.gov/

And SpiceyPy is the Python interface to SPICE.

[brandon-rhodes]

@drbitboy — Thank you for pointing out that that folks can get HORIZONS-compatible results produced locally from their own computer, without needing to interact with the online service.

Do you know what calls someone would make who wanted to produce, from SPICE, the same value for lunar magnitude that HORIZONS produces?

@mkbrewer — So, the Moon’s angular size is, as you suggest, not based on its mean radius. Instead, XEphem returns its largest diameter by using its equatorial radius:

#define MRAD    (1.740e6)   /* moon equitorial radius, m */

I’ll update the docs to mention this. I think it’s a reasonable decision for an amateur astronomy library, since if you want to photograph the Moon and catch the whole disc, then it’s the maximum diameter you need, not the average.

I’ll also update the docs to make it clearer that elongation is geocentric.

The Moon magnitude does look to be a known problem. The code says:

    /* TODO: improve mag; this is based on a flat moon model. */
    i = -12.7 + 2.5*(log10(PI) - log10(PI/2*(1+1.e-6-cos(el)))) 
                    + 5*log10(edistau/.0025) /* dist */;
    set_smag (op, i);

So I would look favorably on a pull request that replaced that line of code with a more accurate formula — and, please go ahead and remove the TODO while you’re at it. 🙂

[mkbrewer]

I think it’s a reasonable decision for an amateur astronomy library, since if you want to photograph the Moon and catch the whole disc, then it’s the maximum diameter you need, not the average.

I'll grant you that, but I'll stick with my definition. Coming from the field of radio astronomy, I tend to be more interested in the Moon's solid angle than its maximum diameter. We use this to calculate beam dilution for beams larger than the Moon.

So I would look favorably on a pull request that replaced that line of code with a more accurate formula

Thanks! I'll work on this.

[drbitboy]

Coming from the field of radio astronomy, ...

Then you should convert to SPICE; the learning curve is steep, but worth it.

Do you know what calls someone would make who wanted to produce, from SPICE, the same value for lunar magnitude that HORIZONS produces?

SPKEZR, VNORM, VADD and such. It's just vector math. We would still need the formula for magnitude, but I assume that is available.

Also, there is a telnet-ish interface to Horizons (port 6775, iirc - yeah, see here), or there used to be. I once wrote some eXcel macros to access it, so it would not be impossible to do it programmatically, especially with Python.