Calculating the Moon's age & phase
Author Message
Calculating the Moon's age & phase

G'day all,

A small donation to the public domain...

The program (below) computes the % illumination and "age" of the Moon;
to an accuracy sufficient for fishermen, artists, romantics, etc....
It's also good enough for general astronomical scheduling - although I
wouldn't predict eclipses with it <g>.

SWAG note : this is a little more accurate than SWAG's existing Moon
phase algorithms.

cheers,
Fraser Farrell

=== cut here =====================================
program moon;

{
Compute illumination and "age" of Moon.  Illumination is the
% of visible lunar surface in sunlight (0 = new moon, 100 = full).
Age is the number of days since new moon (range 0 to ~29.53 days).

Uses low-precision formulae from Chapters 45-46 of Jean Meeus'
"Astronomical Algorithms" (Willmann-Bell, 1991).
Errors, in comparison to the definitive ELP-2000/82 Lunar
Theory, during the years 1900-2100:
Illumination - up to 2%
Age - up to 0.3 day

You will need Extended math to retain enough digits in calculations.

Written by Fraser Farrell.   USE OR MODIFY THIS CODE AS YOU WISH... :)

Quote:
}

{\$N+,E+}

uses  dos,crt;  {only needed for ClrScr, GetDate, GetTime}

var   T, phase_angle : extended;
moon_age, time_zone : real;
moon_light : integer;
year,month,day,dow,hour,minute,second,sec100 : word;

{-------------------------------------------------------------------}
{     first of all, a few functions useful in astronomy...          }
{-------------------------------------------------------------------}

function unwindangle (a : extended) : extended;
{converts very large angles to the range 0-360 degrees}
begin
while a >= 360 do  a:=a-360;
while a < 0 do  a:=a+360;
unwindangle:=a;
end;

function degtorad(d : extended) : extended;
begin
end;

function radtodeg(d : extended) : extended;
begin
end;

{-------------------------------------------------------------------}
{          the Local Time -> Julian Day converter...                }
{-------------------------------------------------------------------}

function julianday(yr,mo,d,h,min : word; zone:real): extended;
{convert the local date & time to Julian Day}

var     a,b,jd : extended;

begin
if mo < 3 then a:=int((yr-1)/100)
else a:=int(yr/100);
b:=2-a+int(a/4);
if mo < 3 then jd:=int(365.25*(yr-1)) + int(30.6001*(mo+13))
else jd:=int(365.25*yr) + int(30.6001*(mo+1));
julianday:=jd + d + (h+min/60)/24 + 1720994.5 + b - zone/24;
end; {julianday}

{-------------------------------------------------------------------}
{          The Moon's MEAN longitude, anomaly & elongation.         }
{   Ignores the Moon's latitude, and 1000's of higher-order terms!  }
{-------------------------------------------------------------------}

function moon_mean_longitude(T:extended):extended;
begin
moon_mean_longitude:=
+T*(481267.88134236 +T*(-1.3268e-3
+T*(1/538841-T/65194000)))));
end;

function moon_mean_anomaly(T:extended):extended;
begin
moon_mean_anomaly:=
+T*(477198.8676313 +T*(8.997e-3
+T*(1/69699-T/14712000)))));
end;

function moon_mean_elongation(T:extended):extended;
begin
moon_mean_elongation:=
+T*(445267.1115168 +T*(-1.63e-3
+T*(1/545868-T/113065000)))));
end;

{-------------------------------------------------------------------}
{     The Sun's MEAN anomaly, ignores many higher-order terms!      }
{-------------------------------------------------------------------}

function sunanomaly(T:extended):extended;
var M,C : extended;
begin
M:= 357.5291 + T*(35999.0503 - T*(1.559e-4 + T*4.8e-7));
C:= sin(M)*(1.9146-T*(0.004817+T*1.4e-5))
+ sin(2*M)*(0.019993-T*1.01e-4)
+ sin(3*M)*2.9e-4;
+T*(35999.0503-T*(1.559e-4 +T/24490000))));
end;

{-------------------------------------------------------------------}
{  The main program                                                 }
{-------------------------------------------------------------------}
BEGIN
clrscr;
writeln('MOON shows approximate illumination & age of moon');
writeln('using computer''s clock, and the Time Zone you input...');
writeln;
write('What is your Time Zone in hours (negative if behind UT) ? ');
writeln; writeln;
getdate(year,month,day,dow);
gettime(hour,minute,second,sec100);
writeln('The date is : ',year:5,month:3,day:3);
writeln('The time is : ',hour:3,minute:3,second:3);
writeln; writeln;
writeln('The Julian Day is ',
julianday(year,month,day,hour,minute,time_zone):14:4);
writeln; writeln;

T:= (julianday(year,month,day,hour,minute,time_zone)
- 2451545) / 36525;

phase_angle:= 180
- 6.289*sin(moon_mean_anomaly(T))
+ 2.100*sin(sunanomaly(T))
- 1.274*sin(2*moon_mean_elongation(T)
-moon_mean_anomaly(T))
- 0.658*sin(2*moon_mean_elongation(T))
- 0.214*sin(2*moon_mean_anomaly(T))
-0.110*sin(moon_mean_elongation(T));

moon_light:= trunc(100 * (1+cos(phase_angle))/2);
moon_age:= 29.530588*moon_mean_elongation(T)/(2*pi);

writeln('The Moon is approximately ',moon_light:3,'% illuminated');
writeln('The Moon''s age is approximately ',moon_age:3:1,' days');

writeln; writeln; writeln('press ENTER to continue...');
END.

=== cut here ===================================================

Wed, 18 Jun 1902 08:00:00 GMT
Calculating the Moon's age & phase

Quote:

>A small donation to the public domain...

>The program (below) computes the % illumination and "age" of the Moon;
>to an accuracy sufficient for fishermen, artists, romantics, etc....
>It's also good enough for general astronomical scheduling - although I
>wouldn't predict eclipses with it <g>.

I can't test it now; the moon should have just risen, but it is a wet
day here.  Might it be easy to make the program also calculate the
moonrise/moonset (or the time when the moon is over the central meridian
of the time-zone, which should need no extra data)?  Duffett-Smith has
published algorithms.

It is quite clear that the program is intended for use with the CURRENT
DOS date, which must now lie within the range 1996/12/20 .. about 2100.
Nevertheless, I would like to see an indication of the range of dates
over which the internal algorithms are satisfactory.

The astronomical constants in the calculation cannot be infinitely
precise, so accuracy must degrade to some extent away from the present
day.

The Julian Date routine evidently requires the Gregorian calendar at
least; and surely understands quadrennial leap-years; but it is not
clear to me as yet whether it implements the 100 & 400 year rules
(which, in 2000, cancel each other).  It agrees with the (quite
different) one I have over 99.997 years each side of 2000/02/29, and
that is all that mine is good for.

Actually, I'd quite like a really good JD routine (no float arithmetic;
wide range of validity; efficient) to link or put on my WWW site
(below).
--

http://www.merlyn.demon.co.uk/

Wed, 18 Jun 1902 08:00:00 GMT
Calculating the Moon's age & phase

>>The program (below) computes the % illumination and "age" of the
>>Moon; to an accuracy sufficient for fishermen, artists, romantics,

Maybe I should have mentioned I wanted modest code size too?  This
program is a small part of a (much larger) variable star data analysis
package I'm upgrading; and approx Moon info is adequate for my purposes.

>wet day here.  Might it be easy to make the program also calculate
>the moonrise/moonset (or the time when the moon is over the central
>meridian of the time-zone, which should need no extra data)?
>Duffett-Smith has published algorithms.

I don't have a copy of Duffet-Smith's book on hand and my local library
have managed to lose their copy.  Jean Meeus also presents moonrise/set
algorithms but a fair bit of math needs to be done.  At the moment I am
experimenting with abbreviations of his routines to get a tradeoff
between code size & accuracy.

As for the "central meridian of the time-zone", UT+9h30m where I live
has the (AFAIK) unique distinction of a "central meridian" entirely
_outside_ its territory!  Even when Summer Time operates.  It's a very
chronological lunacy....

>It is quite clear that the program is intended for use with the
>CURRENT DOS date, which must now lie within the range 1996/12/20 ..
>about 2100. Nevertheless, I would like to see an indication of the
>range of dates over which the internal algorithms are satisfactory.

The GetDate/GetTime was included simply to demonstrate the program.
There is no reason why the user cannot be prompted to enter the
year,month,day,etc.  Sorry, I should have explained this better.

I have tested my algorithm against predictions (ELP-2000/82 Lunar
Theory) for many dates/times in the range 1900-2100 and got good (for my
purposes) results. I also tested a lot of pre-1900 eclipses to check
that my algorithm predicted a new (or full) moon; and found that the
accuracy decreased noticeably for years <1000.

>The astronomical constants in the calculation cannot be infinitely
>precise, so accuracy must degrade to some extent away from the
>present day.

Indeed it does.  You will have noticed that the calculations use
polynomial functions of time, which is standard practice for lunar &
planetary computations.  These _approximations_ are only valid for a
limited number of centuries either side of 2000; so using them eg: to
compute phenomena for the year 50000 BC is invalid and will produce
silly results.  This is one trap that pseudo-scientists fall into when
they "prove" the planets cannot have existed in their present orbits
until recently....

Even the best approximations of planetary motions have limits.  Last
time I asked, the stability of the solar system had been mathematically
proven only for the last 840 million years - and this was on a custom
built supercomputer running algorithms that made the current VSOP87

>The Julian Date routine evidently requires the Gregorian calendar at
>least; and surely understands quadrennial leap-years; but it is not
>clear to me as yet whether it implements the 100 & 400 year rules

The Julian Day routine works for all Gregorian dates (1582 Oct 15 and
later) and does use the century rules.  To get it working for Julian
calendar dates instead (pre- 1582 Oct 4), this small modification is
needed:

replace     b:= 2-a+int(a/4);     with     b:=0;

I didn't bother to code this as I was only interested in Gregorian
calendar dates.  To cope with both calendars you would need to examine
the inputted date and compute "b" appropriately.  And for the truly
rigorous, check the country as well since the Gregorian calendar was
adopted at different times in each one!

A couple more traps to watch out for:  astronomers recognise a Year 0,
historians do not; and Julian Days begin at 12:00 UT....

>Actually, I'd quite like a really good JD routine (no float
>arithmetic; wide range of validity; efficient) to link or put on my
>WWW site (below).

Meeus offers an equation which calculates the JD for January 0.0 of any
year in the range 1901-2099 inclusive:

JD_Jan_0 := 1721409.5 + int(365.25*(year-1));

Add the day number to this to get the JD of a date.  The day number
[1..365 or 366] of any day in any year (according to Meeus) is:

day_number := int(275*month/9) - K*int((month+9)/12) + day - 30;
where K = 1 for leap years, K = 2 for other years

There is also the Modified Julian Day which is defined:

MJD:= JD - 2400000.5;

These begin at 00:00 UT (note the 0.5) and should be calculable using
Reals/Integers. I hope to see your definitive Julian Day function posted
in here before JD 2450500...!  ;)

cheers,
Fraser Farrell

As easy as 3.141592653589793238462643383279502884197169399375105...

Wed, 18 Jun 1902 08:00:00 GMT

 Page 1 of 1 [ 3 post ]

Relevant Pages