fred@inuxe.UUCP (Fred Mendenhall) (08/09/85)
*** REPLACE THIS LINE WITH YOUR MESSAGE ***
program moon;
var {global variables}
time,mm,dd,yy,jdn,jd,t,theta,lambda,beta,mmanom,msanom : real;
mra,mdec : real;
continue : boolean;
answer : char;
procedure start;
var
i : integer;
stz : real;
begin
clrscr;
writeln;
writeln;
writeln(' PROGRAM MOONBEAM ');
writeln(' version 1.01');
writeln(' Copyright 1985 by');
writeln(' F. T. Mendenhall Jr.');
for i :=1 to 6 do
writeln;
writeln('Please supply the date of interest:');
write('What is the month (1-12)? ');
readln(mm);
write('What is the day (1-31)?');
readln(dd);
write('What is the year?');
readln(yy);
clrscr;
writeln(' PLEASE ENTER APPROPRIATE TIME ZONE');
writeln;
writeln(' Enter 0 for Universal Time (GMT)');
writeln;
writeln(' Enter 1 for Eastern Standard Time');
writeln;
writeln(' Enter 2 for Central Standard Time');
writeln;
writeln(' Enter 3 for Mountain Standard Time');
writeln;
writeln(' Enter 4 for Pacific Standard Time');
writeln;
write(' What time zone do you wish to use (0-4)? ');
readln(stz);
clrscr;
writeln;
writeln(' USING THE 24 HOUR CLOCK WHAT IS THE TIME OF OBSERVATION?');
writeln;
writeln('NOTE: 1pm=1300 2pm=1400 3pm=1500');
writeln(' 4pm=1600 5pm=1700 6pm=1800');
writeln(' 7pm=1900 8pm=2000 9pm=2100');
writeln(' 10pm=2200 11pm=2300 12am=0000');
writeln;
writeln(' for example 8:45pm would entered as 2045 hours');
writeln;
write(' What is the time of your observation? ');
readln(time);
clrscr;
writeln;
write('For the date ',mm:2:0,'/',dd:2:0,'/',yy:4:0);
writeln(' At ',time:4:0,' Local Time' );
if stz <> 0 then
time := time + stz*100 +400;
end;{end of procedure start}
procedure julian;
var
m,y,a,b : real;
begin { begin julian day procedure}
dd := dd + Int(time/100.0)/24.0 + Frac(time/100.0)/24.0;
if mm > 2 then
begin
y := yy;
m := mm;
end;{end of then clause}
if (mm=1) or (mm=2) then
begin
y := yy-1;
m := mm +12;
end; {end of then clause}
if yy < 1583 then
writeln('WARNING PROGRAM NOT VALID PRIOR TO 1583');
a := int(y/100);
b := 2-a+int(a/4);
jd := int(365.25*y)+int(30.6001*(m+1))+dd+1720994.5+b;
end; {end of procedure julian}
procedure julian_new;
{ This procedure calculates the previous date of the new moon
from the date of interest}
var
jds,ms,ds,k,y,tt,m,m1,f,cor : real;
begin { begin procedure julian_new}
ms := mm;
{ store the original date and julian day}
ds := dd;
mm := 1;
dd := 1;
jds := jd;
{ find the julian day of the first of the year}
julian;
tt := jds -jd;
{ tt is the number of days since the first of the year}
y := yy+tt/365.;
jd := jds;
{restore the original date and jd}
mm := ms;
dd := ds;
k := (y-1900)*12.3685;
k := int(k);
tt := k/1236.85;
jdn := 2415020.75933+29.53058868*k+0.0001178*tt*tt
-0.000000155*tt*tt*tt+0.00033*sin((166.56+132.87*tt
-0.009173*tt*tt)*pi/180.0);
m := (359.2242+29.10535608*k-0.0000333*tt*tt
-0.00000347*tt*tt*tt)*pi/180.0;
m1 := (306.0256+385.81691806*k+0.0107306*tt*tt
+0.00001236*tt*tt*tt)*pi/180.0;
f := (21.2964+390.67050646*k-0.0016528*tt*tt
-0.00000239*tt*tt*tt)*pi/180.0;
cor := (0.1734-0.000393*tt)*sin(m)
+0.0021*sin(2*m)
-0.4068*sin(m1)
+0.0161*sin(2*m1)
-0.0004*sin(3*m1)
+0.0104*sin(2*f)
-0.0051*sin(m+m1)
-0.0074*sin(m-m1)
+0.0004*sin(2*f+m)
-0.0004*sin(2*f-m)
-0.0006*sin(2*f+m1)
+0.0010*sin(2*f-m1)
+0.0005*sin(m+2*m1);
jdn := jdn +cor;
end;
{ jdn is the julian of the previous new moon }
procedure phaseout;
{ This procedure prints out the moon phase }
var
dpn : real;
begin
writeln;
dpn := jd-jdn;
writeln(' The Moon will be ',dpn:3:1,' days passed new');
writeln;
if dpn < 7.5 then
writeln(' The Moon is between New and First Quarter');
if (dpn > 7.5) and (dpn <14) then
writeln (' The Moon is between First Quarter and Full');
If (dpn >14) and (dpn < 21.5) then
writeln(' The Moon is between Full and Last Quarter');
if dpn > 21.5 then
writeln(' The Moon is between Last Quarter and New');
writeln;
writeln('STANDBY: Lunar position being calculated');
end;{end of phaseout}
procedure sunlat;
{This procedure finds the mean sun anomaly and the suns true
longitude theta}
var
l,c : real;
begin {sunlat}
t := (jd-2415020.0)/36525.0;
l := 279.69668 + 36000.76892*t + 0.0003025*t*t;
msanom := 358.475833 + 35999.04975*t - 0.000150*t*t - 0.0000033*t*t*t;
c := +(1.919460-0.004789*t-0.000014*t*t)*sin(msanom*pi/180.0)
+(0.020094 - 0.000100*t)*sin(2.0*msanom*pi/180.0)
+ 0.000293*sin(3.0*msanom*pi/180.0);
theta := l + c;
end; {end procedure sunlat}
procedure moonpos;
{This procedure calculates the moons mean anomaly, and its
longitude, lambda, and latitude, beta}
var
l,d,f,omeg,epsio,addcon,addcon2,rad : real;
begin {moonpos, note that t,msanom are previously calculated in sunlat}
rad := pi/180.0;
epsio := 1.0 - 0.002495*t - 0.00000752*t*t;
addcon := 0.003964*sin((346.560+132.870*t-0.0091731*t*t)*rad);
addcon2 := sin((51.2+20.2*t)*rad);
l := 270.434164 + 481267.8831*t -0.001133*t*t + 0.0000019*t*t*t;
mmanom := 296.104608 + 477198.8491*t +0.009192*t*t +0.0000144*t*t*t;
d := 350.737486 +445267.1142*t -0.001436*t*t +0.0000019*t*t*t;
f := 11.250889 +483202.0251*t -0.003211*t*t -0.0000003*t*t*t;
omeg := 259.183275 -1934.1420*t +0.002078*t*t +0.0000022*t*t*t;
{ now adding periodic variations}
l := l +addcon +0.001964*sin(omeg*rad) +0.000233*addcon2;
msanom := msanom -0.001778*addcon2;
mmanom := mmanom +addcon +0.000817*addcon2 +0.002541*sin(omeg*rad);
d := d +addcon +0.002011*addcon2 +0.001964*sin(omeg*rad);
f := f+addcon -0.024691*sin(omeg*rad)
-0.004328*sin((omeg+275.05-2.3*t)*rad);
{calculate the coordinates}
lambda := l + 6.288750*sin(mmanom*rad)
+1.274018*sin((2*d-mmanom)*rad)
+0.658309*sin(2*d*rad)
+0.213616*sin(2*mmanom*rad)
-0.185596*epsio*sin(msanom*rad)
-0.114336*sin(2*f*rad)
+0.058793*sin((2*d-2*mmanom)*rad)
+0.057212*epsio*sin((2*d-msanom-mmanom)*rad)
+0.053320*sin((2*d+mmanom)*rad)
+0.045874*epsio*sin((2*d-msanom)*rad)
+0.041024*epsio*sin((mmanom-msanom)*rad)
-0.034718*sin(d*rad)
-0.030465*epsio*sin((mmanom+msanom)*rad)
+0.015326*sin((2*d-2*f)*rad)
-0.012528*sin((2*f+mmanom)*rad)
-0.010980*sin((2*f-mmanom)*rad)
+0.010674*sin((4*d-mmanom)*rad)
+0.010034*sin(3*mmanom*rad)
+0.008548*sin((4*d-2*mmanom)*rad)
-0.007910*epsio*sin((msanom-mmanom+2*d)*rad)
-0.006783*epsio*sin((2*d+msanom)*rad)
+0.005162*sin((mmanom-d)*rad)
+0.005000*epsio*sin((msanom+d)*rad)
+0.004049*epsio*sin((mmanom-msanom+2*d)*rad)
+0.003996*sin((2*mmanom+2*d)*rad)
+0.003862*sin(4*d*rad)
+0.003665*sin((2*d-3*mmanom)*rad)
+0.002695*epsio*sin((2*mmanom-msanom)*rad)
+0.002602*sin((mmanom-2*f-2*d)*rad)
+0.002396*epsio*sin((2*d-msanom-2*mmanom)*rad)
-0.002349*sin((mmanom+d)*rad)
+0.002249*epsio*epsio*sin((2*d-2*msanom)*rad)
-0.002125*epsio*sin((2*mmanom+msanom)*rad)
-0.002079*epsio*epsio*sin(2*msanom*rad)
+0.002059*epsio*epsio*sin((2*d-mmanom-2*msanom)*rad)
-0.001773*sin((mmanom+2*d-2*f)*rad)
-0.001595*sin((2*f+2*d)*rad)
+0.001220*epsio*sin((4*d-msanom-mmanom)*rad)
-0.001110*sin((2*mmanom+2*f)*rad)
+0.000892*sin((mmanom-3*d)*rad)
-0.000811*epsio*sin((msanom+mmanom+2*d)*rad)
+0.000761*epsio*sin((4*d-msanom-2*mmanom)*rad)
+0.000717*epsio*epsio*sin((mmanom-2*msanom)*rad)
+0.000704*epsio*epsio*sin((mmanom-2*msanom-2*d)*rad)
+0.000693*epsio*sin((msanom-2*mmanom+2*d)*rad)
+0.000598*epsio*sin((2*d-msanom-2*f)*rad)
+0.000550*sin((mmanom+4*d)*rad)
+0.000538*sin(4*mmanom*rad)
+0.000521*epsio*sin((4*d-msanom)*rad)
+0.000486*sin((2*mmanom-d)*rad);
while lambda > 360.0 do
lambda := lambda -360.0;
while lambda < 0.0 do
lambda := lambda + 360.0;
{ now for Beta }
beta := +5.128189*sin(f*rad)
+0.280606*sin((mmanom+f)*rad)
+0.277693*sin((mmanom-f)*rad)
+0.173238*sin((2*d-f)*rad)
+0.055413*sin((2*d+f-mmanom)*rad)
+0.046272*sin((2*d-f-mmanom)*rad)
+0.032573*sin((2*d+f)*rad)
+0.017198*sin((2*mmanom+f)*rad)
+0.009267*sin((2*d+mmanom-f)*rad)
+0.008823*sin((2*mmanom-f)*rad)
+0.008247*epsio*sin((2*d-msanom-f)*rad)
+0.004323*sin((2*d-f-2*mmanom)*rad)
+0.004200*sin((2*d+f+mmanom)*rad)
+0.003372*epsio*sin((f-msanom-2*d)*rad)
+0.002472*epsio*sin((2*d+f-msanom-mmanom)*rad)
+0.002222*epsio*sin((2*d+f-msanom)*rad)
+0.002072*epsio*sin((2*d-f-msanom-mmanom)*rad)
+0.001877*epsio*sin((f-msanom+mmanom)*rad)
+0.001828*sin((4*d-f-mmanom)*rad)
-0.001803*epsio*sin((f+msanom)*rad)
-0.001750*sin(3*f*rad)
+0.001570*epsio*sin((mmanom-msanom-f)*rad)
-0.001487*sin((f+d)*rad)
-0.001481*epsio*sin((f+msanom+mmanom)*rad)
+0.001417*epsio*sin((f-msanom-mmanom)*rad)
+0.001350*epsio*sin((f-msanom)*rad)
+0.001330*sin((f-d)*rad)
+0.001106*sin((f+3*mmanom)*rad)
+0.001020*sin((4*d-f)*rad)
+0.000833*sin((f+4*d-mmanom)*rad)
+0.000781*sin((mmanom-3*f)*rad)
+0.000670*sin((f+4*d-2*mmanom)*rad)
+0.000606*sin((2*d-3*f)*rad)
+0.000597*sin((2*d+2*mmanom-f)*rad)
+0.000492*epsio*sin((2*d+mmanom-msanom-f)*rad)
+0.000450*sin((2*mmanom-f-2*d)*rad)
+0.000439*sin((3*mmanom-f)*rad)
+0.000423*sin((f+2*d+2*mmanom)*rad)
+0.000422*sin((2*d-f-3*mmanom)*rad)
-0.000367*epsio*sin((msanom+f+2*d-mmanom)*rad)
-0.000353*epsio*sin((msanom+f+2*d)*rad)
+0.000331*sin((f+4*d)*rad)
+0.000317*epsio*sin((2*d+f-msanom+mmanom)*rad)
+0.000306*epsio*epsio*sin((2*d-2*msanom-f)*rad)
-0.000283*sin((mmanom+3*f)*rad);
beta := beta * (1.0 - 0.0004664*cos(omeg*rad)
- 0.0000754*cos((omeg+275.05-2.30*t)*rad));
end; {end of moonpos}
procedure illum; { this procedure calculates the percent of
the moon's disk that is illuminated}
var
d,i,k,ang,y : real;
begin {start of illum}
ang := cos((lambda-theta)*pi/180.0)*cos(beta*pi/180.0);
{ The line solves for the arccos. i.e
the governing equation is
cos(d)=ang }
d := ArcTan(sqrt(1-ang*ang)/ang);
d := d*180.0/pi;
if d < 0. then
d := d + 180.0;
i := 180.0-d-0.1468*
(1-0.0549*sin(mmanom*pi/180.0))/(1-0.0167*sin(msanom*pi/180.0))
*sin(d*pi/180.0);
k := (1+cos(i*pi/180.0))/2.0;
k := k*100.0;
writeln('The Moon`s Disk will be ',k:3:0,' percent illuminated tonight.');
writeln;
end; {end of illum}
procedure equator; {This procedure converts the eclipetal lunar
coordinates to equatorial coordinates
epoch 2000}
var
a,y,RA,DEC,rad,hr,min,sec,deg : real;
const
epsi = 23.4392911;
begin {start of equator}
rad := pi/180.0;
y := (sin(lambda*rad)*cos(epsi*rad)-(sin(beta*rad)/cos(beta*rad))
*sin(epsi*rad));
a := y/cos(lambda*rad);
RA := ArcTan(a)/rad;
if (a > 0.0) and (y > 0.0) and (RA < 0.0) then
RA := RA + 90.0;
if (a > 0.0) and (y < 0.0) then
begin
if RA > 0.0 then
RA := RA + 180.0 else
RA := RA + 270.0;
end;
if (a < 0.0) and (y > 0.0) then
begin
if RA > 0.0 then
RA := RA + 90.0 else
RA := RA + 180.0;
end;
if (a< 0.0) and (y < 0.0) then
begin
if RA > 0.0 then
RA := RA + 270.0 else
RA := RA + 360.0;
end;
RA := RA/15.0; {convert to hours from degrees}
mra := RA;
writeln(' The Moon`s Right Ascension is:');
hr := Int(RA);
min := Int(Frac(RA)*60.0);
sec := Int(Frac(Frac(RA)*60.0)*60.0);
writeln(' ',hr:2:0,' hours ',min:2:0,' min ',sec:2:0,' seconds');
writeln;
writeln(' The Moons Declination is:');
a := sin(beta*rad)*cos(epsi*rad) + cos(beta*rad)*sin(epsi*rad)*
sin(lambda*rad);
{the following equation needs to be solved:
sin(DEC)=a
the next line finds the arcsin using the ArcTan function.}
DEC := ArcTan(a/sqrt(1-a*a))/rad;
mdec := DEC;
deg := Int(DEC);
min := Int(Frac(DEC)*60.0);
sec := Int(Frac(Frac(DEC)*60.0)*60.0);
writeln(' ',deg:3:0,' degrees ',min:2:0,' min ',sec:2:0,' seconds');
end; {end of equator}
{insert moonplot here}
begin { start of moon}
continue := true;
while continue do
begin
start;
julian;
writeln ('julian date of the observation is ',jd:10:3);
julian_new;
phaseout;
sunlat;
moonpos;
writeln;
illum;
equator;
writeln;
write('Do you want to check another date (y/n)? ');
readln(answer);
if (answer) <> 'y' then continue := false;
end; {of while}
clrscr;
writeln;
writeln;
writeln('Astronomy is a beautiful adventure! I hope this program helps you');
writeln('obtain more pleasure from your investigation of the night skies.');
writeln('This program is Copyrighted. However, the author grants everyone ');
writeln('the right to copy and distribute this program as long as they');
writeln('do not charge money or services in exchange for the program.');
writeln;
writeln('If you find this program useful and wish to encourage the ');
writeln('development of other free programs, the author requests a ');
writeln('donation be sent to:');
writeln;
writeln(' Fred Mendenhall');
writeln(' 2209 Tam-O-Shanter Ct.');
writeln(' Carmel, Ind. 46032');
writeln;
writeln('Copies of the Pascal com file and star plotting routines are');
writeln(' available by sending a SASE and blank disk.');
writeln;
writeln(' Happy Observing!');
end. {end of moon}