 |
 |
| View previous topic :: View next topic |
| Author |
Message |
carvell
Joined: 20 Sep 2016
Posts: 4
|
| Sunset time calculation |
 Posted: Tue Sep 20, 2016 1:42 pm |
 |
|
There was a sunrise/sunset time calculator on this forum, but it was intertwined with a lot of other functions and was fairly huge, so I implemented the following.
I followed this method, implementing it step by step:
http://williams.best.vwh.net/sunrise_sunset_algorithm.htm
Confirmed it against a variety of sources.
Parameters:
int8 offset - Sunset/sunrise times are calculated in UTC. The offset given here will be applied to the UTC result before returning a value.
int1 riseorset - 1 to calculate rise time, 0 to calculate set time
int16 day - Day of the month, e.g. 31.
int16 month - Month, e.g. 12.
int16 year - Year, e.g. 2016.
float lat - Latitude position, e.g. 51.5074 for London.
float lon - Longitude position, e.g. -0.1278 for London.
The time is returned in decimal format. I.e., for 18:30 it will return 18.5, for 09:45 it'll return 9.75, etc.
| Code: |
#pragma device ANSI //this is required in order that variables default to signed instead of unsigned
#include <math.h>
#define PI 3.141592654
#define ZENITH 90.833333333
float calc_suntime(int8 offset, int1 riseorset, int16 day, int16 month, int16 year, float lat, float lon)
{
//offset: offset to apply from UTC (in hours)
//riseorset: 1 for rise time, 0 for set time
float lngHour, t, M, L, RA, Lquadrant, RAquadrant, sinDec;
float cosDec, cosH, H, T, UT, localT;
int16 N1, N2, N3, N;
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;
lngHour = lon / 15;
if (riseorset)
{
//sunrise time
t = N + ((6 - lngHour) / 24);
}
else
{
//sunset time
t = N + ((18 - lngHour) / 24);
}
M = (0.9856 * t) - 3.289;
L = M + (1.916 * sin((PI/180) * M)) + (0.020 * sin((PI/180) * 2 * M)) + 282.634;
//need to adjust L to be in range (0,360)
while ((L > 360) || (L < 0))
{
if (L > 360) L -= 360;
if (L < 0) L += 360;
}
RA = (180/PI) * atan(0.91764 * tan((PI/180) * L));
//need to adjust RA to be in range (0,360)
while ((RA > 360) || (RA < 0))
{
if (RA > 360) RA -= 360;
if (RA < 0) RA += 360;
}
Lquadrant = (floor( L/90)) * 90;
RAquadrant = (floor(RA/90)) * 90;
RA = RA + (Lquadrant - RAquadrant);
RA = RA / 15;
sinDec = 0.39782 * sin((PI/180) * L);
cosDec = cos((PI/180) * ((180/PI) * asin(sinDec)));
cosH = (cos((PI/180) * ZENITH) - (sinDec * sin((PI/180) * lat))) / (cosDec * cos((PI/180) * lat));
/*
if (cosH > 1)
the sun never rises on this location (on the specified date)
if (cosH < -1)
the sun never sets on this location (on the specified date)
*/
if (riseorset)
{
//sunrise time
H = 360 - ((180/PI) * acos(cosH));
}
else
{
//sunset time
H = (180/PI) * acos(cosH);
}
H = H / 15;
T = H + RA - (0.06571 * t) - 6.622;
UT = T - lngHour;
while ((UT > 24) || (UT < 0))
{
if (UT > 360) UT -= 24;
if (UT < 0) UT += 24;
}
localT = UT + offset;
return localT;
} |
|
|
 |
Ttelmah
Joined: 11 Mar 2010
Posts: 19510
|
|
| Posted: Mon Oct 10, 2016 1:44 am |
|
|
One comment.
You need to think in terms of using a smaller unit for 'offset'. Quite a lot of countries/places use 0.5 hour timezones. (Adelaide in Australia, India, and Newfoundland for example).
So I'd suggest changing 'offset' to be in integer 0.5hour steps, then the offset line just becomes:
localT = UT + (offset*0.5);
Which can then handle these locations. 
Best Wishes |
|
|
|
|
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|
FAQ 
Search
Register
Profile 
Log in to check your private messages 
Log in
