Category : Science and Education
Archive   : HC21DEMO.ZIP
Filename : SUNTIMES.PAS

{==================================================================
This program demonstrates how to share the values in System and User
Variables in an HC2000 User Program with a foreground program that
has been executed by the DOS Shell program "HCDOS.EXE".

The HC200 Interface unit (HCIF.TPU) provides several functions and
procedures, two of which (GetVariable, PutVariable) are used to access
the User Program Variables by matching the spelling and case exactly.
A third procedure (HC2000Access) must be called first to initialize
the interface.

The program computes the rise and set times (in minutes after midnight)
of the Sun for a specific longitude and latitude any date after 1990.

The following is derived from the worked examples in the book:-

"Practical Astronomy with your Calculator"
third edition by PETER DUFFET-SMITH.

==================================================================}


program Suntimes;

uses
crt,
HCIF; { HC2000 Interface Unit }

var
PIby180, BigReal, epsilon: Real;
Ra,De,GSTr,GSTs,GEL:Real;
UTr, UTs: Real;
SunRise,SunSet,GMToffset: Integer;
Date,Month,Year: integer;


procedure WriteTime(T:Real);

var
V: longInt;

begin
V := TRUNC(T);
write(V:2,'h ');
T := FRAC(T) * 60;

V := TRUNC(T);
write(V:2,'m ');
T := FRAC(T) * 60;

write(ROUND(T):2,'s ');
end;


function Sign(R: Real):ShortInt;

begin
Sign := 1 - 2 * ord(R < 0);
end;


function SinD(A:Real):Real;

begin
SinD := Sin(A * PIby180);
end;


function CosD(A:Real):Real;

begin
CosD := Cos(A * PIby180);
end;


function TanD(A:Real):Real;

var
T: Real;

begin
T := CosD(A);
if abs(A) < epsilon then
TanD := Sign(T) * Sign(SinD(A)) * BigReal
else TanD := SinD(A) / T;
end;


function AtanD(A:Real):Real;

begin
AtanD := ArcTan(A) / PIby180;
end;


function HMStoDecimalHours(H,M:integer; S:Real):Real;

begin
HMStoDecimalHours := ((S / 60) + M) / 60 + H;
end;



function DMStoDecimalDegrees(D,M,S:integer):Real;

begin
DMStoDecimalDegrees := ((s / 60) + M) /60 + D;
end;


function AsinD(V:Real):Real;

begin
if 1 - V < epsilon then AsinD := 90
else if V + 1 < epsilon then AsinD := - 90
else AsinD := ArcTan(V / Sqrt(1 - Sqr(V))) / PIby180;
end;


function AcosD(V:Real):Real;

begin
if abs(V) < epsilon then AcosD := PI / 2
else AcosD := (ArcTan(Sqrt(1 - Sqr(V)) / V) + PI * ord(V < 0)) / PIby180;
end;


function LSTtoGST(LST,Longitude:Real): Real;

var
D: Real;

begin
D := Longitude / 15;
D := LST + D;
if D > 24 then D := D - 24;
if D < 0 then D := D + 24;
LSTtoGST := D;
end;


function JulianDay(Year,Month:integer; Day:Real):Real;

var
Y1,M1: Real;
A, B, C, D: Real;

begin
Y1 := Year;
M1 := Month;
if Month < 3 then begin
Y1 := Year - 1;
M1 := Month + 12;
end;
A := INT(Y1 / 100);
B := 2 - A + INT(A / 4);
C := INT(365.25 * Y1);
D := INT(30.6001 * (M1 + 1));
JulianDay := B + C + D + Day + 1720994.5;
end;


function GSTtoUT(Y,M,D:integer;GST: real):Real;

var
S, T, X, T0: Real;

begin
S := JulianDay(Y,M,D) - 2451545.0;
T := S / 36525.0;
T0 := 6.697374558 + (2400.051336 * T) + (0.000025862 * T * T);
X := GST - T0;
while X > 24 do X := X - 24;
while X < 0 do X := X + 24;
GSTtoUT := X * 0.9972695663;
end;


procedure RiseAndSeting(Y, M, D:integer;
RA, DE: Real;
LA, LO: Real;
var GSTr, GSTs: Real);

var
Ar,As: Real;
CosAr,H: Real;
LSTr,LSTs: Real;

begin
CosAr := SinD(DE)/CosD(LA);
Ar := AcosD(CosAr);

H := 1 / 15 * AcosD(0 - TanD(LA) * TanD(DE));
LSTr := 24 + RA - H;
if LSTr > 24 then LSTr := LSTr - 24;
LSTs := RA + H;
if LSTs > 24 then LSTs := LSTs - 24;

GSTr := LSTtoGST(LSTr,LO);
GSTs := LSTtoGST(LSTs,LO);
end;


function DateToDayNumber(Y,M,D:integer):integer;
var
X, YT: integer;
begin
if ((Y MOD 4) = 0) and (true) then
YT := 62
else YT := 63;

if M > 2 then
X := TRUNC((M + 1) * 30.6) - YT
else
X := TRUNC(((M - 1) * YT) div 2);
DateToDayNumber := X + D;
end;


procedure GELtoRaDe(GElo,GEla:Real;
var Ra,De:Real);
var
e, SinDe, x, y, ap, M: Real;
begin
e := 23.441884;
SinDe := SinD(GEla) * CosD(e) + CosD(GEla) * SinD(e) * SinD(GElo);
De := AsinD(SinDe);
y := SinD(GElo) * CosD(e) - TanD(GEla) * SinD(e);
x := CosD(GElo);
ap := AtanD(y / x);
if y > 0 then
if x > 0 then M := 90
else M := 180
else if x > 0 then M := 360
else M := 270;
while ap > M do ap := ap - 180;
while ap < (M - 90) do ap := ap + 180;
Ra := ap / 15;
end;


function UTtoGST(Y,M,D:integer; UT:Real):Real;
var
X,JD,S,T,T0:Real;
begin
JD := JulianDay(Y,M,D);
S := JD - 2451545.0;
T := S / 36525.0;
T0 := 6.697374558 + 2400.051336 * T + 0.000025862 * sqr(T);
X := UT * 1.002737909 + T0;
while X < 0 do X := X + 24;
while X > 24 do X := X - 24;
UTtoGST := X;
end;


procedure GetSunRaDe(Y,M,D:integer; var GEL,Ra,De:Real);
var
Dx,N,Msun,Ec:Real;
Bias: integer;
begin
Bias := trunc((Y - 1990) * 365.24 + 0.5);

Dx :=DateToDayNumber(Y,M,D) + Bias;
N := 360.0 / 365.242191 * Dx;
while N > 360 do N := N - 360;
while N < 0 do N := N + 360;
Msun := N + 279.403303 - 282.768422;
if Msun < 0 then Msun := Msun + 360;
Ec := 360.0 / PI * 0.016713 * SinD(Msun);
GEL := N + Ec + 279.403303;
if GEL > 360 then GEL := GEL - 360;
GELtoRaDe(GEL,0.0,Ra,De);
end;


procedure SunRiseAndSet(Y,M,D:integer;
LA, LO: Real;
var UTr, UTs: Real);

var
GSTr,GSTs,GST1r,GST1s,GST2r,GST2s:Real;
T0,T0p,Ra,De:Real;

begin
{Sun Coords midnight this day}
GetSunRaDe(Y,M,D, GEL,Ra,De);

{Rise and set of this position}
RiseAndSeting(Y,M,D,Ra,De,LA,LO,GST1r,GST1s);

{Sun Coords midnight next day}
GetSunRaDe(Y,M,D + 1, GEL,Ra,De);

{Rise and set of this position}
RiseAndSeting(Y,M,D + 1,Ra,De,LA,LO,GST2r,GST2s);

T0 := UTtoGST(Y,M,D,0.0);
T0p := T0 - (LO / 15.0 * 1.002738);
if T0p < 0 then T0p := T0p + 24;

if GST1r < T0p then begin
GST1r := GST1r + 24;
GST2r := GST2r + 24;
end;

if GST1s < T0p then begin
GST1s := GST1s + 24;
GST2s := GST2s + 24;
end;

GSTr := (24.07 * GST1r - T0 * (GST2r - GST1r)) / (24.07 + GST1r - GST2r);
GSTs := (24.07 * GST1s - T0 * (GST2s - GST1s)) / (24.07 + GST1s - GST2s);

GSTr := GSTr - 0.075319;
GSTs := GSTs + 0.075319;

UTr := GSTtoUT(Y,M,D,GSTr);
UTs := GSTtoUT(Y,M,D,GSTs);
end;


function LocalTimeInMinutes(X,Zone:Real):integer;

var
Y:Real;

begin
Y := X + Zone;
if Y > 24 then Y := Y - 24;
if Y < 0 then Y := Y + 24;
LocalTimeInMinutes := round(Y * 60);
end;


begin {Main}
writeln;
writeln('**************** SunTimes calculation begins ****************');

epsilon := 1E-100; {Very small number}
BigReal := 1E2+40; {Very big number}
PIby180 := PI / 180;


{ Initialize HC Interface so that Variables can be accessed in the
HC2000 User Program }
AccessHC2000;

{ Get Variables from User Program. The Names must match exactly. }

Year := GetVariable('Year');
Month := GetVariable('Month');
Date := GetVariable('Date');
GMToffset := GetVariable('Time Zone Code'); {for Ottawa:-
winter .... 5
summer .... 4 }

{=====================================
CITY Map coordinate

Ottawa, Can 45.5 73.5
Toronto, Can 43.5 79.5
New York, NY 40.75 74.0
Los Angeles, CA 34.0 118.3
London, UK 51.5 0.0
Sydney, Aus -33.9 -151.0
=====================================}

{ Compute Sunrise and Sunset times (for Ottawa, Canada) }
SunRiseAndSet(Year,Month,Date, 45.5, 73.5, UTr, UTs);

write(Month,'/',Date,'/',Year,' Sunrise and Sunset (GMT) ');
WriteTime(UTr);
WriteTime(UTs);
writeln;

SunRise := LocalTimeInMinutes(GMToffset,UTr);
SunSet := LocalTimeInMinutes(GMToffset,UTs);

writeln('Local Sunrise & Sunset in minutes after midnight:- ',
SunRise,', ', SunSet);

{ Return results to User Variables in User Program }

PutVariable('Sunrise time in minutes', SunRise);
PutVariable('Sunset time in minutes', SunSet);

writeln('**************** SunTimes calculation ends ******************');
writeln;
end.