Sunrise / Sunset Algorithm Implemented in Free Pascal / Lazarus

I needed to know sunrise / sunset times and wasn’t happy with any Pascal implementations I found. This code was implemented using the algorithm found at

http://williams.best.vwh.net/sunrise_sunset_algorithm.htm

This is a simple unit that can be easily implemented. There are no O/S restrictions (that I know of), so it should work anywhere. I’ve run it on Windows and Raspberry Pi platforms successfully.

The call is:

time := sunRiseSet(sunMode, date, latitude, longitude, zenithMode, localOffset);

Where sunmode is

    sunrising                   = 0;
    sunsetting                  = 1;

Date is the date in question.

latitude,longitude are the fractional coordinates of the location in question.

ZenithMode is

    zenithOfficial              = 0;
    zenithCivil                 = 1;
    zenithNautical              = 2;
    zenithAstronomical          = 3;

and localOffset is the adjustment to convert the time from UTC time to local time (example -8 for US, West Coast).

unit sunRiseSetUnit;

{$mode objfpc}{$H+}{$R+}{$I+}

//-----------------------------------------------------------------------------
//
//  Sunrise/Sunset Time module
//
// 1.0 Apr 2015:
//      Original Version:
//
// ----------------------------------------------------------------------------

interface

uses
    classes,
    dateUtils;

function sunRiseSet(
    sunMode                     : integer;
    date                        : TDateTime;
    latitude                    : double;
    longitude                   : double;
    zenithMode                  : integer;
    localOffset                 : real
    )                           : TDateTime;

// ----------------------------------------------------------------------------

implementation

uses
    math,
    sysUtils;

function sunRiseSet(
    sunMode                     : integer;
    date                        : TDateTime;
    latitude                    : double;
    longitude                   : double;
    zenithMode                  : integer;
    localOffset                 : real
    )                           : TDateTime;

const
    sunrising                   = 0;
    sunsetting                  = 1;

    zenithOfficial              = 0;
    zenithCivil                 = 1;
    zenithNautical              = 2;
    zenithAstronomical          = 3;

    zenithDeg                   : array [0 .. 3] of double = (
                                    90.833333,
                                    96.0,
                                    102.0,
                                    108.0
                                    );

var
    cosDec                      : double;
    cosH                        : double;
    day                         : integer;
    h                           : double;
    l                           : double;
    lQuadrant                   : double;
    lngHour                     : double;
    m                           : double;
    month                       : integer;
    n                           : integer;
    n1                          : integer;
    n2                          : integer;
    n3                          : integer;
    oldDateSeparator            : char;
    oldShortDateFormat          : string;
    ra                          : double;
    raQuadrant                  : double;
    sinDec                      : double;
    t                           : double;
    tLocal                      : double;
    tLocalMean                  : double;
    tUTC                        : double;
    year                        : integer;

    // the algorith requires degrees NOT radians. These functions are
    // wrappers for the various trig functions which convert degrees to radians
    function degACOS(
            n                   : double
            )                   : double;
    begin
    degACOS := (180/pi)*arccos(n);
    end; // degACOS
    function degASIN(
            n                   : double
            )                   : double;
    begin
    degASIN := (180/pi)*arcsin(n);
    end; // degASIN
    function degATAN(
            n                   : double
            )                   : double;
    begin
    degATAN := (180/pi)*arctan(n);
    end; // degATAN
    function degCOS(
            angle               : double
            )                   : double;
    begin
    degCOS := cos((pi/180)*angle);
    end; // degCOS
    function degSIN(
            angle               : double
            )                   : double;
    begin
    degSin := sin((pi/180)*angle);
    end; // degSIN
    function degTAN(
            angle               : double
            )                   : double;
    begin
    degTAN := tan((pi/180)*angle);
    end; //degTAN

begin

if (not(sunMode in [0,1])) then begin
    writeln('sunMode must be 0-1.');
    halt(1);
    end;
if (latitude < -90) or (latitude > 90) then begin
    writeln('latitude outside of -90-90 range.');
    halt(1);
    end;
if (longitude < -180) or (longitude > 180) then begin
    writeln('longitude outside of -180-180 range.');
    halt(1);
    end;
if (not(zenithMode in [0..4])) then begin
    writeln('zenithMode must be 0..3.');
    halt(1);
    end;

month   := monthof(date);
day     := dayof(date);
year    := yearof(date);

// Calc day of year
n1  := floor(275 * month / 9);
n2  := floor((month + 9) / 12);
n3  := (1 + floor((year - 4 * floor(year / 4) + 2) / 3));
n   := n1 - (n2 * n3) + day - 30;

// Convert longitude to hour and calc approximate time
lngHour := longitude / 15;

if sunMode = sunRising then
    t := n + ((6 - lngHour) / 24)
else
    t := n + ((18 - lngHour) / 24);

// Calc the Sun's Mean anomaly
m := (0.9856 * t) - 3.289;

// Calc Sun's true longitude
l := m + (1.916 * degSIN(m)) + (0.020 * degSIN(2*M)) + 282.634;
if (l < 0) then
    l := l + 360
else if (l > 360) then
    l := l - 360;

// Calc the sun's right ascension
ra := degATAN(0.91764 * degTAN(L));
if (ra < 0) then
    ra := ra + 360
else if (ra > 360) then
    ra := ra - 360;

// Right ascension needs to be in the same quadrant as L
lQuadrant := (floor(l / 90)) * 90;
raQuadrant := (floor(ra / 90)) * 90;
ra := ra + (lQuadrant - raQuadrant);

// Convert Right ascension into hours
ra := ra / 15;

// Calc the Sun's declination
sinDec := 0.39782 * degSIN(l);
cosDec := degCOS(degASIN(sinDec));

// Calc the Sun's local hour angle
cosH := (degCOS(zenithDeg[zenithMode]) - (sinDec * degSIN(latitude))) /
            (cosDec * degCOS(latitude));
if (cosH > 1) then begin
    // Sun doesn't rise at this location on this date. We will just
    // set it to 12:00
    result := strToDateTime(format('%2d/%2d/%2d',[month, day, year]) + ' ' +
                '12:00');
    end
else if (cosH < -1) then begin
    result := strToDateTime(format('%2d/%2d/%2d',[month, day, year]) + ' ' +
                '12:00');
    end;

// Finish calculating H and convert into hours
if sunMode = sunRising then
    h := 360 - degACOS(cosH)
else
    h := degACOS(cosH);
h := h / 15;

// Calculate local mean time of rising/setting
tLocalMean := h + ra - (0.06571 * t) - 6.622;

// adjust back to UTC
tUTC   := tLocalMean - lngHour;

// adjust back to local time
tLocal := tUTC + localOffset;
if (tLocal < 0) then
    tLocal := tLocal + 24
else if (tLocal > 24) then
    tLocal := tLocal - 24;

// if minutes would round up to 60, go to next hour:
if (round((tLocal-floor(tLocal))*60) = 60) then
    tLocal := floor(tLocal) + 1;

oldShortDateFormat  := shortDateFormat;
oldDateSeparator    := dateSeparator;
shortDateFormat     := 'm/d/y';
dateSeparator       := '/';
result              := strToDateTime(format('%2D/%2D/%4D',[month, day, year]) + ' ' +
                        inttostr(floor(tLocal)) + ':' +
                        inttostr(round((tLocal-floor(tLocal))*60)));
shortDateFormat     := oldShortDateFormat;
dateSeparator       := oldDateSeparator;

end; //sunRiseSet

end.
This entry was posted in c-lazarus, c-pcos and tagged . Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s