Category : Science and Education
Archive   : DAYLIGHT.ZIP
Filename : CALCULAT.PAS

unit Calculate;

interface

uses Globals;

type

SunTime = record Dawn,Dusk:Float; end;

Daytime = record
Astronomical,
Nautical,
Civil,
Actual:SunTime;
end;

procedure Phenomina(lat,long,zone:float; y:longint;m,d:integer;
var Day:DayTime);


implementation

const
Sunup = pi/2;
Sundown = pi*3/2;


function Normalize(z:float):float;
begin
z := 2*pi * Frac(z/(2*pi));
if z<0 then
z := z + 2*pi;
Normalize := z;
end;

procedure CalcSolar(SunAt,lat,long:float; y:longint; m,d:integer;
var ApproxTime,
SolarRightAscension, SolarDeclination:float);
var
SolarMeanAnomaly,
SolarTrueLong
:Float;
n :integer;
begin
ApproxTime := DayOfYear(y,m,d) + (SunAt+long)/(2*pi);
SolarMeanAnomaly := ApproxTime * 0.017202 - 0.0574039;
SolarTrueLong := Normalize(
SolarMeanAnomaly
+ 0.0334405 * sin(SolarMeanAnomaly)
+ 3.49066E-4 * sin(2*SolarMeanAnomaly)
+ 4.93289
);
{Quadrant Determination}
if Frac(SolarTrueLong*2/pi) = 0.0 then
SolarTrueLong := SolarTrueLong + 4.84814e-6;
n := 2;
if SolarTrueLong <= pi/2 then
n := 0
else if SolarTrueLong <= pi*3/2 then
n := 1;
SolarRightAscension := ArcTan(0.91746 * Tan(SolarTrueLong));
{Quadrant Adjustment}
if n = 1 then
SolarRightAscension := SolarRightAscension + pi
else if n = 2 then
SolarRightAscension := SolarRightAscension + 2*pi;
{Solar Declination}
SolarDeclination := 0.39782 * Sin(SolarTrueLong);
SolarDeclination := ArcTan(SolarDeclination
/ Sqrt(1-Sqr(SolarDeclination))
);
end;

procedure ConvertCoords(r,q,lat:float; var s:float);
begin
s := (r-(sin(q)*sin(lat))) / (cos(q)*cos(lat));
end;

function SettingAdjustment(var SolarDeclination:float):float;
begin
SettingAdjustment := pi/2 - ArcTan( SolarDeclination
/ Sqrt(1-sqr(SolarDeclination))
);
end;

function RisingAdjustment(var SolarDeclination:float):float;
begin
RisingAdjustment := 2*pi - SettingAdjustment(SolarDeclination);
end;

procedure GetTime(long,zone:float; y,m,d: integer;
ApproxTime,SolarRightAscension, SolarDeclination:float;
var t:float);
begin
t := {Loval Apparent Time}
SolarDeclination + SolarRightAscension
- ApproxTime * 2*pi/365.2425 - 1.73364
{Universal Time}
+ long
{Wall Clock Time}
- zone;
t := Normalize(t) * 24{hours}/(2*pi{radians});

end;


procedure Phenomina(lat,long,zone:float; y:longint;m,d:integer;
var Day:DayTime);
var time,ascension,declination:float;
procedure RisingTime(r:float; var t:SunTime);
var s:float;
i:integer;
begin;
ConvertCoords(r,declination,lat,s);
if abs(s)>1.0 then
t.Dawn := 0
else
GetTime(long,zone,y,m,d,time,ascension,RisingAdjustment(s),t.Dawn);
end;
procedure SettingTime(r:float; var t:SunTime);
var s:float;
i:integer;
begin;
ConvertCoords(r,declination,lat,s);
if abs(s)>1.0 then
t.Dusk := 0
else
GetTime(long,zone,y,m,d,time,ascension,SettingAdjustment(s),t.Dusk);
end;
const
Astro = -0.309017;
Nautic = -0.207912;
Civil = -0.104528;
Actual = -0.0145439;
begin
lat := Radians(lat);
long := Radians(long);
zone := Radians(zone)*15.0; {15 degrees/Time Zone}

CalcSolar(Sunup,lat,long, y,m,d, time,ascension,declination);
RisingTime(Astro ,Day.Astronomical);
RisingTime(Nautic ,Day.Nautical);
RisingTime(Civil ,Day.Civil);
RisingTime(Actual ,Day.Actual);

CalcSolar(Sundown,lat,long, y,m,d, time,ascension,declination);
SettingTime(Actual ,Day.Actual);
SettingTime(Civil ,Day.Civil);
SettingTime(Nautic ,Day.Nautical);
SettingTime(Astro ,Day.Astronomical);
end;

end.