In order to make quite comfortable input possible the source code is somewhat longer than necessary. The following items are to be noticed:
program ParallaxOfTheSun;
uses Crt;
const deg = Pi/180.0;
type real = extended; (* highest precision is necessary! *)
var name: String; (* name of the observed object *)
day, month, year: integer; (* date of the observations *)
UT1, UT2, (* observation times in UT *)
long1, lat1, h1, long2, lat2, h2: real;(* observer's geogr. coordinates *)
alpha1, delta1, d1, alpha2, delta2, d2, (* geocentric coordinates *)
alpha11, delta11, alpha21, delta21, (* topocentric coordinates *)
dist, parallaxe, dmin: real;
geozC, key: Char;
geoz: boolean;
dayStr, monthStr, yearStr: String; (* strings for comfortable input *)
UT1Str, UT2Str: String;
alpha1Str, delta1Str, alpha11Str, delta11Str, d1Str: String;
alpha2Str, delta2Str, alpha21Str, delta21Str, d2Str: String;
lat1Str, long1Str, h1Str, lat2Str, long2Str, h2Str: String;
function hmin2hdec(hmin: real): real;
(* ... convertes hours with minutes into hours with decimals *)
begin
hmin2hdec:=Trunc(hmin)+(hmin-Trunc(hmin))/0.6
end;
function String2Angle(s: String): real;
(* ... converts a string into an angle in rad. The must have one of the
following forms:
- degree with decimals XXX.YYYY, for instance 35.134,
- hours with minutes and seconds XXhYYmZZ.ZZZ, for instance 3h56m5.03,
- degrees with arcmin and arcsec XXdYYmZZ.ZZZ, for instance 25d35m28.33,
- degrees with arcmin and arcsec XXgYYmZZ.ZZZ, for instance 25g35m28.33.*)
var pd, ph, code, sgn: integer;
r, r1: real;
s1: String;
begin
while (s[1]=' ') do Delete(s, 1, 1);
while Pos(',', s)>0 do s[Pos(',', s)]:='.';
if s[1]='-'
then
begin
sgn:=-1; Delete(s, 1, 1);
end
else sgn:=1;
pd:=Pos('d', s); ph:=Pos('h', s);
if (pd=0) and (ph=0)
then
begin
Val(s, r, code);
if Code>0 then begin Write(#7); Halt end;
end
else
if pd>0
then
begin
s1:=Copy(s, 1, pd-1); Delete(s, 1, pd);
Val(s1, r, code);
pd:=Pos('m', s);
if pd=0 then pd:=Pos(#39, s);
s1:=Copy(s, 1, pd-1); Delete(s, 1, pd);
Val(s1, r1, code);
r:=r+r1/60.0;
Val(s, r1, code);
r:=r+r1/3600.0;
end
else
begin
s1:=Copy(s, 1, ph-1); Delete(s, 1, ph);
Val(s1, r, code);
ph:=Pos('m', s);
s1:=Copy(s, 1, ph-1); Delete(s, 1, ph);
Val(s1, r1, code);
r:=r+r1/60.0;
Val(s, r1, code);
r:=(r+r1/3600.0)*15.0;
end;
String2Angle:=sgn*r*deg
end;
procedure readAngle(s: String; var sw: String; var w: real);
(* quite comfortable input routine for angles: The prechoosen value is
displayed. Pressing "return" only will keep the old value unchanged. *)
var sh: String;
Code: integer;
begin
Write(s+'('+sw+') '); ReadLn(sh);
if sh<>'' then sw:=sh;
w:=String2Angle(sw)
end;
procedure readReal(s: String; var sr: String; var r: real);
(* quite comfortable input routine for reals, similar to "readAngle" *)
var sh: String;
Code: integer;
begin
Write(s+'('+sr+') '); ReadLn(sh);
if sh<>'' then sr:=sh;
Val(sr, r, Code);
end;
procedure readInteger(s: String; var si: String; var i: integer);
(* quite comfortable input routine for integers, similar to "readAngle" *)
var sh: String;
Code: integer;
begin
Write(s+'('+si+') '); ReadLn(sh);
if sh<>'' then si:=sh;
Val(si, i, Code);
end;
procedure polar2cart(alpha, delta, r: real; var x, y, z: real);
(* ... converts polar coordinates in rectangular coordinates *)
var rcst: real;
begin
rcst:=r*cos(delta);
x:=rcst*cos(alpha);
y:=rcst*sin(alpha);
z:=r*sin(delta);
end;
function julDate(day, month, year: integer; UT: real): real;
(* ... calculates the Julian date according to normal date and time
(in h.dec UT) *)
var t, k, b, y, m: real; a: integer;
begin
t:=day+UT/24.0;
k:=10000.0*year+100.0*month+t;
b:=-63.5;
y:=year+4712.0;
m:=month+1.0;
if m<=3.0 then
begin
y:=y-1.0; m:=m+12.0
end;
if k>=15821015.0 then
begin
a:=Trunc((y+88.0)/100.0);
b:=b+38.0-a+Trunc(a/4.0)
end;
t:=Trunc(365.25*y)+Trunc(30.6001*m)+t+b;
julDate:=t
end;
function sideralTime(day, month, year: integer; UT, longitude: real): real;
(* ... calculates the sideral time for a moment and a site with
geographical "longtude" (in rad; east of Greenwich: >0!) *)
const ST2UT = 1.0027379093;
var s, j: real;
begin
j:=julDate(day, month, year, 0.0);
s:=6.656306+0.0657098242*(j-2445700.5)+longitude*12.0/Pi;
s:=s+ST2UT*UT;
while s>24.0 do s:=s-24.0;
while s< 0.0 do s:=s+24.0;
sideralTime:=s/12.0*Pi;
end;
procedure ObserversCoord(day, month, year: integer;
ut, longitude, latitude, height: real;
var alpha, delta, rho: real);
(* ... calculates the equitorial coordinates of an observer for one moment,
takes into account that the earth is no perfect sphere. Latitude>0 means:
east of Greenwich! ( see J. Meeus, Astronomical Algorithms, for instance) *)
const exEarth = 0.08181922; (* excentricity of the earth's body *)
var denum, sinl, cosl, rhosin, rhocos: real;
begin
sinl:=sin(latitude); cosl:=cos(latitude);
denum:=sqrt(1.0-sqr(exEarth*sinl));
rhocos:=cosl/denum+height/6378140.0*sinl;
rhosin:=(1.0-exEarth*exEarth)*sinl/denum+height/6378140.0*cosl;
rho:=sqrt(sqr(rhosin)+sqr(rhocos));
delta:=arctan(rhosin/rhocos);
alpha:=sideraltime(day, month, year, ut, longitude);
end;
procedure initialisation;
(* presetting of the input variables, allowing a test run by pressing
"return" only *)
begin
lat1Str:= '52.25417'; long1Str:='8.3244'; h1Str:='200.0';
lat2Str:= '52.25417'; long2Str:='8.3244'; h2Str:='200.0';
UT1Str:='19.1695'; UT2Str:='20.1002';
dayStr:='14'; monthStr:='10'; yearStr:='1996';
alpha1Str:='2h42m47.518'; delta1Str:='33d9m6.54';
alpha11Str:='2h42m47.729'; delta11Str:='33d9m1.332'; d1Str:='1.028148';
alpha2Str:='2h42m45.377'; delta2Str:='33d9m11.67';
alpha21Str:='2h42m45.54'; delta21Str:='33d9m6.732'; d2Str:='1.028072';
end;
procedure input;
begin
ClrScr;
WriteLn('This program calculates the sun'+#39+'s parallax from');
WriteLn
('- two simult. topocentric positions and the corresp. geocentric distance,');
WriteLn
('- two sequ. topocentric positions and the corresp. geocentric proper motion,');
WriteLn
('- one topocentric position and corresponding geocentric coordinates.');
WriteLn;
Write(' Asteroid'+#39+'s name : '); ReadLn(name);
Write(' Has '+name+' been observed twice (y/n) ? '); ReadLn(geozC);
geoz:=geozC in ['n','N'];
WriteLn;
WriteLn(' Date:');
readInteger(' Day: ', dayStr, day);
readInteger(' Month: ', monthStr, month);
readInteger(' Year: ', yearStr, year);
WriteLn('geographic coordinates of observer 1:');
readAngle(' latitude = ', lat1Str, lat1);
readAngle(' longitude (east of Greenwich: >0) = ', long1Str, long1);
readReal(' height [meters over sea level] = ', h1Str, h1);
WriteLn;
readReal(' observation time 1 (h.min UT) = ', UT1Str, UT1);
UT1:=hmin2hdec(UT1);
if not geoz then
begin
readReal(' observation time 2 (h.min MEZ) = ', UT2Str, UT2);
UT2:=hmin2hdec(UT2);
end;
WriteLn;
WriteLn(' topocentric coordinates:');
readAngle(' longitude = ', alpha11Str, alpha11);
readAngle(' latitude = ', delta11Str, delta11);
WriteLn(' geocentric coordinates:');
if UT1<>UT2 then
begin
readAngle(' longitude = ', alpha1Str, alpha1);
readAngle(' latitude = ', delta1Str, delta1);
end;
readReal(' distance [AU] = ', d1Str, d1);
if not geoz then
begin
WriteLn;
WriteLn(' geographic coordinates of observer 2:');
readAngle(' latitude = ', lat2Str, lat2);
readAngle(' longitude (east of Greenwich: >0) = ', long2Str, long2);
readReal(' height [meters over sea level] = ', h2Str, h2);
WriteLn;
WriteLn(' topocentric coordinates:');
readAngle(' longitude = ', alpha21Str, alpha21);
readAngle(' latitude = ', delta21Str, delta21);
WriteLn(' geocentric coordinates:');
if UT1<>UT2 then
begin
readAngle(' longitude = ', alpha2Str, alpha2);
readAngle(' latitude = ', delta2Str, delta2);
end;
readReal(' distance [AU] = ', d2Str, d2);
end;
end;
function arccos(x: real): real;
begin
if x>0.0
then arccos:=arctan(sqrt(1.0/sqr(x)-1.0))
else if x<0.0
then arccos:=-arctan(sqrt(1.0/sqr(x)-1.0))+Pi
else arccos:=Pi/2.0 end;
function distance(alpha1, delta1, alpha2, delta2: real): real;
(*... calculates the angular distance of two points of the celestial sphere *)
begin
distance:=arccos(sin(delta1)*sin(delta2)+
cos(delta1)*cos(delta2)*cos(alpha1-alpha2))
end;
procedure calculation(var dist, parallaxe, dmin: real);
(* the heart of this program *)
var ao1, do1, ro1, ao2, do2, ro2: real;
xo1, yo1, zo1, xo2, yo2, zo2: real;
xe1, ye1, ze1, xe2, ye2, ze2: real;
lpm, lmm, lambda, mu, r, AU: real;
begin
ObserversCoord(day, month, year, UT1, long1, lat1, h1, ao1, do1, ro1);
polar2cart(ao1, do1, ro1, xo1, yo1, zo1);
polar2cart(alpha11, delta11, 1.0, xe1, ye1, ze1);
if not geoz
then
begin
ObserversCoord(day, month, year, UT2, long2, lat2, h2, ao2, do2,ro2);
polar2cart(ao2, do2, ro2, xo2, yo2, zo2);
if UT1<>UT2 then
begin (* reduction of topoc. position no. 2: *)
alpha21:=alpha21-(alpha2-alpha1);
delta21:=delta21-(delta2-delta1);
end;
polar2cart(alpha21, delta21, 1.0, xe2, ye2, ze2);
end
else
begin
polar2cart(alpha1, delta1, 1.0, xe2, ye2, ze2);
xo2:=0.0; yo2:=0.0; zo2:=0.0;
end;
(* lambda+mu: *)
lpm:=((xo2-xo1)*(xe1-xe2)+(yo2-yo1)*(ye1-ye2)+(zo2-zo1)*(ze1-ze2))/
(1.0-(xe1*xe2+ye1*ye2+ze1*ze2));
(* lambda-mu: *)
lmm:=((xo2-xo1)*(xe1+xe2)+(yo2-yo1)*(ye1+ye2)+(zo2-zo1)*(ze1+ze2))/
(1.0+(xe1*xe2+ye1*ye2+ze1*ze2));
lambda:=0.5*(lpm+lmm);
mu:=0.5*(lpm-lmm);
(* shortest distance of the lines of view: *)
dmin:=sqrt(sqr(xo1+lambda*xe1-xo2-mu*xe2)+
sqr(yo1+lambda*ye1-yo2-mu*ye2)+
sqr(zo1+lambda*ze1-zo2-mu*ze2));
(* geocentric distance as multiple of the earth's radius: *)
r:=sqrt(sqr(xo1+lambda*xe1)+sqr(yo1+lambda*ye1)+sqr(zo1+lambda*ze1));
(* the Astronomical Unit as multiple of the earth's radius: *)
if lambda<0.0 then dmin:=-dmin;
if geoz
then AU:=r/d1
else AU:=r/(0.5*(d1+d2));
parallaxe:=1.0/AU;
dist:=distance(alpha11, delta11, alpha21, delta21);
end;
procedure output(var key: Char);
(* ... displays input data and results *)
begin
ClrScr;
Write('Two topocentric observations of ', name);
WriteLn(', ', month, '/', day, '/', year);
WriteLn('=============================================================');
WriteLn('Observation 1:');
WriteLn('--------------');
WriteLn('UT: ', UT1Str);
Write(' geographic position:');
Write(' latitude = ', lat1Str, ',');
Write(' longitude = ', long1Str, ',');
WriteLn(' height = ', h1Str);
Write('topocentric position:');
Write(' alpha = ', alpha11Str, ',');
WriteLn(' delta = ', delta11Str);
Write(' geocentric position:');
Write(' alpha = ', alpha1Str, ',');
WriteLn(' delta = ', delta1Str);
WriteLn(' r = ', d1Str, ' AU');
if not geoz then
begin
WriteLn;
WriteLn('Observation 2:');
WriteLn('--------------');
WriteLn('UT: ', Ut2Str);
Write(' geographic position:');
Write(' latitude = ', lat2Str, ',');
WriteLn(' longitude = ', long2Str, ',');
WriteLn(' height = ', h1Str);
Write('topocentric position:');
Write (' alpha = ', alpha21Str, ',');
WriteLn(' delta = ', delta21Str);
Write(' geocentric position:');
Write (' alpha = ', alpha2Str, ',');
WriteLn(' delta = ', delta2Str);
WriteLn(' r = ', d2Str, ' AU');
end;
WriteLn;
WriteLn('angle of parallax: ', dist/deg*3600.0:6:3, '"');
WriteLn('calculated parallax of the sun: ', parallaxe/deg*3600.0:6:3, '"');
WriteLn('=======================================');
WriteLn('minimal distance of the lines of view:', dmin:6:3, ' rE');
WriteLn; WriteLn;
Write(' once more (y/n) ? '); ReadLn(key);
end;
begin (* main program: *)
initialisation;
repeat
input;
calculation(dist, parallaxe, dmin);
output(key);
until key='n';
end.