Lunatics

sometimes ask for moon computations for the AVR. I generally shy away from that, because of the constraints imposed. To do it well requires a lot of high precision floating point math, involving a bazillion transendental function calls. With the limitations of the avr 'double', the results are un-impressive.

Some of those operations can be done away with by the use of interpolation tables. That works well, for a limited range of time. But tables grow quickly as the range increases.

Using either long form formulas or huge lookup tables, in the general avr context, is... lunacy.

However, if we can be satisfied with an approximation, and allow an error of a few degrees, we can write functions that do not consume lots of cpu resources.
Elongation and Moon Rise
At solar noon, on the day of the 'new moon', Earth, Moon, and Sun are lined up... the angle between moon and sun is zero. If we check again at noon the following day, we will notice that the angle between the Sun and Moon has increased , and a tiny portion of the Moons Eastern hemisphere is visible.

Continuing as the moon 'waxes', the angle has increased to 90 degrees at the time of the First Quarter, when the Moons entire Eastern Hemisphere is lit.

This angle is called the Elongation of the Moon. It's easy to see that the Phase and Elongation are related. AVR-Libc provides a function, moon_phase(), that returns the phase as a percentage. We can get the Elongation easily from moon_phase()...

Given the elongation of the moon El in hours, we can say that the Moon reaches its highest point in the sky El hours before the Sun does. Indeed, for any lunar 'event', that event will occur El hours before the equivalent solar event... and that includes the rising and setting of the moon. Since AVR-Libc provides the functions sun_rise() and sun_set(), we can proceed to determine the time of the equivalent lunar events...
time_t MoonRise(time_t * point_in_time){

  time_t sunrise = sun_rise(point_in_time);
  long elongation =  moon_phase(&sunrise) * 432L;
  return  sunrise + elongation;

}

Wasn't that too easy?! Lets do MoonSet()...
time_t MoonSet(time_t * point_in_time){

  time_t sunset = sun_set(point_in_time);
  long elongation =  moon_phase(&sunset) * 432L;
  return  sunset + elongation;

}
Hour Angle and Moon Position
A typical lunatic is satisfied with the moons phase and time of rising / setting. The most outrageously lunatic, however, want to keep track of where the moon is at all times. To obtain the altitude and azimuth of the moon, we must know the moons Hour Angle and its Declination.

The declination is easy. The Moon orbits the Earth almost on the same plane that the Earth orbits the Sun, the 'Ecliptic' plane, tilted by about 5 degrees. We can approximate the declination of the moon simply by substituting the declination of the sun. Since AVR-Libc provides the solar_declination() function, we have that part of it in the bag.

Now, about the Hour Angle of the Moon, HAmoon. We said previously that, for every Solar event, the equivalent Lunar event lags or leads by the amount of time of the Elongation, El. So, if we know HAsun sun, we just subtract El to get HAmoon. We learned how to get the hour angle of the sun in a previous article, so lets write some code to derive the HourAngle of the Moon...

#include <math.h>
#include <time.h>
#define RAD_TO_DEGREE 57.29577951308233
double MoonHourAngle(time_t * point_in_time){

  // compute hour angle of sun (in seconds)
  time_t noon = solar_noon(point_in_time) % ONE_DAY;      // time of noon mod one day
  time_t tday = *point_in_time % ONE_DAY;                 // current time mod one day
  long HAsun = tday - noon;                               // offset from noon, -43200 to + 43200 seconds

  // compute elongation of moon in seconds
  long elongation =  moon_phase(point_in_time) * 432L;

  // compute moon hour angle in seconds
  double HAmoon = HAsun - elongation;

  // convert to radians
  HAmoon *= M_PI;
  return HAmoon /  43200.0;
}


Given the Moons Hour Angle and Declination, we are now in a position to obtain the altitude and azimuth. We will re-use the AltAz() function described in the previous article, and reproduced here...
void AltAz(double HourAngle, double Declination, double * Azimuth, double * Altitude){
  extern long  __latitude;
  double Latitude =  __latitude / ( ONE_DEGREE * RAD_TO_DEGREE);          // latitude in radians

  double CosHourAngle = cos(HourAngle);
  double SinLatitude = sin(Latitude);
  double CosLatitude = cos(Latitude);

  *Azimuth = atan ( sin(HourAngle) / ( CosHourAngle * SinLatitude - tan(Declination ) * CosLatitude ) );
  *Altitude = asin( CosHourAngle * cos(Declination) * CosLatitude + sin(Declination) * SinLatitude );
}

So finally, to obtain the moons position...

void MoonPosition(time_t * point_in_time, double * Azimuth, double * Altitude){

  double HAmoon = MoonHourAngle(point_in_time);

  double Declination = solar_declination(point_in_time);

  AltAz( HAmoon, Declination,  Azimuth,  Altitude );

}
Questions and comments can be sent to