Strange sunset/sunrise calulation

[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.

[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...

[dali]

I've ported this php version (which is based onwilliams.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...

[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.

[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.

[dali]

spierie,

Ron has added the code to the svn, just download the latest revision.

[peevee]

Code works fine, as far as i can tell.
Nice job dali!

[dali]

Thanks, but all I did was translating php code to gambas code.

[spierie]

Yep, installed the 509 code yesterday, and it seems to work fine now.
Sunrise 08:48 Sunset 16:28

Thx Dali!

Rgds,
Michel.

[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??

[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?

[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???

[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).

[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.

[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.