Just came across this, and thought I would share some code I modified to run on a PAC many years ago that I’ve been using for astronomical clock use. The math doesn’t seem to be too much for the PAC. I don’t remember how accurate it is, but I use it for turning off and on lights and it works great for that and it’s more reliable than photo cells.

I have this running in a subroutine with the following parameters:

fLatitude - The Latitude where you want the sunrise/sunset time.

fLongitude

nDayOfYear - The number of days since 12/31 last year (So Jan 1 = 1 and Feb 1 = 32, etc.)

fUTCOffset - Adjust the result into the local time zone

nSunset - This is the output in seconds since midnight for sunset

nSunrise - This is the output in seconds since midnight for sunrise

The rest are all local subroutine variables.

The code:

```
//Calculate Sunset given the latitude, longitude and day of year
//Adapted from http://williams.best.vwh.net/sunrise_sunset_algorithm.htm
//which came from the Almanac for Computers, 1990 USNO
//Define Pi and Degree and Radian conversion factors.
fPi = 3.14159265;
fDegToRad = .0174532925;
fRadToDeg = 57.2957795;
//Define Zenith (official sunset zenith)
fZenith = 90.83;
//Get the Longitude Hour
fLngHour = fLongitude / 15.0;
//Calculate an approximate sunset time for the longitude
fApproxTime = nDayOfYear + ((18.0 - fLngHour)/24.0);
//Sun Mean Anomaly
fSunMeanAnomaly = (.9856 * fApproxTime) - 3.289;
//Sun True Longitude
fSunTrueLongitude = fSunMeanAnomaly + (1.916 * Sine(fSunMeanAnomaly * fDegToRad)) + (.02 * Sine(2 * fSunMeanAnomaly * fDegToRad)) + 282.634;
//Make sure fSunTrueLongitude is between 0 and 360, adjust by 360.
if (fSunTrueLongitude>360.0) then
fSunTrueLongitude = fSunTrueLongitude - 360.0;
elseif (fSunTrueLongitude < 0.0) then
fSunTrueLongitude = fSunTrueLongitude + 360.0;
endif
//Sun Right Ascension (RA)
fSunRightAscension = Arctangent(0.91764 * Tangent(fSunTrueLongitude * fDegToRad)) * fRadToDeg + 360;
//Make sure fSunRightAscension is between 0 and 360, adjust by 360
if (fSunRightAscension>360.0) then
fSunRightAscension = fSunRightAscension - 360.0;
elseif (fSunRightAscension < 0.0) then
fSunRightAscension = fSunRightAscension + 360.0;
endif
//Put RA and L in the same quadrant
fSunTrueLongitudeQuadrant = Truncate(fSunTrueLongitude / 90.0) * 90.0;
fSunRAQuadrant = Truncate(fSunRightAscension / 90.0) * 90.0;
fSunRightAscension = fSunRightAscension + fSunTrueLongitudeQuadrant - fSunRAQuadrant;
//RA converted to Hours
fSunRAHours = fSunRightAscension / 15.0;
//Calc suns Declination
fSunSinDeclination = 0.39782 * Sine(fSunTrueLongitude * fDegToRad);
fSunCosDeclination = Cosine(Arcsine(fSunSinDeclination));
//Calc suns local hour angle
fSunLocalHourAngle = (Cosine(fZenith * fDegToRad) - (fSunSinDeclination * Sine(fLatitude * fDegToRad))) / (fSunCosDeclination * Cosine(fLatitude * fDegToRad));
if (fSunLocalHourAngle < -1.0) then //Check if the sun never sets on this date
nSunset=-1; //No sunset
else
fSunH = Arccosine(fSunLocalHourAngle) * fRadToDeg / 15.0;
fLocalMeanTime = fSunH + fSunRAHours - (0.06571 * fApproxTime) - 6.622;
fUTCTime = fLocalMeanTime - fLngHour;
fSunsetTime = fUTCTime + fUTCOffset;
//Adjust to be within 0 and 24
if (fSunsetTime > 24.0) then
fSunsetTime = fSunsetTime - 24.0;
elseif (fSunsetTime < 0.0) then
fSunsetTime = fSunsetTime + 24.0;
endif
nSunset = fSunsetTime * 3600; //Put the result in seconds since midnight;
endif
/////////////Sunrise//////////////
//Calculate an approximate sunrise time for the longitude
fApproxTime = nDayOfYear + ((6.0 - fLngHour)/24.0);
//Sun Mean Anomaly
fSunMeanAnomaly = (.9856 * fApproxTime) - 3.289;
//Sun True Longitude
fSunTrueLongitude = fSunMeanAnomaly + (1.916 * Sine(fSunMeanAnomaly * fDegToRad)) + (.02 * Sine(2 * fSunMeanAnomaly * fDegToRad)) + 282.634;
//Make sure fSunTrueLongitude is between 0 and 360, adjust by 360.
if (fSunTrueLongitude>360.0) then
fSunTrueLongitude = fSunTrueLongitude - 360.0;
elseif (fSunTrueLongitude < 0.0) then
fSunTrueLongitude = fSunTrueLongitude + 360.0;
endif
//Sun Right Ascension (RA)
fSunRightAscension = Arctangent(0.91764 * Tangent(fSunTrueLongitude * fDegToRad)) * fRadToDeg + 360;
//Make sure fSunRightAscension is between 0 and 360, adjust by 360
if (fSunRightAscension>360.0) then
fSunRightAscension = fSunRightAscension - 360.0;
elseif (fSunRightAscension < 0.0) then
fSunRightAscension = fSunRightAscension + 360.0;
endif
//Put RA and L in the same quadrant
fSunTrueLongitudeQuadrant = Truncate(fSunTrueLongitude / 90.0) * 90.0;
fSunRAQuadrant = Truncate(fSunRightAscension / 90.0) * 90.0;
fSunRightAscension = fSunRightAscension + fSunTrueLongitudeQuadrant - fSunRAQuadrant;
//RA converted to Hours
fSunRAHours = fSunRightAscension / 15.0;
//Calc suns Declination
fSunSinDeclination = 0.39782 * Sine(fSunTrueLongitude * fDegToRad);
fSunCosDeclination = Cosine(Arcsine(fSunSinDeclination));
//Calc suns local hour angle
fSunLocalHourAngle = (Cosine(fZenith * fDegToRad) - (fSunSinDeclination * Sine(fLatitude * fDegToRad))) / (fSunCosDeclination * Cosine(fLatitude * fDegToRad));
if (fSunLocalHourAngle > 1.0) then //Check if the sun never rises on this date
nSunrise=-1; //No sunrise
else
fSunH = (360 - Arccosine(fSunLocalHourAngle) * fRadToDeg) / 15.0;
fLocalMeanTime = fSunH + fSunRAHours - (0.06571 * fApproxTime) - 6.622;
fUTCTime = fLocalMeanTime - fLngHour;
fSunriseTime = fUTCTime + fUTCOffset;
//Adjust to be within 0 and 24
if (fSunriseTime > 24.0) then
fSunriseTime = fSunriseTime - 24.0;
elseif (fSunriseTime < 0.0) then
fSunriseTime = fSunriseTime + 24.0;
endif
nSunrise = fSunriseTime * 3600; //Put the result in seconds since midnight;
endif
```