This is a Tcl implementation of sunrise/sunset calculation according 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
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 observers 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.
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]
#//#
# <h2>Package Noon</h2>
# This is a Tcl implementation of sunrise/sunset calculation according 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 observers 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> (90°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° available with <var>$Noon::CIVILIAN</var></dd>
# <dt>nautical Sunrise/Sunset</dt>
# <dd>Since on sea the horizon is earlier and longer visible there is another defined zenith angle for this situation.
# 102° / <var>$Noon::NAUTICAL</var></dd>
# <dt>astronomical Sunrise/Sunset</dt>
# <dd>When the sun is 108° 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