Moon parallactic, position, and zenith angles

[sbooth]

(Related to #253, #278)

Thank you for making this library available!

I wanted to calculate the position and zenith angles for the Moon's illuminated limb and wasn't aware of the camera demo which does essentially that. I wrote an experimental demo to perform the calculations based on Meeus' Astronomical Algorithms chapter 48.

/**
 * @brief Calculates the parallactic angle `q` for a body.
 *
 * The parallactic angle `q` is the angle between a body's zenith (the uppermost point of the
 * body in the sky as seen by an observer) and the northern celestial pole.
 * `q` is the angle between a body's vertical and its hour circle.
 *
 * @param body
 *     The celestial body to be observed. Not allowed to be `BODY_EARTH`.
 * @param time
 *     The date and time at which the observation takes place.
 * @param observer
 *     A location on or near the surface of the Earth.
 * @param aberration
 *     Selects whether or not to correct for aberration.
 *
 * @return
 *     The body's parallactic angle.
 */
astro_angle_result_t Astronomy_ParallacticAngle(astro_body_t body, astro_time_t *time, astro_observer_t observer, astro_aberration_t aberration)
{
    astro_equatorial_t body_equator_of_date;
    astro_func_result_t body_hour_angle;
    double H, phi, delta;
    double q;
    astro_angle_result_t result;

    if (time == NULL)
        return AngleError(ASTRO_INVALID_PARAMETER);

    body_equator_of_date = Astronomy_Equator(body, time, observer, EQUATOR_OF_DATE, aberration);
    if (body_equator_of_date.status != ASTRO_SUCCESS)
        return AngleError(body_equator_of_date.status);

    /*
     Calculate the parallactic angle (q) using
     Meeus equation 14.1 (p. 98).
     */
    body_hour_angle = Astronomy_HourAngle(body, time, observer);
    if (body_hour_angle.status != ASTRO_SUCCESS)
        return AngleError(body_hour_angle.status);

    H = body_hour_angle.value * HOUR2RAD;
    phi = observer.latitude * DEG2RAD;
    delta = body_equator_of_date.dec * DEG2RAD;

    q = atan2(sin(H), tan(phi) * cos(delta) - sin(delta) * cos(H));

    result.angle = q * RAD2DEG;
    result.status = ASTRO_SUCCESS;

    return result;
}

/**
 * @brief Calculates the position angle `chi` of a body's bright limb.
 *
 * The angle `chi` is the position angle of the midpoint of the illuminated limb
 * of the body reckoned eastward from the north point of the disk (not from the axis
 * of rotation of the globe). The position angles of the cusps are `chi` ± 90°.
 *
 * @param body
 *     The celestial body to be observed. Not allowed to be `BODY_EARTH`.
 * @param time
 *     The date and time at which the observation takes place.
 * @param observer
 *     A location on or near the surface of the Earth.
 * @param aberration
 *     Selects whether or not to correct for aberration.
 *
 * @return
 *     The body's bright limb position angle.
 */
astro_angle_result_t Astronomy_BrightLimbAngle(astro_body_t body, astro_time_t *time, astro_observer_t observer, astro_aberration_t aberration)
{
    astro_equatorial_t body_equator_of_date;
    double delta, alpha;
    astro_equatorial_t sun_equator_of_date;
    double delta0, alpha0;
    double chi;
    double body_position_angle;
    astro_angle_result_t result;

    if (time == NULL)
        return AngleError(ASTRO_INVALID_PARAMETER);

    /*
     Calculate the position angle of the body's bright limb (chi) using
     Meeus equation 48.5 (p. 346).
     */
    body_equator_of_date = Astronomy_Equator(body, time, observer, EQUATOR_OF_DATE, aberration);
    if (body_equator_of_date.status != ASTRO_SUCCESS)
        return AngleError(body_equator_of_date.status);

    delta = body_equator_of_date.dec * DEG2RAD;
    alpha = body_equator_of_date.ra * HOUR2RAD;

    sun_equator_of_date = Astronomy_Equator(BODY_SUN, time, observer, EQUATOR_OF_DATE, aberration);
    if (sun_equator_of_date.status != ASTRO_SUCCESS)
        return AngleError(sun_equator_of_date.status);

    delta0 = sun_equator_of_date.dec * DEG2RAD;
    alpha0 = sun_equator_of_date.ra * HOUR2RAD;

    chi = atan2(cos(delta0) * sin(alpha0 - alpha), sin(delta0) * cos(delta) - cos(delta0) * sin(delta) * cos(alpha0 - alpha));

    body_position_angle = chi * RAD2DEG;
    while (body_position_angle < 0)
        body_position_angle += 360;

    result.angle = body_position_angle;
    result.status = ASTRO_SUCCESS;

    return result;
}

Their use in my demo:

moon_parallactic_angle = Astronomy_ParallacticAngle(BODY_MOON, &time, observer, ABERRATION);
moon_position_angle = Astronomy_BrightLimbAngle(BODY_MOON, &time, observer, ABERRATION);

/*
    The angle χ is not measured from the direction of the observer's zenith.
    The zenith angle of the bright limb is χ - q.
*/
moon_zenith_angle = moon_position_angle.angle - moon_parallactic_angle.angle;

I'm new to astronomy so I'm not sure if the above is completely correct. However, the resulting zenith angle is similar to the result from camera although the angles have different signs. As an example, my luna demo calculated the zenith angle as 310.466 degrees while camera calculated -49.533 degrees.

Any corrections or suggestions for improvement are appreciated!

[cosinekitty]

These functions look great. About correctness, I guess the best thing would be to find some source of test data from a reliable source, and see if we can calculate the same numbers using these functions. Did the Meeus book give any example problems with answers?

[sbooth]

The paper Position of the bright limb of the Moon derives the same formulas for the parallactic and position angles and gives several worked examples for 2030 on pages 9-10. I haven't had the chance to run comparisons for each of their data points but they generally seem to agree within 2 ° of the calculations using the functions above.

[sbooth]

I compared the values in that paper labeled Parallactic angle geocentric, Geocentric angle with respect to the hour circle, and Topocentric angle with respect to the vertical with the values from my luna demo called Parallactic angle, Position angle, and Zenith angle and they generally agreed within 1°.

[sbooth]

The book doesn't give any examples for parallactic angle calculations. It does say the following about the parallactic angle for the Sun or Moon:

C is the center of the disk. The direction of the zenith and that of the celestial North are indicated. Z is the zenith point of the disk; it is the uppermost point of the disk at the sky as seen by the observer at the given instant. N is the north point of the disk; the direction CN points towards the northern celestial pole.

The angle ZCN is called the parallactic angle and is generally designed by q.

This parallactic angle has absolutely nothing to do with the parallax! The name arises from the fact that the celestial body moves along a parallel circle. Compare with the "parallactic" mounting of a telescope.

By convention, the angle q is negative before, and positive after the passage through the southern meridian. Exactly on the meridian, we have q = 0°.

There is one example for calculating χ. It says for the Moon on 1992 April 12 at 0h Dynamical Time, χ = 285.0 °.

The TT-UT1 delta for April 1992 was 58.5401 seconds so I believe that 0h TT was approximately 23:59:01 UTC the day prior.

The luna demo calculated χ = 284.99 degrees for 1992-04-11T23:59:01.4599Z so for the example it's close.

[cosinekitty]

Your functions ParallacticAngle and PositionAngle are nicely done, and I wonder if you would consider donating them to Astronomy Engine? I still have a rule that all functions need some kind of unit test to exercise them before I publish them in the master branch. If you want to take a whack at any kind of initial test, that would be great. Or I could take a look at it next week. This week I'm super focused on another project.

The TT-UT1 delta for April 1992 was 58.5401 seconds so I believe that 0h TT was approximately 23:59:01 UTC the day prior.

It might be helpful to know that you can call the function Astronomy_TerrestrialTime and pass it what Meeus calls dynamical time as what Astronomy Engine calls tt or Terrrestrial Time. That function calculates the ΔT value to deduce the correct civil UTC in ut.

[sbooth]

I'm definitely happy to contribute something to the project! I haven't delved much into the build process but I'm willing to give it a shot and submit a PR. I'd like to clean up the code slightly first, though, and add a few more comments. Is UTF-8 allowed in source files or should variables like δ0 be renamed to delta_nought?

[cosinekitty]

I'm trying to keep the C code as portable as possible. It would be a good idea to use ASCII symbols only in the variable names.

[sbooth]

I've edited the code in the original post to add comments and remove any non-ASCII characters. I'm still working on getting the build tools installed so I can submit a PR.

One thing I wonder is whether the name Astronomy_PositionAngle is a bit generic. Would Astronomy_BrightLimbPositionAngle be more descriptive?

[cosinekitty]

I agree the PositionAngle name is generic now that you mention it. How about Astronomy_BrightLimbAngle?

[sbooth]

Done!

Also, you're welcome to include the luna demo or a portion thereof if you're interested. It doesn't add much above the current demos other than the bright limb and zenith angle calculations (and some sweet Unicode box-drawing characters).