Page 1 of 2
Strange sunset/sunrise calulation
Posted: Wed Dec 15, 2010 1:16 pm
by dali
Hi,
I'm using the "Dark" variable to turn on the outside lights at home.
The calculation of sunset/sunrise seems to be wrong. According to DomotiGa, the sunset today is at 11:47 (which is definitely wrong). According to
timeanddate.com/worldclock/astronomy.ht ... tml?n=1382 the sunset is at 15:15, which is correct.
I'll try to take a look at the calculations, but the function isn't that well documented.
Re: Strange sunset/sunrise calulation
Posted: Wed Dec 15, 2010 1:30 pm
by RDNZL
The sunrise/set code is ported from the links in the source code Astro.module.
It's a known issue and on my todo list, I have tried to look for better code, but haven't find it yet.
Some code takes daylight saving time in account other not, not clear why...
Re: Strange sunset/sunrise calulation
Posted: Wed Dec 15, 2010 10:48 pm
by dali
I've ported this php version (which is based on
williams.best.vwh.net/sunrise_sunset_al ... orithm.htm:
Code: Select all
define("SUNRISE", 0);
define("SUNSET", 1);
class SUN {
var $timezone;
var $latitude;
var $longitude;
var $zenith;
var $date;
var $type;
var $coords;
function sin($d) {
return sin(deg2rad($d));
}
function cos($d) {
return cos(deg2rad($d));
}
function tan($d) {
return tan(deg2rad($d));
}
function atan($x) {
return rad2deg(atan($x));
}
function asin($x) {
return rad2deg(asin($x));
}
function acos($x) {
return rad2deg(acos($x));
}
function set_timezone($timezone) {
$this->timezone = $timezone;
}
function set_latitude($latitude) {
$this->latitude = $this->coords2dec($latitude);
}
function set_longitude($longitude) {
$this->longitude = $this->coords2dec($longitude);
}
function set_zenith($zenith) {
$this->zenith = $zenith;
}
function set_coords($coords) {
if(preg_match("/([^\s]+) ([^\s]+)/", $coords, $matches)) {
$this->set_latitude($this->coords2dec($matches[1]));
$this->set_longitude($this->coords2dec($matches[2]));
$this->coords = $coords;
} else return;
}
function get_timezone() {
return $this->timezone;
}
function get_latitude() {
return $this->latitude;
}
function get_longitude() {
return $this->longitude;
}
function get_zenith() {
return $this->zenith;
}
function get_coords() {
return $this->coords;
}
function coords2dec($coords) {
if(preg_match("/([0-9]+)°([0-9]+)'([0-9\.]+)\"([NnWwEeSs])/", $coords, $matches)) {
$const = (strtolower($matches[4]) == "e") || (strtolower($matches[4]) == "n") ? 1 : -1;
return $const * ($matches[1] + $matches[2] / 60 + $matches[3] / 3600);
} else return $coords;
}
function SUN() {
$this->zenith = 90 + (50 / 60);
$this->timezone = date("O") / 100;
}
function calculation() {
$n = (int)date("z", $this->date);
$lnghour = $this->longitude / 15;
if($this->type == SUNSET) $t = $n + ((18 - $lnghour) / 24);
else $t = $n + ((6 - $lnghour) / 24);
$m = (0.9856 * $t) - 3.289;
$l = $m + (1.916 * $this->sin($m)) + (0.020 * $this->sin(2 * $m)) + 282.634;
$l = fmod($l, 360);
$ra = $this->atan(0.91764 * $this->tan($l));
$ra = fmod($ra, 360);
$lquadrant = floor($l / 90) * 90;
$raquadrant = floor($ra / 90) * 90;
$ra = $ra + ($lquadrant - $raquadrant);
$ra = $ra / 15;
$sindec = 0.39782 * $this->sin($l);
$cosdec = $this->cos($this->asin($sindec));
$cosh = ($this->cos($this->zenith) - ($sindec * $this->sin($this->latitude))) / ($cosdec * $this->cos($this->latitude));
if(($this->type == SUNRISE) && ($cosh > 1)) return NULL;
if(($this->type == SUNSET) && ($cosh < -1)) return NULL;
if($this->type == SUNSET) $h = $this->acos($cosh);
else $h = 360 - $this->acos($cosh);
$h = $h / 15;
$t = $h + $ra - (0.06571 * $t) - 6.622;
$ut = $t - $lnghour;
$ut = fmod($ut, 24);
$summertime = date("I", $this->date + 24 * 60 * 60) ? 1 : 0;
$localt = $ut + $this->timezone + $summertime;
return $this->date + $localt * 60 * 60;
}
function rise($date = 0) {
$this->type = SUNRISE;
$this->date = $date ? $date : mktime(0, 0, 0, date("m"), date("d"), date("Y"));
return date("H:i", $this->calculation());
}
function set($date = 0) {
$this->type = SUNSET;
$this->date = $date ? $date : mktime(0, 0, 0, date("m"), date("d"), date("Y"));
return date("H:i", $this->calculation());
}
function dawn($date = 0) {
$this->type = SUNRISE;
$this->date = $date ? $date : mktime(0, 0, 0, date("m"), date("d"), date("Y"));
return date("H:i", $this->calculation() - 30 * 60);
}
function dusk($date = 0) {
$this->type = SUNSET;
$this->date = $date ? $date : mktime(0, 0, 0, date("m"), date("d"), date("Y"));
return date("H:i", $this->calculation() + 30 * 60);
}
}
What I came up with is this:
Code: Select all
PUBLIC FUNCTION CalcSunTimes(lon AS Float, lat AS Float, timezone AS Float, isRise AS Boolean, twilight AS Integer) AS String
DIM n, lnghour, m, l, ra, lquadrant, raquadrant, sindec, cosdec, coshour, zenith, h, t, ut, localt, ihour, imin AS Float
' twilight setting
IF (twilight = 0) THEN zenith = 90 + (50 / 60) ' effective sunrise/sunset
IF (twilight = 1) THEN zenith = 96 ' civil twilight (brightest)
IF (twilight = 2) THEN zenith = 102 ' nautical twilight
IF (twilight = 3) THEN zenith = 108 ' astronomical twilight (darkest)
n = DateDiff("01/01/" & Year(Now()), Now(), gb.day)
lnghour = lon / 15
IF (isRise) THEN
t = n + ((6 - lnghour) / 24)
ELSE
t = n + ((18 - lnghour) / 24)
ENDIF
m = (0.9856 * t) - 3.289
l = m + (1.916 * Sin(Rad(m))) + (0.020 * Sin(Rad(2 * m))) + 282.634
IF (l < 360) THEN
l += 360
ELSE IF (l > 360) THEN
l -= 360
ENDIF
ra = Deg(ATan(0.91764 * Tan(Rad(l))))
IF (ra < 360) THEN
ra += 360
ELSE IF (ra > 360) THEN
ra -= 360
ENDIF
lquadrant = Fix(l / 90) * 90
raquadrant = Fix(ra / 90) * 90
ra = ra + (lquadrant - raquadrant)
ra = ra / 15
sindec = 0.39782 * Sin(Rad(l))
cosdec = Cos(Rad(Deg(ASin(sindec))))
coshour = (Cos(Rad(zenith)) - (sindec * Sin(Rad(lat)))) / (cosdec * Cos(Rad(lat)))
IF ((isRise) & (coshour > 1)) THEN PRINT "Midnight sun"
IF ((NOT isRise) & (coshour < -1)) THEN PRINT "Polar night"
IF (isRise) THEN
h = 360 - Deg(ACos(coshour))
ELSE
h = Deg(ACos(coshour))
ENDIF
h = h / 15
t = h + ra - (0.06571 * t) - 6.622
ut = t - lnghour
IF (ut < 24) THEN
ut += 24
ELSE IF (ut > 24) THEN
ut -= 24
ENDIF
localt = ut + timezone
ihour = Hour(localt)
imin = Minute(localt)
RETURN Format(ihour, "0#") & ":" & Format(imin, "0#")
It sort of works, but it returns sunset when asking for sunrise, and vice versa.
I can't figure out why! Of course I could just swap the isRise variable, but that would be cheating. I want to know why it's wrong...
Re: Strange sunset/sunrise calulation
Posted: Thu Dec 16, 2010 12:21 am
by dali
The problem was unix time versus gambas time.
It is now returning the same result as
timeanddate.com/worldclock/astronomy.ht ... tml?n=1382
This version is working, I'd be glad if someone else could try it out.
Code: Select all
PUBLIC FUNCTION CalcSunTimes(lon AS Float, lat AS Float, timezone AS Float, isRise AS Boolean, twilight AS Integer) AS String
DIM n, lnghour, m, l, ra, lquadrant, raquadrant, sindec, cosdec, coshour, zenith, h, t, ut, ihour, imin, gTime AS Float
DIM localt AS Integer
' twilight setting
IF (twilight = 0) THEN zenith = 90 + (50 / 60) ' effective sunrise/sunset
IF (twilight = 1) THEN zenith = 96 ' civil twilight (brightest)
IF (twilight = 2) THEN zenith = 102 ' nautical twilight
IF (twilight = 3) THEN zenith = 108 ' astronomical twilight (darkest)
' first calculate the day of the year
n = DateDiff("01/01/" & Year(Now()), Now(), gb.day)
'convert the longitude to hour value and calculate an approximate time
lnghour = lon / 15
IF (isRise) THEN
t = n + ((6 - lnghour) / 24) ' looking for sunrise
ELSE
t = n + ((18 - lnghour) / 24) ' looking for sunset
ENDIF
' calculate the Sun's mean anomaly
m = (0.9856 * t) - 3.289
' calculate the Sun's true longitude
l = m + (1.916 * Sin(Rad(m))) + (0.020 * Sin(Rad(2 * m))) + 282.634
' L potentially needs to be adjusted into the range [0,360] by adding/subtracting 360
IF (l < 360) THEN
l += 360
ELSE IF (l > 360) THEN
l -= 360
ENDIF
' calculate the Sun's right ascension
ra = Deg(ATan(0.91764 * Tan(Rad(l))))
' RA potentially needs to be adjusted into the range [0,360] by adding/subtracting 360
IF (ra < 360) THEN
ra += 360
ELSE IF (ra > 360) THEN
ra -= 360
ENDIF
' right ascension value needs to be in the same quadrant as L
lquadrant = Fix(l / 90) * 90
raquadrant = Fix(ra / 90) * 90
ra = ra + (lquadrant - raquadrant)
' right ascension value needs to be converted into hours
ra = ra / 15
' calculate the Sun's declination
sindec = 0.39782 * Sin(Rad(l))
cosdec = Cos(Rad(Deg(ASin(sindec))))
' calculate the Sun's local hour angle
coshour = (Cos(Rad(zenith)) - (sindec * Sin(Rad(lat)))) / (cosdec * Cos(Rad(lat)))
' this is not working properly, trying to figure out why...
'IF ((isRise) & (coshour > 1)) THEN PRINT "Midnight sun"
'IF ((NOT isRise) & (coshour < -1)) THEN PRINT "Polar night"
' finish calculating H and convert into hours
IF (isRise) THEN
h = 360 - Deg(ACos(coshour)) ' looking for sunrise
ELSE
h = Deg(ACos(coshour)) ' looking for sunset
ENDIF
h = h / 15
' calculate local mean time of rising/setting
t = h + ra - (0.06571 * t) - 6.622
' adjust back to UTC
ut = t - lnghour
' UT potentially needs to be adjusted into the range [0,24] by adding/subtracting 24
IF (ut < 24) THEN
ut += 24
ELSE IF (ut > 24) THEN
ut -= 24
ENDIF
'convert UT value TO local Time zone OF latitude / longitude
'convert hours to seconds
'add seconds from Jan 1st, 1970 until todays date to convert value to unix time (gambas doesnt calculate time in the unix way)
localt = Int((CFloat(Date(Year(Now()), Month(Now()), Day(Now))) - CFloat(Date(1970, 1, 1))) * 86400) + (ut + timezone) * 60 * 60
' add unix time to Jan 1st, 1970 to get gambas time
gTime = DateAdd(Date(1970, 1, 1), gb.Second, localt)
' get hour and minute for sunset/sunrise
ihour = Hour(gTime)
imin = Minute(gTime)
RETURN Format(ihour, "0#") & ":" & Format(imin, "0#")
END
The code isn't good looking, I will comment everything as soon as I've verified it.
EDIT: Added comments to the code.
Re: Strange sunset/sunrise calulation
Posted: Thu Dec 16, 2010 4:57 pm
by spierie
Great you are putting some time into this, I have always been wondering if I was the only one with this problem.....
I will test this weekend to see if it works with me too.
Rgds,
Michel.
Re: Strange sunset/sunrise calulation
Posted: Thu Dec 16, 2010 5:13 pm
by dali
spierie,
Ron has added the code to the svn, just download the latest revision.
Re: Strange sunset/sunrise calulation
Posted: Mon Dec 20, 2010 2:12 pm
by peevee
Code works fine, as far as i can tell.
Nice job dali!
Re: Strange sunset/sunrise calulation
Posted: Mon Dec 20, 2010 2:51 pm
by dali
Thanks, but all I did was translating php code to gambas code.
Re: Strange sunset/sunrise calulation
Posted: Mon Dec 20, 2010 2:54 pm
by spierie
Yep, installed the 509 code yesterday, and it seems to work fine now.
Sunrise 08:48 Sunset 16:28
Thx Dali!
Rgds,
Michel.
Re: Strange sunset/sunrise calulation
Posted: Wed Feb 16, 2011 10:59 pm
by structor
First of all thanks for posting this, very helpfull
.
But some parts of the code don't feel right e.g.:
Code: Select all
IF (ra < 360) THEN
ra += 360
ELSE IF (ra > 360) THEN
ra -= 360
ENDIF
The key idea is to map in the range 0..360. So suppose we feed an ra of 350. A value of 360 will be added.... So i think the first check should be ra<0. Any ideas on this??
Re: Strange sunset/sunrise calulation
Posted: Wed Feb 16, 2011 11:57 pm
by dali
You're right, of course.
It should be like this instead:
Code: Select all
IF (ra < 360) THEN
ra += 360
ELSE IF (ra > 360) THEN
ra -= 360
ENDIF
The same goes for the "l" part:
Code: Select all
IF (l < 0) THEN
l += 360
ELSE IF (l > 360) THEN
l -= 360
ENDIF
You said that "some parts" don't feel right, have you found more errors? Please tell me if you have.
Ron, will you update the source code?
Re: Strange sunset/sunrise calulation
Posted: Thu Feb 17, 2011 10:10 am
by structor
You're correct in my code I also corrected these. At the moment I don't have my code, but the red part also is a bit strange..
cosdec = Cos(Rad(Deg(ASin(sindec))))
I will post my c# version here (for comparison purposes) in a week or so...
What does this code do with daysaving time is that taken into account???
Re: Strange sunset/sunrise calulation
Posted: Mon Feb 21, 2011 2:14 pm
by RDNZL
I have updated code with changes above.
Been a bit busy with non-domotiga stuff lately.
At home with mounting my LCD tv and soundbar to the living room wall, and installing glashart FttH services.
And at work with migrating from one MPLS provider to another (European wide).
Re: Strange sunset/sunrise calulation
Posted: Wed Mar 30, 2011 7:26 pm
by wwolkers
Digging up an old post..
Could it be that the sunset/sunrise calculations don't take summertime into account?
It seems my sunset/sunrise is off by one hour lately.
Re: Strange sunset/sunrise calulation
Posted: Sun Apr 03, 2011 10:30 am
by RDNZL
Yes, daylight saving time is not taken into account for sunset/rise.
Anyone have an idea to add it?
Then we also need to have the dates on which DST is active.