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;

{degrees to radians conversion}

begin

degtorad:=d*pi/180;

end;

function radtodeg(d : extended) : extended;

{radians to decimal degrees conversion}

begin

radtodeg:=d*180/pi;

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:=

degtorad(unwindangle(218.3164591

+T*(481267.88134236 +T*(-1.3268e-3

+T*(1/538841-T/65194000)))));

end;

function moon_mean_anomaly(T:extended):extended;

begin

moon_mean_anomaly:=

degtorad(unwindangle(134.9634114

+T*(477198.8676313 +T*(8.997e-3

+T*(1/69699-T/14712000)))));

end;

function moon_mean_elongation(T:extended):extended;

begin

moon_mean_elongation:=

degtorad(unwindangle(297.8502042

+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));

M:= degtorad(M);

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;

sunanomaly:= degtorad(unwindangle(C+357.5291

+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) ? ');

readln(time_zone);

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

- radtodeg(moon_mean_elongation(T))

- 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));

phase_angle:= degtorad(phase_angle);

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...');

readln;

END.