This is an empty page.
Enter page contents here, upload content using the button above, or click cancel to leave it empty.
#! /bin/bash # \ exec tclsh "$0" "$@"
#//# # <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