**Calculation of Sunrise and Sunset**

This is a Tcl implementation of sunrise/sunset calculation according to the receipe

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

which is derived from:
 Almanac for computers, 1990
 published by Nautical Almanac Office
 US Naval Observatory
 Washington, DC 20392

[AM] (3 august 2017) The above link fails for me. A Google search didn't bring up a useful alternative. Maybe this is temporary, but otherwise a working alternative would be useful.
Hm, this one perhaps: http://www.edwilliams.org/sunrise_sunset_algorithm.htm



***Sunrise and Sunset Times***

There are 4 different categories of sunrise and sunset times and appropriate twilight intervals.
They differ in the threshold value that denotes how far the solar disk is below the horizon.
The threshold value can be expressed by the zenith angle.
From the position of an observer the zenith can be seen as a pointer from the ground up in the sky
and an arrow can be drawn from the observer's location to the horizon forming this angle.
It is quite clear that this definition depends solely on the coordinates of the location of the observer
and does not consider true horizon.


   '''official Sunrise/Sunset''':   The official Sunrise/Sunset denotes the time when the edge of the solar disk seems to touch the horizon. The zenith value for this kind of sunrise/sunset is available with `$Noon::OFFICIAL` (i.e. 90° 50'). This is also the default value for the procedures in this package.
  
   '''civilian Sunrise/Sunset''':   The civilian Sunrise/Sunset depicts the situation where the horizon is just visible (i.e. the sun is clearly below the horizon but lights up the atmosphere). The defined zenith value is 96° available with `$Noon::CIVILIAN`

   '''nautical Sunrise/Sunset''':   Because on sea the horizon is visible earlier and longer another zenith angle is defined for this situation. 102° / `$Noon::NAUTICAL`


   '''astronomical Sunrise/Sunset''':   When the sun is 108° below the horizon the stars become visible. This zenith angle is defined by `$Noon::ASTRONOMICAL`


Example:

======
set lat   53.918085
set long  12.27929

set y  2016
set m  9
set d  26

set sr  [Noon::sunriseUTC $y $m $d $lat  $long]

======
(Note: east, north are positive values) 
This prints out the time of sunrise (official sunrise) in UTC:  `05:05`

As of 09/10/2017 there seems to be an issue with the code below. This doesn't seem to work:
======
set lat   53.918085
set long  12.27929

set y  2016
set m  9
set d  08

set sr  [Noon::sunriseUTC $y $m $d $lat  $long]

======

Setting day to 08 or 09 doesn't seem to work and gives the error
 missing operator at _@_
 in expression "int(30.0 - (0.0*2.0) + 0_@_8 - 30)";

The current solution I have is to check for preceding 0's.:

======
set systemTime [clock seconds]
    set lat   51.522416
    set long  0

    set y  [clock format $systemTime -format %Y]
    set m  [clock format $systemTime -format %m]
    set d  [clock format $systemTime -format %d]
    
    #Silly loop to see if the first character of the date is 0 - Because
    # for some reason 08 and 09 fails.
    if { [string first "0" $d] == 0 } {
      set dstrip [string trimleft $d 0]
      puts $dstrip
      set sr  [Noon::sunriseUTC $y $m $dstrip $lat $long]
      set ss  [Noon::sunsetUTC $y $m $dstrip $lat $long]
    } else {
      set sr  [Noon::sunriseUTC $y $m $d $lat $long]
      set ss  [Noon::sunsetUTC $y $m $d $lat $long]
    }

======

(2017-10-10): It doesn't work because 08 and 09 are interpreted as (invalid) octal numbers (more or less the same thing happens in Javascript, Perl, Ruby, Python etc). Tcl 9 is supposed to drop the "leading 0 means octal" convention. See [Tcl and octal numbers].

As a slightly less messy workaround, one could use `clock format $systemTime -format %e` to get a day number without leading 0 (it is returned as blank + 1..9 and 10..31, which `expr` can handle).


----
======
#//#
# <h2>Package Noon</h2>
# This is a Tcl implementation of sunrise/sunset calculation according to the <b>receipe</b><br>
# <a id="RECEIPE"></a>
# <blockquote>
# {@link http://williams.best.vwh.net/sunrise_sunset_algorithm.htm http://williams.best.vwh.net/sunrise_sunset_algorithm.htm}
# </blockquote>
# <p>
# which is derived from:<br>
# <blockquote>
# <em>
# Almanac for computers, 1990<br>
# published by Nautical Almanac Office<br>
# US Naval Observatory<br>
# Washington, DC 20392<br>
# </em>
# </blockquote>
# </p>
#
# <h2>Sunrise and Sunset Times</h2>
# There are 4 different categories of sunrise and sunset times and appropriate twilight intervals.
# They differ in the threshold value that denotes how far the solar disk is below the horizon.
# The threshold value can be expressed by the zenith angle.
# From the position of an observer the zenith can be seen as a pointer from the ground up in the sky
# and an arrow can be drawn from the observer's location to the horizon forming this angle.
# It is quite clear that this definition depends solely on the coordinates of the location of the observer
# and does not consider true horizon.
# 
# <dl>
# <dt>official Sunrise/Sunset</dt>
#   <dd>The official Sunrise/Sunset denotes the time when the edge of the solar disk seems to touch the horizon.
#       The zenith value for this kind of sunrise/sunset is available with <var>$Noon::OFFICIAL</var> (i.e. 90&deg;50').
#       This is also the default value for the procedures in this package</dd>
#   <dt>civilian Sunrise/Sunset</dt>
#   <dd>The civilian Sunrise/Sunset depicts the situation where the horizon is just visible (i.e. the sun is clearly below the horizon
#       but lights up the atmosphere).
#       The defined zenith value is 96&deg; available with <var>$Noon::CIVILIAN</var></dd>
#   <dt>nautical Sunrise/Sunset</dt>
#   <dd>Because on sea the horizon is visible earlier and longer another zenith angle is defined for this situation.
#       102&deg; / <var>$Noon::NAUTICAL</var></dd>
#   <dt>astronomical Sunrise/Sunset</dt>
#   <dd>When the sun is 108&deg; below the horizon the stars become visible. This zenith angle is defined by 
#       <var>$Noon::ASTRONOMICAL</var></dd>
# </dl>
# @author T. Heinlein
# @version 1.0
#//#

package provide noon 1.0;

namespace eval Noon  {

  namespace export sunriseUTC;
  namespace export sunriseDecUTC;
  namespace export sunsetUTC;
  namespace export sunsetDecUTC;
  namespace export formatDecTime;

  # zenith controls kind of sunrise/sunset times
  variable OFFICIAL       90.8333;        # solar disk touches horizon
  variable CIVILIAN       96.0;                # horizon is still/already visible
  variable NAUTICAL      102.0;                # horizon is visible on sea
  variable ASTRONOMICAL  108.0;                # the stars become visible / will disappear

  variable D2R  0.01745329252;                # factor degree to rad i.e. PI / 180


  # formats a time given as decimal float value (0 ... 24) in HH:MM presentation
  # Regardless of the timezone offset given the resulting time is garanteed in 00 ... 24 interval.
  # @param dec  time of day as decimal number   e.g. 23.5
  # @param offset  time offset from GMT in hours (decimal) default is 0.
  # @return time in Hours/Minutes format HH:MM  e.g. 23:30
  # {@link #RECEIPE  receipe}
  proc formatDecTime { dec {offset 0} }  {
    set decTime [expr $dec + $offset]

    if { $decTime >= 24.0 }  {
      set $decTime [expr $decTime - 24]
    }
    if { $decTime < 0.0 }  {
      set $decTime [expr $decTime + 24]
    }
    set hours [expr floor($decTime)]
    set frac  [expr $decTime - $hours]
    set mins  [expr round($frac * 60.0)]
    return [format "%02d:%02d" [expr int($hours)] [expr int($mins)]]
  }



  # sinus with argument in degrees
  # @param deg  angle in degrees 0 ... 360
  # @return sinus value -1 ... +1
  proc _sin { deg }  {
    variable D2R;
    return [expr sin($deg * $D2R)]
  }

  # cosinus with argument in degrees
  # @param deg  angle in degrees 0 ... 360
  # @return cosinus value -1 ... +1
  proc _cos { deg }  {
    variable D2R;
    return [expr cos($deg * $D2R)]
  }

  # tangens with argument in degrees
  # @param deg  angle in degrees 0 ... 360
  # @return tangens value 
  proc _tan { deg }  {
    variable D2R;
    return [expr tan($deg * $D2R)]
  }

  # invers tangens returning degrees
  # @param val  value for looking up tangens 
  # @return tangens value in degrees 0 ... 360
  proc _atan { val }  {
    variable D2R;
    return [expr atan($val) / $D2R]
  }

  # invers sinus returning degrees
  # @param val  value for looking up sinus 
  # @return sinus value in degrees 0 ... 360
  proc _asin { val }  {
    variable D2R;
    return [expr asin($val) / $D2R]
  }

  # invers cosinus returning degrees
  # @param val  value for looking up cosine  -1 ... +1
  # @return cosinus value in degrees 0 ... 360
  proc _acos { val }  {
    variable D2R;
    return [expr acos($val) / $D2R]
  }


  # given a date calculates day of the year 1 ... 365 or 366 
  # step 1 from (**)
  # @param y  year (4 digits)
  # @param m  month where January is 1 .... December is 12
  # @param d  day of month  (1 ... 31)
  # @return day of the year (count from 1st January)
  proc _dayOfYear { y m d }  {
    set n1 [expr floor(275 * $m / 9.0) ];
    set n2 [expr floor(($m + 9) / 12.0)];
    set n3 [expr (1 + floor(($y -4 * floor($y / 4.0) + 2) / 3.0))];
    set N  [expr int($n1 - ($n2*$n3) + $d - 30)];
    return $N
  }


  # calculate a time base for sunrise
  # step 2 from (**)
  # @param y  year (4 digits)
  # @param m  month where January is 1 .... December is 12
  # @param d  day of month  (1 ... 31)
  # @param longitude  longitude
  # @return base time of day for calculating sunrise (decimal)
  proc _sunriseBaseTime { y m d longitude }  {
    set N [_dayOfYear $y $m $d]
    set longH  [expr $longitude / 15.0];        # from step 2: Longitude in hours (24h -> 360 deg)
    set bt [expr $N + ((6 - $longH) / 24.0)]
    return $bt
  }

  # calculate a time base for sunset
  # step 2 from (**)
  # @param y  year (4 digits)
  # @param m  month where January is 1 .... December is 12
  # @param d  day of month  (1 ... 31)
  # @param longitude  longitude
  # @return base time of day for calculating sunset (decimal)
  proc _sunsetBaseTime { y m d longitude }  {
    set N [_dayOfYear $y $m $d]
    set longH  [expr $longitude / 15.0];        # from step 2: Longitude in hours (24h -> 360 deg)
    set bt [expr $N + ((18 - $longH) / 24.0)]
    return $bt
  }



  # calculate sun's true longitude
  # step 5a from (**)
  # @param timeBase a time base for sunrise or sunset
  #        (@see #_sunriseBaseTime _sunriseBaseTime
  #         @see #_sunsetBaseTime _sunsetBaseTime)
  # @return the true longitude of the sun
  proc _trueLongitude { timeBase }  {
    set M [expr (0.9856 * $timeBase) - 3.289]; # anomaly of the sun
    set trueL [expr $M + (1.916 * [_sin $M]) + (0.020 * [_sin [expr 2*$M]]) + 282.634]
  
    # align
    if { $trueL > 360.0 }  {
      set trueL [expr $trueL - 360.0]
    } elseif { $trueL < 0.0 }  {
      set trueL [expr $trueL + 360.0]
    }
    return $trueL;
  }




  # compute sun's right ascension
  # steps 5b and 5c from (**)
  # @param trueLong the true longitude of the sun
  # @return the right ascension of the sun in hours
  proc _rightAscension { trueLong }  {
    set ra [ expr [_atan [expr 0.91764 * [_tan $trueLong]]] ]
    # align
    if { $ra > 360.0 }  {
      set ra [expr $ra - 360.0]
    } elseif { $ra < 0.0 }  {
      set ra [expr $ra + 360.0]
    }
    set Lq  [expr (floor($trueLong/90.0) * 90.0)]
    set RAq [expr (floor($ra/90.0) * 90.0)]
    set RA  [expr $ra + ($Lq - $RAq)]
    return [expr $RA / 15.0]
  }


  # compute hour angle for sunrise
  # step 7a from (**)
  # @param lat the latitude
  # @param trueLong_sr the true longitude of the sun at sunrise
  # @param zenith the zenith in degrees
  # @return the hour angle of the sun
  proc _hourAngleSunrise { lat trueLong_sr zenith }  {
    # step 6: 
    set sinDec [expr 0.39782 * [_sin $trueLong_sr]]
    set cosDec [expr [_cos [_asin $sinDec]]]
  
    # step 7:
    set cosH [ expr ([_cos $zenith] - ($sinDec * [_sin $lat])) / ($cosDec * [_cos $lat]) ]
  
    if { $cosH > 1 }  {
      # sun does not rise this day at this location
      return "none"
    }
    set Hrise [expr (360 - [_acos $cosH]) / 15.0]
    return $Hrise
  }


  # compute hour angle for sunset
  # step 7a from (**)
  # @param lat the latitude
  # @param trueLong_ss the true longitude of the sun at sunset
  # @param zenith the zenith in degrees
  # @return the hour angle of the sun
  proc _hourAngleSunset { lat trueLong_ss  zenith }  {
    # step 6: 
    set sinDec [expr 0.39782 * [_sin $trueLong_ss]]
    set cosDec [expr [_cos [_asin $sinDec]]]

    # step 7:
    set cosH [ expr ([_cos $zenith] - ($sinDec * [_sin $lat])) / ($cosDec * [_cos $lat]) ]

    if { $cosH < -1 }  {
      # sun does not set this day at this location
      return "none"
    }
    set Hset [expr [_acos $cosH] / 15.0]
    return $Hset
  }


  # calculates the time of sunrise. Date and Time are UTC.
  # @param y  year (4 digits)
  # @param m  month where January is 1 .... December is 12
  # @param d  day of month  (1 ... 31)
  # @param lat     latitude
  # @param long    longitude
  # @param zenith  the zenith in degrees (default is Noon::OFFICIAL)
  # @return the UTC time of sunrise as a decimal value 0 ... 24
  #    or the string 'none' if there is no sunrise at this location
  proc sunriseDecUTC { y m d lat long {zenith $Noon::OFFICIAL} }  {
    set longH     [expr $long / 15.0];                        # from step 2: Longitude in hours (24h -> 360 deg)
    set baseTime  [_sunriseBaseTime  $y $m $d $long]
    set trueLong  [_trueLongitude    $baseTime]
    set ra        [_rightAscension   $trueLong]
    set ha        [_hourAngleSunrise $lat $trueLong $zenith]

    if { $ha == "none" }  {
      return $ha;
    }

    set t         [expr $ha + $ra - (0.06571 * $baseTime) - 6.622]
    set utc       [expr $t - $longH]

    if { $utc > 24 }  {
      set utc [expr $utc - 24]
    }  elseif { $utc < 0 }  {
      set utc [expr $utc + 24]
    }
    return $utc
  }


  # calculates the time of sunrise. Date and Time are UTC.
  # @param y  year (4 digits)
  # @param m  month where January is 1 .... December is 12
  # @param d  day of month  (1 ... 31)
  # @param lat     latitude
  # @param long    longitude
  # @param zenith  the zenith in degrees (default is Noon::OFFICIAL)
  # @return the UTC time of sunrise in hours and minutes / HH:MM
  #    or the string 'none' if there is no sunrise at this location
  proc sunriseUTC { y m d lat long {zenith $Noon::OFFICIAL} }  {
    set utc [sunriseDecUTC $y $m $d $lat $long $zenith]
    if { $utc == "none" }  {
      return $utc;
    }
    return [formatDecTime $utc]
  }



  # calculates the time of sunset
  # @param y  year (4 digits)
  # @param m  month where January is 1 .... December is 12
  # @param d  day of month  (1 ... 31)
  # @param lat     latitude
  # @param long    longitude
  # @param zenith  the zenith in degrees (default is Noon::OFFICIAL)
  # @return the UTC time of sunrise as a decimal value 0 ... 24
  #    or the string 'none' if there is no sunset at this location
  proc sunsetDecUTC { y m d lat long {zenith $Noon::OFFICIAL} }  {
    set longH     [expr $long / 15.0];                        # from step 2: Longitude in hours (24h -> 360 deg)
    set baseTime  [_sunsetBaseTime  $y $m $d $long]
    set trueLong  [_trueLongitude    $baseTime]
    set ra        [_rightAscension   $trueLong]
    set ha        [_hourAngleSunset  $lat $trueLong $zenith]

    if { $ha == "none" }  {
      return $ha;
    }

    set t         [expr $ha + $ra - (0.06571 * $baseTime) - 6.622]
    set utc       [expr $t - $longH]

    if { $utc > 24 }  {
      set utc [expr $utc - 24]
    }  elseif { $utc < 0 }  {
      set utc [expr $utc + 24]
    }
    return $utc
  }


  # calculates the time of sunset. Date and Time are UTC.
  # @param y  year (4 digits)
  # @param m  month where January is 1 .... December is 12
  # @param d  day of month  (1 ... 31)
  # @param lat     latitude
  # @param long    longitude
  # @param zenith  the zenith in degrees (default is Noon::OFFICIAL)
  # @return the UTC time of sunset in hours and minutes / HH:MM
  #    or the string 'none' if there is no sunset at this location
  proc sunsetUTC { y m d lat long {zenith $Noon::OFFICIAL} }  {
    set utc [sunsetDecUTC $y $m $d $lat $long $zenith]
    if { $utc == "none" }  {
      return $utc;
    }
    return [formatDecTime $utc]
  }


}
# --- end namespace
======

<<categories>>Enter Category Here