Calculating moon face rotation as a function of Earth coordinates

Question

I'm writing an Android app that shows moon phases for any date. Currently I am showing a generalized view from the northern hemisphere (sunlit part of the moon moves from right to left with 0º viewing rotation). However, this motion is reversed when viewing the moon from the southern hemisphere, and near the equator, something else happens entirely. Some of this information can be found on this webpage on moon phases.

Instead of using just three simplified cases (northern, southern, and equator), I would ideally like to be able to calculate the exact angle of moon viewing rotation as a function of a latitude and longitude on Earth:

double viewingRotationAngle = calculateAngle(earthLat, earthLng);

Does anyone know how to do this math, or where I could find more information to potentially figure it out myself? I have searched for several days and haven't yet found what I'm looking for. Thanks in advance.

Just to clarify, the angle of rotation I seek is the 2D rotation of the visible disk of the moon.

Update 3/17/14

I did more research and realized this questions may have more to do with the celestial sphere in general rather than the moon specifically. All celestial bodies visible from Earth, including our moon, are fixed on the celestial sphere, and this is what rotates relative to our position on Earth. So instead of focusing on a lunar-centric solution, I need to figure out how to calculate the rotation of the celestial sphere itself. Still haven't wrapped my head around this one.

Update 4/13/14

I found an article on Crescent Moon Tilt which may be the key to solving my problem. The math is obviously out there/can be derived, since most astronomical simulation software can show the exact location and perspective from earth for any celestial body. I think the hard part is that the single piece I want is wrapped into much more complicated equations that cover overall positioning. I will work to derive a simpler version by reading and re-reading this helpful information.

Answer 1

You need 4 things to get Moon light direction as seen from Earth:

1- Sun altitude. 2- Sun azimuth. 3- Moon altitude. 4- Moon azimuth. [All altitudes are airless (without refraction effect), and all angles in degrees.]

I use this PHP function to get the light's rotation from these things:

function getAngle($sunalt, $sunaz,$moonalt, $moonaz) {
$dLon = ($sunaz - $moonaz);
$y = sin(deg2rad($dLon)) * cos(deg2rad($sunalt));
$x = cos(deg2rad($moonalt)) * sin(deg2rad($sunalt)) - sin(deg2rad($moonalt)) * cos(deg2rad($sunalt)) * cos(deg2rad($dLon));
$brng = atan2($y, $x);
$brng = rad2deg($brng);  
return $brng;
}

This is the direction of the crescent in degrees clockwise.

Answer 2

To make future research easier, I believe what you're looking for is known as an ephemeris, which gives "the positions of naturally occurring astronomical objects as well as artificial satellites in the sky at a given time or times".

In Python, there are at least two good choices pyephem and ephem, both of which use the XEphem c library, but alternative exist for most languages including perl.

There also appear to be some related questions about moon / lunar phases already on stackoverflow

Answer 3

I've been experimenting with the JPL DE430/431 ephemerides and NOVAS-C. My work is available on GitHub:

https://github.com/brhamon/astro

Included in the source code is a small application (planets/planets.c) that calculates the positions of the planets as viewed by an observer on Earth.

I've added a call to moon_phase, implemented in ephutil/moon.c. It finds the moon phase in equatorial spherical coordinates that provide the phase latitude and longitude. The phase latitude is the number of degrees north of the Earth transit point on the Moon. The phase longitude is the number of degrees east of the Earth transit point on the Moon. Dividing phase longitude by 45 degrees produces the common phase name, like "last quarter" or "waxing crescent".

The Solar termination on the Moon is a hemisphere centered on the Solar transit point on the Moon. You should be able to convert these phase coordinates into a line representing the Solar termination on a model of the Moon.

Finding the phase coordinates required first calling NOVAS to get the position vectors of the Sun and Moon with respect to an observer at the center of Earth's mass. The difference of these two vectors yields the position vector of the Sun with respect to an observer at the center of Moon's mass. Multiplying the Earth-Moon vector by -1 reverses it, so it becomes a position vector from Moon to Earth.

Next, I convert these two vectors into Right Ascension and Declination; which are spherical coordinates in the "equinox of date" (although the exact origin point becomes irrelevant when we subtract them). The result is the phase coordinates, which provide the angular distance north and east of the Earth's transit point. From Earth, this point will be approximately in the center of the Moon's disk.

/*
 * Calculate the equatorial spherical coordinates of the solar transit point with
 * respect to the center of the Earth-facing surface of the Moon.
 */
short int moon_phase(time_parameters_t* tp,
        object* sun, object* moon, short int accuracy,
        double* phlat, double* phlon, int* phindex)
{
    short int error = 0;
    observer geo_ctr;
    sky_pos t_place;
    double lon1, lat1;
    double lon2, lat2;
    double earth_sun[3];
    double moon_earth[3];
    double moon_sun[3];

    make_observer(0, NULL, NULL, &geo_ctr);
    if ((error = place(tp->jd_tt, sun, &geo_ctr, tp->delta_t,
                    coord_equ, accuracy, &t_place)) != 0) {
        printf("Error %d from place.", error);
        return error;
    }
    radec2vector(t_place.ra, t_place.dec, t_place.dis, earth_sun);
    if ((error = place(tp->jd_tt, moon, &geo_ctr, tp->delta_t,
                    coord_equ, accuracy, &t_place)) != 0) {
        printf("Error %d from place.", error);
        return error;
    }
    radec2vector(t_place.ra, t_place.dec, t_place.dis, moon_earth);
    /* The vector points from Earth to Moon. Reverse it. */
    moon_earth[0] *= -1.0;
    moon_earth[1] *= -1.0;
    moon_earth[2] *= -1.0;

    /* Calculate the position vector of the Sun w/r/t Moon */
    moon_sun[0] = earth_sun[0] + moon_earth[0];
    moon_sun[1] = earth_sun[1] + moon_earth[1];
    moon_sun[2] = earth_sun[2] + moon_earth[2];

    vector2radec(moon_earth, &lon1, &lat1);
    lon1 *= 15.0;
    if (lon1 > 180.0) {
        lon1 -= 360.0;
    }
    /* {lat1, lon1} is now the equatorial spherical coordinates of the
     * Earth transit point on the Moon. */

    vector2radec(moon_sun, &lon2, &lat2);
    lon2 *= 15.0;
    if (lon2 > 180.0) {
        lon2 -= 360.0;
    }
    /* {lat2, lon2} is now the equatorial spherical coordinates of the
     * Solar transit point on the Moon. */

    *phlon = normalize(lon2 - lon1, 360.0);
    *phindex = (int)floor(normalize(lon2 - lon1 + 22.5, 360.0)/45.0);
    if (*phlon > 180.0) {
        *phlon -= 360.0;
    }
    *phlat = lat2 - lat1;
    return error;
}