Solar Position UDF

tim292stro · March 21, 2011

Hi,

I've been building up an AutoIt home automation script over the past months and one script I just finished this weekend I thought I'd share. This function allows you to calculate your sun-rise, sun-set, and solar-noon times for a given latitude and longitude - plus it allows you to calculate the current relative position of the sun to your position.

The original code was JavaScrip from a NOAA website, they have graciously allowed me to convert it to AutoIt Script, so long as I reference the website that I obtained the code from which is: http://www.esrl.noaa.gov/gmd/grad/solcalc/

I will make the same disclaimer that NOAA makes regarding accuracy: “this calculator function is for ‘entertainment only’ – the author’s cannot be held accountable for accuracy issues and bugs in critical systems”

Here is the UDF (select all and save as "SunPosition.au3" in your include folder...):

#include-once
#include <array.au3>
#include <Math.au3>
#include <Date.au3>


; #INDEX# =======================================================================================================================
; Title .........: SunPosition
; AutoIt Version : 3.3.6++
; Language ......: English
; Description ...: Function that allows the sun's postion to be calculated relative to a geographical position (lat/long)
;
; Author(s) .....: Tim Strommen (tim292stro), original Javascript source from NOAA "http://www.esrl.noaa.gov/gmd/grad/solcalc/"
; Dll(s) ........: N/A
; ===============================================================================================================================

; #CONSTANTS# ===================================================================================================================
Global Const $Pi = 4 * ATan(1)
Global Const $MonthLength[12] = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
Global Const $MonthFullName[12] = ["January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December"]
Global Const $MonthAbrev[12] = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]
; ===============================================================================================================================

; #CURRENT# =====================================================================================================================
;_GetSolarAstronimicalData
; ===============================================================================================================================

; #FUNCTION# ====================================================================================================================
; Name...........: _GetSolarAstronimicalData
; Description ...: Calculate the position of the sun relative to a position on earth for a given time/date.
; Syntax.........: _GetSolarAstronimicalData ( $Latitude, $Longitude, $Month, $MDay, $Year, $TimeZone, $Hour, $Minute, $Second, $DST )
; Parameters ....: $Latitude     - The Latitude portion of the earth position to calculate the sun's position relative to (float)
;                  $Longitude     - The Latitude portion of the earth position to calculate the sun's position relative to (float)
;                  $Month     - The month of the year you are calculating for (numeric-string: "1-12")
;                  $MDay     - The day of the month you are calculating for (numeric-string: "1-31")
;                  $Year     - The year you are calculating for (numeric-string: "-2000 to 3000")
;                  $TimeZone     - The time zone you are calculating for (numeric-string: "-11 to 12")
;                  $Hour     - Optional!  The hour you are calculating for in 24-hrs (numeric-string: "0-23")
;                  $Minute     - Optional!  The minute you are calculating for (numeric-string: "00-59")
;                  $Second    - Optional!  The second you are calculating for (numeric-string: "00-59")
;                  $DST     - Optional!  Is Daylight Saving's Time in effect? (boolean)
; Return values .: Success - Returns an array with the following values:
;                            0 = Element count. If count = 3, then only items 1-3 available.  If count = 6, then all elements available)
;                            1 = Sunrise time "07:12" 24-hr time, or "07:12 Jul 17" if position is near international date line
;                            2 = Solar Noon "13:16:46" 24-hr time with seconds.  Higest point of sun in sky.
;                            3 = Sunset time "19:22" 24-hr time, or "19:22 Jul 19" if position is near international date line
;                            4 = Sun Azimuth in degrees relative to computed position (0deg = true north, 180deg = true south, 90deg = East, 270 = West)
;                            5 = Sun Elevation is degrees relative to computed position at sea-level
;                            6 = Ambient light disposition (possibilities: Day, Civil Twilight, Nautical Twilight, Astronomical Twilight, Night)
; Author ........: Tim Strommen (tim292stro)
; Modified.......:
; Remarks .......: Again, special thanks to NOAA for allowing re-use of their code!!  See: "http://www.esrl.noaa.gov/"
; Related .......:
; Link ..........:
; Example .......: Yes
; ===============================================================================================================================
Func _GetSolarAstronimicalData ( $AstroLat = "", $AstroLong = "", $AstroMonth = "", $AstroDay = "", $AstroYear = "", $AstroTimeZone = "", $AstroHour = "", $AstroMinute = "", $AstroSecond = "", $AstroDST = False )
    Local $AstroBeginAT, $AstroBeginNT, $AstroBeginCT, $AstroSunrise, $AstroSolarNoon, $AstroSunset, $AstroEndCT, $AstroEndNT, $AstroEndAT, $AstroSolarElevation, $AstroSolarAzimuth
    If ( $AstroLat = "" ) OR ( $AstroLong = "" ) Then
        ; Return current Greenwich, UK basic times with current solar position
        $AstroLat = "51.48"
        $AstroLong = "0.000000001"
        $AstroTimeZone = "0"
        $AstroMonth = @MON
        $AstroDay = @MDAY
        $AstroYear = @YEAR
        $AstroHour = @HOUR
        $AstroMinute = @MIN
        $AstroSecond = @SEC
        $RetunedArray = _SolarCalculate ( $AstroLat, $AstroLong, $AstroMonth, $AstroDay, $AstroYear, $AstroTimeZone, $AstroHour, $AstroMinute, $AstroSecond, $AstroDST, True )
        Return $RetunedArray
    ElseIf $AstroHour = "" Then
        ; Just return the specified day's basic times, no current solar position
        $AstroHour = "12"
        $AstroMinute = "00"
        $AstroSecond = "00"
        $RetunedArray = _SolarCalculate ( $AstroLat, $AstroLong, $AstroMonth, $AstroDay, $AstroYear, $AstroTimeZone, $AstroHour, $AstroMinute, $AstroSecond, $AstroDST, False )
        Return $RetunedArray
    Else
        ; Return both the basic times, plus the current solar position
        $RetunedArray = _SolarCalculate ( $AstroLat, $AstroLong, $AstroMonth, $AstroDay, $AstroYear, $AstroTimeZone, $AstroHour, $AstroMinute, $AstroSecond, $AstroDST, True )
        Return $RetunedArray
    EndIf
EndFunc ;==>_GetSolarAstronimicalData




#cs
********** Developer Notice!!! **********
The following are the original NOAA
Javascripts with minor modifications to
support programmatic query instead of
web-form query.  The original Javascript
text preceeds the AutoIt-v3 conversion
and are provided with gracious permission
from NOAA.

The original webform and calculator is located at: "http://www.esrl.noaa.gov/gmd/grad/solcalc/"

Thank you again NOAA!!

-Tim S.



<script type="text/javascript">

function calcTimeJulianCent(jd)
{
  var T = (jd - 2451545.0)/36525.0
  return T
}
#ce
Func _CalcTimeJulianCent( $jd = 0 )
    Local $Time = ( $jd - 2451545.0 ) / 36525.0
    Return $Time
EndFunc
#cs




function calcJDFromJulianCent(t)
{
  var JD = t * 36525.0 + 2451545.0
  return JD
}
#ce
Func _CalcJDFromJulianCent ( $t = 0 )
    Local $JD = $t * 36525.0 + 2451545.0
    Return $JD
EndFunc
#cs




function isLeapYear(yr)
{
  return ((yr % 4 == 0 && yr % 100 != 0) || yr % 400 == 0);
}
#ce
Func _isLeapYear ( $yr = 0 )
    Return ( ( Mod ( $yr, 4 ) And Mod ($yr, 100 ) ) Or ( Mod ( $yr, 400 ) = 0 ) )
EndFunc
#cs




function calcDoyFromJD(jd)
{
  var z = Math.floor(jd + 0.5);
  var f = (jd + 0.5) - z;
  if (z < 2299161) {
    var A = z;
  } else {
    alpha = Math.floor((z - 1867216.25)/36524.25);
    var A = z + 1 + alpha - Math.floor(alpha/4);
  }
  var B = A + 1524;
  var C = Math.floor((B - 122.1)/365.25);
  var D = Math.floor(365.25 * C);
  var E = Math.floor((B - D)/30.6001);
  var day = B - D - Math.floor(30.6001 * E) + f;
  var month = (E < 14) ? E - 1 : E - 13;
  var year = (month > 2) ? C - 4716 : C - 4715;

  var k = (isLeapYear(year) ? 1 : 2);
  var doy = Math.floor((275 * month)/9) - k * Math.floor((month + 9)/12) + day -30;
  return doy;
}
#ce
Func _calcDoyFromJD ( $jd = 0 )
    Local $z = Floor ( $jd + 0.5 )
    Local $f = ( $jd + 0.5 ) - $z
    Local $A, $alpha, $B, $C, $D, $doy, $E, $day, $month, $year, $K
    If $z < 2299161  Then
        $A = $z
    Else
        $alpha = Floor ( ( $z - 1867216.25 ) / 36524.25 )
        $A = $z + 1 + $alpha - Floor ( $alpha / 4 )
    EndIf
    $B = $A + 1524
    $C = Floor ( ( $B - 122.1 ) / 365.25 )
    $D = Floor ( 365.25 * $C )
    $E = Floor ( ( $B - $D ) / 30.6001 )
    $day = $B - $D - Floor ( 30.6001 * $E ) + $f
    If $E < 14 Then
        $month = $E - 1
    Else
        $month = $E - 13
    EndIf
    If $month > 2 Then
        $year = $C - 4716
    Else
        $year = $C - 4715
    EndIf
    If _isLeapYear ( $year ) Then
        $K = 1
    Else
        $K = 2
    EndIf
    $doy = Floor ( ( 275 * $month ) / 9 ) - $K * Floor ( ( $month + 9 ) / 12 ) + $day - 30
    Return $doy
EndFunc
#cs




function radToDeg(angleRad)
{
  return (180.0 * angleRad / Math.PI);
}
#ce
Func _radToDeg ( $AngleRad = 0 )
    Return 180.0 * $AngleRad / $Pi
EndFunc
#cs




function degToRad(angleDeg)
{
  return (Math.PI * angleDeg / 180.0);
}
#ce
Func _degToRad ( $AngleDeg = 0 )
    Return $Pi * $AngleDeg / 180.0
EndFunc
#cs




function calcGeomMeanLongSun(t)
{
  var L0 = 280.46646 + t * (36000.76983 + t*(0.0003032))
  while(L0 > 360.0)
  {
    L0 -= 360.0
  }
  while(L0 < 0.0)
  {
    L0 += 360.0
  }
  return L0     // in degrees
}
#ce
Func _calcGeomMeanLongSun ( $t = 0 )
    Local $L0 = 280.46646 + $t * ( 36000.76983 + $t * (0.0003032) )
    While $L0 > 360.0
        $L0 -= 360.0
    WEnd
    While $L0 < 0.0
        $L0 += 360.0
    WEnd
    Return $L0 ; in degrees
EndFunc
#cs




function calcGeomMeanAnomalySun(t)
{
  var M = 357.52911 + t * (35999.05029 - 0.0001537 * t);
  return M;     // in degrees
}
#ce
Func _calcGeomMeanAnomalySun ( $t = 0 )
    Local $M = 357.52911 + $t * (35999.05029 - 0.0001537 * $t )
    Return $M  ; in degrees
EndFunc
#cs




function calcEccentricityEarthOrbit(t)
{
  var e = 0.016708634 - t * (0.000042037 + 0.0000001267 * t);
  return e;     // unitless
}
#ce
Func _calcEccentricityEarthOrbit ( $t = 0 )
    Local $e = 0.016708634 - $t * ( 0.000042037 + 0.0000001267 * $t )
    Return $e  ; unitless
EndFunc
#cs




function calcSunEqOfCenter(t)
{
  var m = calcGeomMeanAnomalySun(t);
  var mrad = degToRad(m);
  var sinm = Math.sin(mrad);
  var sin2m = Math.sin(mrad+mrad);
  var sin3m = Math.sin(mrad+mrad+mrad);
  var C = sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) + sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289;
  return C;     // in degrees
}
#ce
Func _calcSunEqOfCenter ( $t = 0 )
    Local $m = _calcGeomMeanAnomalySun ( $t )
    Local $mrad = _degToRad ( $m )
    Local $sinm = Sin ( $mrad )
    Local $sin2m = Sin ( $mrad + $mrad )
    Local $sin3m = Sin ( $mrad + $mrad + $mrad )
    Local $C = $sinm * ( 1.914602 - $t * ( 0.004817 + 0.000014 * $t ) ) + $sin2m * ( 0.019993 - 0.000101 * $t ) + $sin3m * 0.000289
    Return $C
EndFunc
#cs




function calcSunTrueLong(t)
{
  var l0 = calcGeomMeanLongSun(t);
  var c = calcSunEqOfCenter(t);
  var O = l0 + c;
  return O;     // in degrees
}
#ce
Func _calcSunTrueLong ( $t = 0 )
    Local $lo = _calcGeomMeanLongSun ( $t )
    Local $c = _calcSunEqOfCenter ( $t )
    Local $O = $lo + $c
    Return $O  ; in degrees
EndFunc
#cs




function calcSunTrueAnomaly(t)
{
  var m = calcGeomMeanAnomalySun(t);
  var c = calcSunEqOfCenter(t);
  var v = m + c;
  return v;     // in degrees
}
#ce
Func _calcSunTrueAnomaly ( $t = 0 )
    Local $m = _calcGeomMeanAnomalySun ( $t )
    Local $c = _calcSunEqOfCenter ( $t )
    Local $v = $m + $c
    Return $v
EndFunc
#cs




function calcSunRadVector(t)
{
  var v = calcSunTrueAnomaly(t);
  var e = calcEccentricityEarthOrbit(t);
  var R = (1.000001018 * (1 - e * e)) / (1 + e * Math.cos(degToRad(v)));
  return R;     // in AUs
}
#ce
Func _calcSunRadVector ( $t = 0 )
    Local $v = _calcSunTrueAnomaly ( $t )
    Local $e = _calcEccentricityEarthOrbit ( $t )
    Local $R = ( 1.000001018 * ( 1 - $e * $e ) ) / ( 1 + $e * cos ( _degToRad ( $v ) ) )
    Return $R  ; in AUs
EndFunc
#cs




function calcSunApparentLong(t)
{
  var o = calcSunTrueLong(t);
  var omega = 125.04 - 1934.136 * t;
  var lambda = o - 0.00569 - 0.00478 * Math.sin(degToRad(omega));
  return lambda;        // in degrees
}
#ce
Func _calcSunApparentLong ( $t = 0 )
    Local $o = _calcSunTrueLong ( $t )
    Local $omega = 125.04 - 1934.136 * $t
    Local $lambda = $o - 0.00569 - 0.00478 * Sin ( _degToRad ( $omega ) )
    Return $lambda  ; in degrees
EndFunc
#cs




function calcMeanObliquityOfEcliptic(t)
{
  var seconds = 21.448 - t*(46.8150 + t*(0.00059 - t*(0.001813)));
  var e0 = 23.0 + (26.0 + (seconds/60.0))/60.0;
  return e0;        // in degrees
}
#ce
Func _calcMeanObliquityOfEcliptic ( $t = 0 )
    Local $seconds = 21.448 - $t * ( 46.8150 + $t * ( 0.00059 - $t * ( 0.001813 ) ) )
    Local $e0 = 23.0 + ( 26.0 + ( $seconds / 60.0 ) ) / 60.0
    Return $e0  ; in degrees
EndFunc
#cs




function calcObliquityCorrection(t)
{
  var e0 = calcMeanObliquityOfEcliptic(t);
  var omega = 125.04 - 1934.136 * t;
  var e = e0 + 0.00256 * Math.cos(degToRad(omega));
  return e;     // in degrees
}
#ce
Func _calcObliquityCorrection ( $t = 0 )
    Local $e0 = _calcMeanObliquityOfEcliptic ( $t )
    Local $omega = 125.04 - 1934.136 * $t
    Local $e = $e0 + 0.00256 * Cos ( _degToRad ( $omega ) )
    Return $e
EndFunc
#cs




function calcSunRtAscension(t)
{
  var e = calcObliquityCorrection(t);
  var lambda = calcSunApparentLong(t);
  var tananum = (Math.cos(degToRad(e)) * Math.sin(degToRad(lambda)));
  var tanadenom = (Math.cos(degToRad(lambda)));
  var alpha = radToDeg(Math.atan2(tananum, tanadenom));
  return alpha;     // in degrees
}
#ce
Func _calcSunRtAscension ( $t = 0 )
    Local $e = _calcObliquityCorrection ( $t )
    Local $lambda = _calcSunApparentLong ( $t )
    Local $tananum = ( Cos ( _degToRad ( $e ) ) * Sin ( _degToRad ( $lambda ) ) )
    Local $tanadenom = ( Cos ( _degToRad ( $lambda ) ) )
    Local $alpha = _radToDeg ( _atan2 ( $tananum, $tanadenom ) )
    Return $alpha  ; in degrees
EndFunc
#cs




function calcSunDeclination(t)
{
  var e = calcObliquityCorrection(t);
  var lambda = calcSunApparentLong(t);
  var sint = Math.sin(degToRad(e)) * Math.sin(degToRad(lambda));
  var theta = radToDeg(Math.asin(sint));
  return theta;     // in degrees
}
#ce
Func _calcSunDeclination ( $t = 0 )
    Local $e = _calcObliquityCorrection ( $t )
    Local $lambda = _calcSunApparentLong ( $t )
    Local $sint = Sin ( _degToRad ( $e ) ) * Sin ( _degToRad ( $lambda ) )
    Local $theta = _radToDeg ( ASin ( $sint ) )
    return $theta  ; in degrees
EndFunc
#cs




function calcEquationOfTime(t)
{
  var epsilon = calcObliquityCorrection(t);
  var l0 = calcGeomMeanLongSun(t);
  var e = calcEccentricityEarthOrbit(t);
  var m = calcGeomMeanAnomalySun(t);
  var y = Math.tan(degToRad(epsilon)/2.0);
  y *= y;
  var sin2l0 = Math.sin(2.0 * degToRad(l0));
  var sinm   = Math.sin(degToRad(m));
  var cos2l0 = Math.cos(2.0 * degToRad(l0));
  var sin4l0 = Math.sin(4.0 * degToRad(l0));
  var sin2m  = Math.sin(2.0 * degToRad(m));
  var Etime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m;
  return radToDeg(Etime)*4.0;   // in minutes of time
}
#ce
Func _calcEquationOfTime ( $t = 0 )
    Local $epsilon = _calcObliquityCorrection ( $t )
    Local $l0 = _calcGeomMeanLongSun ( $t )
    Local $e = _calcEccentricityEarthOrbit ( $t )
    Local $m = _calcGeomMeanAnomalySun ( $t )
    Local $y = Tan ( _degToRad ( $epsilon ) / 2.0 )
    $y *= $y
    Local $sin2l0 = Sin ( 2.0 * _degToRad ( $l0 ) )
    Local $sinm = Sin ( _degToRad ( $m ) )
    Local $cos2l0 = Cos ( 2.0 * _degToRad ( $l0 ) )
    Local $sin4l0 = Sin ( 4.0 * _degToRad ( $l0 ) )
    Local $sin2m = Sin ( 2.0 * _degToRad ( $m ) )
    Local $Etime = $y * $sin2l0 - 2.0 * $e * $sinm + 4.0 * $e * $y * $sinm * $cos2l0 - 0.5 * $y * $y * $sin4l0 - 1.25 * $e * $e * $sin2m
    Return _radToDeg ( $Etime ) * 4.0  ; in minutes of time
EndFunc
#cs




function calcHourAngleSunrise(lat, solarDec)
{
  var latRad = degToRad(lat);
  var sdRad  = degToRad(solarDec);
  var HAarg = (Math.cos(degToRad(90.833))/(Math.cos(latRad)*Math.cos(sdRad))-Math.tan(latRad) * Math.tan(sdRad));
  var HA = Math.acos(HAarg);
  return HA;        // in radians (for sunset, use -HA)
}
#ce
;SunRise = Horizon+0.8333
Func _calcHourAngleSunrise( $lat, $solarDec )
    Local $latRad = _degToRad ( $lat )
    Local $sdRad  = _degToRad ( $solarDec )
    Local $HAarg = ( Cos ( _degToRad ( 90.833 ) ) / ( Cos ( $latRad ) * Cos ( $sdRad ) ) - Tan ( $latRad ) * Tan ( $sdRad ) )
    Local $HA = ACos ( $HAarg )
    Return $HA  ;  in radians for morning (for evening, use -HA)
EndFunc




;CivilTwilight = (Horizon+6) < SunCenter < (Horizon+0.8333)
Func _calcHourAngleCivilTwilight( $lat, $solarDec )
    Local $latRad = _degToRad ( $lat )
    Local $sdRad  = _degToRad ( $solarDec )
    Local $HAarg = ( Cos ( _degToRad ( 96.0 ) ) / ( Cos ( $latRad ) * Cos ( $sdRad ) ) - Tan ( $latRad ) * Tan ( $sdRad ) )
    Local $HA = ACos ( $HAarg )
    Return $HA  ;  in radians for morning (for evening, use -HA)
EndFunc




;NauticalTwilight = (Horizon+12) < SunCenter < (Horizon+6)
Func _calcHourAngleNauticalTwilight( $lat, $solarDec )
    Local $latRad = _degToRad ( $lat )
    Local $sdRad  = _degToRad ( $solarDec )
    Local $HAarg = ( Cos ( _degToRad ( 102.0 ) ) / ( Cos ( $latRad ) * Cos ( $sdRad ) ) - Tan ( $latRad ) * Tan ( $sdRad ) )
    Local $HA = ACos ( $HAarg )
    Return $HA  ;  in radians for morning (for evening, use -HA)
EndFunc




;AstronimicalTwilight = (Horizon+18) < SunCenter < (Horizon+12)
Func _calcHourAngleAstronimicalTwilight( $lat, $solarDec )
    Local $latRad = _degToRad ( $lat )
    Local $sdRad  = _degToRad ( $solarDec )
    Local $HAarg = ( Cos ( _degToRad ( 108.0 ) ) / ( Cos ( $latRad ) * Cos ( $sdRad ) ) - Tan ( $latRad ) * Tan ( $sdRad ) )
    Local $HA = ACos ( $HAarg )
    Return $HA  ;  in radians for morning ( for evening, use -HA )
EndFunc
#cs




function isNumber(inputVal)
{
  var oneDecimal = false;
  var inputStr = "" + inputVal;
  for (var i = 0; i < inputStr.length; i++)
  {
    var oneChar = inputStr.charAt(i);
    if (i == 0 && (oneChar == "-" || oneChar == "+"))
    {
      continue;
    }
    if (oneChar == "." && !oneDecimal)
    {
      oneDecimal = true;
      continue;
    }
    if (oneChar < "0" || oneChar > "9")
    {
      return false;
    }
  }
  return true;
}
#ce
Func _isNumber ( $inputval )
    Return ( IsFloat ( $inputval ) Or IsNumber ( $inputval ) Or IsInt ( $inputval ) )
EndFunc
#cs




function zeroPad(n, digits) {
  n = n.toString();
  while (n.length < digits) {
    n = '0' + n;
  }
  return n;
}
#ce
Func _zeroPad ( $n, $digits )
    $n = String ( $n )
    While StringLen ( $n ) < $digits
        $n = "0" & $n
    WEnd
    Return $n
EndFunc
#cs




function month(name, numdays, abbr)
{
  this.name = name;
  this.numdays = numdays;
  this.abbr = abbr;
}
var monthList = new Array();
var i = 0;
monthList[i++] = new month("January", 31, "Jan");
monthList[i++] = new month("February", 28, "Feb");
monthList[i++] = new month("March", 31, "Mar");
monthList[i++] = new month("April", 30, "Apr");
monthList[i++] = new month("May", 31, "May");
monthList[i++] = new month("June", 30, "Jun");
monthList[i++] = new month("July", 31, "Jul");
monthList[i++] = new month("August", 31, "Aug");
monthList[i++] = new month("September", 30, "Sep");
monthList[i++] = new month("October", 31, "Oct");
monthList[i++] = new month("November", 30, "Nov");
monthList[i++] = new month("December", 31, "Dec");
#ce
; Global variable defined at top of file...
#cs




function getJD()
{
  var docmonth = document.getElementById("mosbox").selectedIndex + 1
  var docday =   document.getElementById("daybox").selectedIndex + 1
  var docyear =  readTextBox("yearbox", 5, 1, 0, -2000, 3000, 2009)
  if ( (isLeapYear(docyear)) && (docmonth == 2) ) {
    if (docday > 29) {
      docday = 29
      document.getElementById("daybox").selectedIndex = docday - 1
    }
  } else {
    if (docday > monthList[docmonth-1].numdays) {
      docday = monthList[docmonth-1].numdays
      document.getElementById("daybox").selectedIndex = docday - 1
    }
  }
  if (docmonth <= 2) {
    docyear -= 1
    docmonth += 12
  }
  var A = Math.floor(docyear/100)
  var B = 2 - A + Math.floor(A/4)
  var JD = Math.floor(365.25*(docyear + 4716)) + Math.floor(30.6001*(docmonth+1)) + docday + B - 1524.5
  return JD
}
#ce
Func _getJD ( $CompJDMon = "", $CompJDDay = "", $CompJDYear = "" )
    If ( ( _isLeapYear ( $CompJDYear ) ) And ( $CompJDMon = 2 ) ) Then
        If $CompJDDay > 29 Then
          $CompJDDay = 29
        EndIf
    Else
        If $CompJDDay > $MonthLength[$CompJDMon] Then
          $CompJDDay = $MonthLength[$CompJDMon]
        EndIf
    EndIf
    If $CompJDMon <= 2 Then
        $CompJDYear -= 1
        $CompJDMon += 12
    EndIf
    Local $A = Floor ( $CompJDYear / 100 )
    Local $B = 2 - $A + Floor ( $A / 4 )
    Local $JD = Floor ( 365.25 * ( $CompJDYear + 4716 ) ) + Floor ( 30.6001 * ( $CompJDMon + 1 ) ) + $CompJDDay + $B - 1524.5
    Return $JD
EndFunc
#cs




function getTimeLocal()
{
  var dochr = readTextBox("hrbox", 2, 1, 1, 0, 23, 12)
  var docmn = readTextBox("mnbox", 2, 1, 1, 0, 59, 0)
  var docsc = readTextBox("scbox", 2, 1, 1, 0, 59, 0)
  var docpm = document.getElementById("pmbox").checked
  var docdst = document.getElementById("dstCheckbox").checked
  if ( (docpm) && (dochr < 12) ) {
    dochr += 12
  }
  if (docdst) {
    dochr -= 1
  }
  var mins = dochr * 60 + docmn + docsc/60.0
  return mins
}
#ce
Func _getTimeLocal ( $CompHour = "", $CompMin = "", $CompSec = "", $CompDST = False )
    If ( $CompDST ) Then
        $CompHour -= 1
    EndIf
    Local $mins = $CompHour * 60 + $CompMin + $CompSec / 60.0
    Return $mins
EndFunc
#cs




function calcAzEl(output, T, localtime, latitude, longitude, zone)
{
  var eqTime = calcEquationOfTime(T)
  var theta  = calcSunDeclination(T)
  if (output) {
    document.getElementById("eqtbox").value = Math.floor(eqTime*100 +0.5)/100.0
    document.getElementById("sdbox").value = Math.floor(theta*100+0.5)/100.0
  }
  var solarTimeFix = eqTime + 4.0 * longitude - 60.0 * zone
  var earthRadVec = calcSunRadVector(T)
  var trueSolarTime = localtime + solarTimeFix
  while (trueSolarTime > 1440)
  {
    trueSolarTime -= 1440
  }
  var hourAngle = trueSolarTime / 4.0 - 180.0;
  if (hourAngle < -180)
  {
    hourAngle += 360.0
  }
  var haRad = degToRad(hourAngle)
  var csz = Math.sin(degToRad(latitude)) * Math.sin(degToRad(theta)) + Math.cos(degToRad(latitude)) * Math.cos(degToRad(theta)) * Math.cos(haRad)
  if (csz > 1.0)
  {
    csz = 1.0
  } else if (csz < -1.0)
  {
    csz = -1.0
  }
  var zenith = radToDeg(Math.acos(csz))
  var azDenom = ( Math.cos(degToRad(latitude)) * Math.sin(degToRad(zenith)) )
  if (Math.abs(azDenom) > 0.001) {
    azRad = (( Math.sin(degToRad(latitude)) * Math.cos(degToRad(zenith)) ) - Math.sin(degToRad(theta))) / azDenom
    if (Math.abs(azRad) > 1.0) {
      if (azRad < 0) {
    azRad = -1.0
      } else {
    azRad = 1.0
      }
    }
    var azimuth = 180.0 - radToDeg(Math.acos(azRad))
    if (hourAngle > 0.0) {
      azimuth = -azimuth
    }
  } else {
    if (latitude > 0.0) {
      azimuth = 180.0
    } else {
      azimuth = 0.0
    }
  }
  if (azimuth < 0.0) {
    azimuth += 360.0
  }
  var exoatmElevation = 90.0 - zenith

// Atmospheric Refraction correction

  if (exoatmElevation > 85.0) {
    var refractionCorrection = 0.0;
  } else {
    var te = Math.tan (degToRad(exoatmElevation));
    if (exoatmElevation > 5.0) {
      var refractionCorrection = 58.1 / te - 0.07 / (te*te*te) + 0.000086 / (te*te*te*te*te);
    } else if (exoatmElevation > -0.575) {
      var refractionCorrection = 1735.0 + exoatmElevation * (-518.2 + exoatmElevation * (103.4 + exoatmElevation * (-12.79 + exoatmElevation * 0.711) ) );
    } else {
      var refractionCorrection = -20.774 / te;
    }
    refractionCorrection = refractionCorrection / 3600.0;
  }

  var solarZen = zenith - refractionCorrection;

  if ((output) && (solarZen > 108.0) ) {
    document.getElementById("azbox").value = "dark"
    document.getElementById("elbox").value = "dark"
  } else if (output) {
    document.getElementById("azbox").value = Math.floor(azimuth*100 +0.5)/100.0
    document.getElementById("elbox").value = Math.floor((90.0-solarZen)*100+0.5)/100.0
    if (document.getElementById("showae").checked) {
      showLineGeodesic("#ffff00", azimuth)
    }
  }
  return (azimuth)
}
#ce
Func _calcAzEl ( $output, $T, $localtime, $latitude, $longitude, $zone )
    Local $SolarReturnAzEl[1] = [0]
    Local $SolarLightStatus, $SolarEquationOfTime, $SolarDeclination
    Local $eqTime = _calcEquationOfTime ( $T )
    Local $theta  = _calcSunDeclination ( $T )
    If $output Then
        $SolarEquationOfTime = Floor ( $eqTime * 100 + 0.5 ) / 100.0
        $SolarDeclination = Floor ( $theta * 100 + 0.5 ) / 100.0
    EndIf
    Local $solarTimeFix = $eqTime + 4.0 * $longitude - 60.0 * $zone
    Local $earthRadVec = _calcSunRadVector ( $T )
    Local $trueSolarTime = $localtime + $solarTimeFix
    While $trueSolarTime > 1440
        $trueSolarTime -= 1440
    WEnd
    Local $hourAngle = $trueSolarTime / 4.0 - 180.0;
    If $hourAngle < -180 Then
        $hourAngle += 360.0
    EndIf
    Local $haRad = _degToRad ( $hourAngle )
    Local $csz = Sin ( _degToRad ( $latitude ) ) * Sin ( _degToRad ( $theta ) ) + Cos ( _degToRad ( $latitude ) ) * Cos ( _degToRad ( $theta ) ) * Cos ( $haRad )
    If $csz > 1.0 Then
        $csz = 1.0
    ElseIf $csz < -1.0 Then
        $csz = -1.0
    EndIf
    Local $zenith = _radToDeg ( ACos ( $csz ) )
    Local $azDenom = ( Cos ( _degToRad ( $latitude ) ) * Sin ( _degToRad ( $zenith ) ) )
    If Abs ( $azDenom ) > 0.001 Then
        $azRad = ( ( Sin ( _degToRad ( $latitude ) ) * Cos ( _degToRad ( $zenith ) ) ) - Sin ( _degToRad ( $theta ) ) ) / $azDenom
        If Abs ( $azRad ) > 1.0 Then
            If $azRad < 0 Then
                $azRad = -1.0
            Else
                $azRad = 1.0
            EndIf
        EndIf
        Local $azimuth = 180.0 - _radToDeg ( ACos ( $azRad ) )
        If $hourAngle > 0.0 Then
            $azimuth = -$azimuth
        EndIf
    Else
        If $latitude > 0.0 Then
            $azimuth = 180.0
        Else
            $azimuth = 0.0
        EndIf
    EndIf
    If $azimuth < 0.0 Then
        $azimuth += 360.0
    EndIf
    Local $exoatmElevation = 90.0 - $zenith
    ; Atmospheric Refraction correction
    If $exoatmElevation > 85.0 Then
        Local $refractionCorrection = 0.0
    Else
        Local $te = Tan ( _degToRad ( $exoatmElevation ) )
        If $exoatmElevation > 5.0 Then
            Local $refractionCorrection = 58.1 / $te - 0.07 / ( $te * $te * $te ) + 0.000086 / ( $te * $te * $te * $te * $te )
        ElseIf $exoatmElevation > -0.575 Then
            Local $refractionCorrection = 1735.0 + $exoatmElevation * ( -518.2 + $exoatmElevation * ( 103.4 + $exoatmElevation * ( -12.79 + $exoatmElevation * 0.711) ) )
        Else
            Local $refractionCorrection = -20.774 / $te
        EndIf
        $refractionCorrection = $refractionCorrection / 3600.0
    EndIf

    Local $solarZen = $zenith - $refractionCorrection

    If $solarZen > 108.0 Then
        $SolarLightStatus = "Night"
    ElseIf ( ( 108.0 > $solarZen ) And ( $solarZen >= 102.0 ) ) Then
        $SolarLightStatus = "Astronomical Twilight"
    ElseIf ( ( 102.0 > $solarZen ) And ( $solarZen >= 96.0 ) ) Then
        $SolarLightStatus = "Nautical Twilight"
    ElseIf ( ( 96.0 > $solarZen ) And ( $solarZen >= 90.8333 ) ) Then
        $SolarLightStatus = "Civil Twilight"
    Else
        $SolarLightStatus = "Day"
    EndIf

    _ArrayAdd ( $SolarReturnAzEl, Floor ( ( ( $azimuth * 100 ) + 0.5 ) / 100.0 ) )
    _ArrayAdd ( $SolarReturnAzEl, Floor ( ( 90.0 - $solarZen ) * 100 + 0.5) / 100.0 )
    _ArrayAdd ( $SolarReturnAzEl, $SolarLightStatus )
    _ArrayDelete ( $SolarReturnAzEl, 0 )

    Return ( $SolarReturnAzEl )
EndFunc
#cs




function calcSolNoon(jd, longitude, timezone, dst)
{
  var tnoon = calcTimeJulianCent(jd - longitude/360.0)
  var eqTime = calcEquationOfTime(tnoon)
  var solNoonOffset = 720.0 - (longitude * 4) - eqTime // in minutes
  var newt = calcTimeJulianCent(jd + solNoonOffset/1440.0)
  eqTime = calcEquationOfTime(newt)
  solNoonLocal = 720 - (longitude * 4) - eqTime + (timezone*60.0)// in minutes
  if(dst) solNoonLocal += 60.0
  document.getElementById("noonbox").value = timeString(solNoonLocal, 3)
}
#ce
Func _calcSolNoon ( $jd, $longitude, $timezone, $dst )
    Local $tnoon = _calcTimeJulianCent ( $jd - $longitude / 360.0 )
    Local $eqTime = _calcEquationOfTime ( $tnoon )
    Local $solNoonOffset = 720.0 - ( $longitude * 4 ) - $eqTime  ; in minutes
    Local $newt = _calcTimeJulianCent ( $jd + $solNoonOffset / 1440.0 )
    $eqTime = _calcEquationOfTime ( $newt )
    $solNoonLocal = 720 - ( $longitude * 4 ) - $eqTime + ( $timezone * 60.0 )  ; in minutes
    If $dst Then $solNoonLocal += 60.0
    Return _timeString ( $solNoonLocal, 3 )
EndFunc
#cs




function dayString(jd, next, flag)
{
// returns a string in the form DDMMMYYYY[ next] to display prev/next rise/set
// flag=2 for DD MMM, 3 for DD MM YYYY, 4 for DDMMYYYY next/prev
  if ( (jd < 900000) || (jd > 2817000) ) {
    var output = "error"
  } else {
  var z = Math.floor(jd + 0.5);
  var f = (jd + 0.5) - z;
  if (z < 2299161) {
    var A = z;
  } else {
    alpha = Math.floor((z - 1867216.25)/36524.25);
    var A = z + 1 + alpha - Math.floor(alpha/4);
  }
  var B = A + 1524;
  var C = Math.floor((B - 122.1)/365.25);
  var D = Math.floor(365.25 * C);
  var E = Math.floor((B - D)/30.6001);
  var day = B - D - Math.floor(30.6001 * E) + f;
  var month = (E < 14) ? E - 1 : E - 13;
  var year = ((month > 2) ? C - 4716 : C - 4715);
  if (flag == 2)
    var output = zeroPad(day,2) + " " + monthList[month-1].abbr;
  if (flag == 3)
    var output = zeroPad(day,2) + monthList[month-1].abbr + year.toString();
  if (flag == 4)
    var output = zeroPad(day,2) + monthList[month-1].abbr + year.toString() + ((next) ? " next" : " prev");
  }
  return output;
}
#ce
Func _dayString ( $jd, $next, $flag )
; returns a string in the form DDMMMYYYY[ next] to display prev/next rise/set
; flag=2 for DD MMM, 3 for DD MM YYYY, 4 for DDMMYYYY next/prev
    If ( $jd < 900000 ) Or ( $jd > 2817000 ) Then
        Local $output = "error"
    Else
        Local $z = Floor ( $jd + 0.5 )
        Local $f = ( $jd + 0.5 ) - $z
        If $z < 2299161 Then
            $A = $z
        Else
            Local $alpha = Floor ( ( $z - 1867216.25 ) / 36524.25 )
            $A = $z + 1 + $alpha - Floor ( $alpha / 4 )
        EndIf
        Local $B = $A + 1524
        Local $C = Floor ( ( $B - 122.1 ) / 365.25 )
        Local $D = Floor ( 365.25 * $C )
        Local $E = Floor ( ( $B - $D ) / 30.6001 )
        Local $day = $B - $D - Floor ( 30.6001 * $E ) + $f
        Local $month
        If $E < 14 Then
            $month  = $E - 1
        Else
            $month = $E - 13
        EndIf
        Local $year
        If $month > 2 Then
            $year = $C - 4716
        Else
            $year = $C - 4715
        EndIf
        If $flag = 2 Then
            Local $output = _zeroPad ( $day, 2 ) & " " & $MonthAbrev[$month-1]
        ElseIf $flag = 3 Then
            Local $output = _zeroPad ( $day, 2 ) & $MonthAbrev[$month-1] & String ( $year )
        ElseIf $flag = 4 Then
            If $next Then
                Local $output = _zeroPad ( $day, 2 ) & $MonthAbrev[$month-1] & String ( $year ) & " next"
            Else
                Local $output = _zeroPad ( $day, 2 ) & $MonthAbrev[$month-1] & String ( $year ) & " prev"
            EndIf
        EndIf
    EndIf
    Return $output
EndFunc
#cs




function timeDateString(JD, minutes)
{
  var output = timeString(minutes, 2) + " " + dayString(JD, 0, 2);
  return output;
}
#ce
Func _TimeDateString ( $JD, $minutes )
    Local $output = _timeString ( $minutes, 2 ) & " " & _dayString ( $JD, 0, 2 )
    Return $output
EndFunc
#cs




function timeString(minutes, flag)
// timeString returns a zero-padded string (HH:MM:SS) given time in minutes
// flag=2 for HH:MM, 3 for HH:MM:SS
{
  if ( (minutes >= 0) && (minutes < 1440) ) {
  var floatHour = minutes / 60.0;
  var hour = Math.floor(floatHour);
  var floatMinute = 60.0 * (floatHour - Math.floor(floatHour));
  var minute = Math.floor(floatMinute);
  var floatSec = 60.0 * (floatMinute - Math.floor(floatMinute));
  var second = Math.floor(floatSec + 0.5);
  if (second > 59) {
    second = 0
    minute += 1
  }
  if ((flag == 2) && (second >= 30)) minute++;
  if (minute > 59) {
    minute = 0
    hour += 1
  }
  var output = zeroPad(hour,2) + ":" + zeroPad(minute,2);
  if (flag > 2) output = output + ":" + zeroPad(second,2);
  } else {
    var output = "error"
  }
  return output;
}
#ce
Func _timeString ( $minutes, $flag )
    ; timeString returns a zero-padded string (HH:MM:SS) given time in minutes
    ; flag=2 for HH:MM, 3 for HH:MM:SS
    If ( ( $minutes >= 0 ) And ( $minutes < 1440 ) ) Then
        Local $floatHour = $minutes / 60.0
        Local $hour = Floor ( $floatHour )
        Local $floatMinute = 60.0 * ( $floatHour - Floor ( $floatHour) )
        Local $minute = Floor ( $floatMinute )
        Local $floatSec = 60.0 * ( $floatMinute - Floor ( $floatMinute ) )
        Local $second = Floor ( $floatSec + 0.5 )
        If ( $second > 59) Then
            $second = 0
            $minute += 1
        EndIf
        If ( ( $flag = 2 ) And ( $second >= 30 ) ) Then $minute += 1
        If ( $minute > 59 ) Then
            $minute = 0
            $hour += 1
        EndIf
        Local $output = _zeroPad ( $hour, 2 ) & ":" & _zeroPad ( $minute, 2 )
        If $flag > 2 Then $output = $output & ":" & _zeroPad ( $second, 2 )
    Else
        $output = "error"
    EndIf
    Return $output
EndFunc
#cs




function calcSunriseSetUTC(rise, JD, latitude, longitude)
{
  var t = calcTimeJulianCent(JD);
  var eqTime = calcEquationOfTime(t);
  var solarDec = calcSunDeclination(t);
  var hourAngle = calcHourAngleSunrise(latitude, solarDec);
  //alert("HA = " + radToDeg(hourAngle));
  if (!rise) hourAngle = -hourAngle;
  var delta = longitude + radToDeg(hourAngle);
  var timeUTC = 720 - (4.0 * delta) - eqTime;   // in minutes
  return timeUTC
}
#ce
Func _calcSunriseSetUTC ( $rise, $JD, $latitude, $longitude )
    Local $t = _calcTimeJulianCent ( $JD )
    Local $eqTime = _calcEquationOfTime ( $t )
    Local $solarDec = _calcSunDeclination ( $t )
    Local $RiseSetAngle = _calcHourAngleSunrise ( $latitude, $solarDec )
    If Not $rise Then $RiseSetAngle = -$RiseSetAngle
    Local $delta = $longitude + _radToDeg ( $RiseSetAngle )
    Local $timeUTCRS = 720 - ( 4.0 * $delta ) - $eqTime
    Return $timeUTCRS  ; in minutes
EndFunc
#cs




function calcSunriseSet(rise, JD, latitude, longitude, timezone, dst)
// rise = 1 for sunrise, 0 for sunset
{
  var id = ((rise) ? "risebox" : "setbox")
  var timeUTC = calcSunriseSetUTC(rise, JD, latitude, longitude);
  var newTimeUTC = calcSunriseSetUTC(rise, JD + timeUTC/1440.0, latitude, longitude);
  if (isNumber(newTimeUTC)) {
    var timeLocal = newTimeUTC + (timezone * 60.0)
    if (document.getElementById(rise ? "showsr" : "showss").checked) {
      var riseT = calcTimeJulianCent(JD + newTimeUTC/1440.0)
      var riseAz = calcAzEl(0, riseT, timeLocal, latitude, longitude, timezone)
      showLineGeodesic(rise ? "#66ff00" : "#ff0000", riseAz)
    }
    timeLocal += ((dst) ? 60.0 : 0.0);
    if ( (timeLocal >= 0.0) && (timeLocal < 1440.0) ) {
      document.getElementById(id).value = timeString(timeLocal,2)
    } else  {
      var jday = JD
      var increment = ((timeLocal < 0) ? 1 : -1)
      while ((timeLocal < 0.0)||(timeLocal >= 1440.0)) {
        timeLocal += increment * 1440.0
    jday -= increment
      }
      document.getElementById(id).value = timeDateString(jday,timeLocal)
    }
  } else { // no sunrise/set found
    var doy = calcDoyFromJD(JD)
    if ( ((latitude > 66.4) && (doy > 79) && (doy < 267)) ||
    ((latitude < -66.4) && ((doy < 83) || (doy > 263))) )
    {   //previous sunrise/next sunset
      if (rise) { // find previous sunrise
        jdy = calcJDofNextPrevRiseSet(0, rise, JD, latitude, longitude, timezone, dst)
      } else { // find next sunset
        jdy = calcJDofNextPrevRiseSet(1, rise, JD, latitude, longitude, timezone, dst)
      }
      document.getElementById(((rise)? "risebox":"setbox")).value = dayString(jdy,0,3)
    } else {   //previous sunset/next sunrise
      if (rise == 1) { // find previous sunrise
        jdy = calcJDofNextPrevRiseSet(1, rise, JD, latitude, longitude, timezone, dst)
      } else { // find next sunset
        jdy = calcJDofNextPrevRiseSet(0, rise, JD, latitude, longitude, timezone, dst)
      }
      document.getElementById(((rise)? "risebox":"setbox")).value = dayString(jdy,0,3)
    }
  }
}
#ce
Func _calcSunriseSet ( $rise, $JD, $latitude, $longitude, $timezone, $dst )
; rise = 1 for sunrise, 0 for sunset
    Local $timeUTC = _calcSunriseSetUTC ( $rise, $JD, $latitude, $longitude )
    Local $newTimeUTC = _calcSunriseSetUTC ( $rise, $JD + $timeUTC / 1440.0, $latitude, $longitude )
    If isNumber ( $newTimeUTC ) Then
        Local $timeLocal = $newTimeUTC + ( $timezone * 60.0 )
        If $rise Then
            Local $riseT = _calcTimeJulianCent ( $JD + $newTimeUTC / 1440.0 )
            Local $riseAz = _calcAzEl ( 0, $riseT, $timeLocal, $latitude, $longitude, $timezone )
        EndIf
        If $dst Then
            $timeLocal += 60.0
        Else
            $timelocal += 0.0
        EndIf
        If ( ( $timeLocal >= 0.0) And ( $timeLocal < 1440.0) ) Then
            Return _timeString ( $timeLocal, 2 )
        Else
            Local $jday = $JD
            Local $increment
            If $timeLocal < 0 Then
                $increment = 1
            Else
                $increment = -1
            EndIf
            While ( ( $timeLocal < 0.0 ) Or ( $timeLocal >= 1440.0 ) )
                $timeLocal += $increment * 1440.0
                $jday -= $increment
            WEnd
            Return _timeDateString ( $jday, $timeLocal )
        EndIf
    Else
        Local $doy = _calcDoyFromJD ( $JD )
        If ( ( ( $latitude > 66.4 ) And ( $doy > 79 ) And ( $doy < 267 ) ) Or ( ( $latitude < -66.4 ) And ( ( $doy < 83 ) Or ( $doy > 263 ) ) ) ) Then
            If ( $rise ) Then
                $jdy = _calcJDofNextPrevRiseSet ( 0, $rise, $JD, $latitude, $longitude, $timezone, $dst )
            Else
                $jdy = _calcJDofNextPrevRiseSet ( 1, $rise, $JD, $latitude, $longitude, $timezone, $dst )
            EndIf
            ;Return _dayString ( $jdy, 0, 3 )
        Else
            If ( $rise ) Then
                $jdy = _calcJDofNextPrevRiseSet ( 1, $rise, $JD, $latitude, $longitude, $timezone, $dst )
            Else
                $jdy = _calcJDofNextPrevRiseSet ( 0, $rise, $JD, $latitude, $longitude, $timezone, $dst )
            EndIf
            ;Return _dayString ( $jdy, 0, 3 )
        EndIf
    EndIf
EndFunc
#cs




function calcJDofNextPrevRiseSet(next, rise, JD, latitude, longitude, tz, dst)
{
  var julianday = JD;
  var increment = ((next) ? 1.0 : -1.0);

  var time = calcSunriseSetUTC(rise, julianday, latitude, longitude);
  while(!isNumber(time)){
    julianday += increment;
    time = calcSunriseSetUTC(rise, julianday, latitude, longitude);
  }
  var timeLocal = time + tz * 60.0 + ((dst) ? 60.0 : 0.0)
  while ((timeLocal < 0.0) || (timeLocal >= 1440.0))
  {
    var incr = ((timeLocal < 0) ? 1 : -1)
    timeLocal += (incr * 1440.0)
    julianday -= incr
  }
  return julianday;
}
#ce
Func _calcJDofNextPrevRiseSet ( $next, $rise, $JD, $latitude, $longitude, $tz, $dst )
    Local $julianday = $JD
    Local $increment
    If $next Then
        $increment = 1.0
    Else
        $increment = -1.0
    EndIf
    Local $time = _calcSunriseSetUTC ( $rise, $julianday, $latitude, $longitude )
    While Not isNumber ( $time )
        $julianday += $increment
        $time = _calcSunriseSetUTC ( $rise, $julianday, $latitude, $longitude )
    WEnd
    Local $timeLocal
    If $dst Then
        $timeLocal = $time + $tz * 60.0 + 60.0
    Else
        $timeLocal = $time + $tz * 60.0
    EndIf
    While ( ( $timeLocal < 0.0 ) Or ( $timeLocal >= 1440.0 ) )
        Local $incr
        If $timeLocal < 0 Then
            $incr = 1
        Else
            $incr = -1
        EndIf
        $timeLocal += ( $incr * 1440.0 )
        $julianday -= $incr
    WEnd
    Return $julianday
EndFunc
#cs




function calculate() {
  //refreshMap()
  //clearOutputs()
  //map.clearOverlays()
  //showMarkers()
  var jday = getJD()
  var tl = getTimeLocal()
  var tz = readTextBox("zonebox", 5, 0, 0, -14, 13, 0)
  var dst = document.getElementById("dstCheckbox").checked
  var total = jday + tl/1440.0 - tz/24.0
  var T = calcTimeJulianCent(total)
  var lat = parseFloat(document.getElementById("latbox").value.substring(0,9))
  var lng = parseFloat(document.getElementById("lngbox").value.substring(0,10))
  calcAzEl(1, T, tl, lat, lng, tz)
  calcSolNoon(jday, lng, tz, dst)
  var rise = calcSunriseSet(1, jday, lat, lng, tz, dst)
  var set  = calcSunriseSet(0, jday, lat, lng, tz, dst)
  //alert("JD " + jday + "  " + rise + "  " + set + "  ")
}
#ce
;( $AstroLat, $AstroLong, $AstroMonth, $AstroDay, $AstroYear, $AstroTimeZone, $AstroHour, $AstroMinute, $AstroSecond, $AstroDST )
Func _SolarCalculate ( $CalcLat, $CalcLong, $CalcMonth, $CalcDay, $CalcYear, $CalcTz, $CalcHour, $CalcMinute, $CalcSecond, $CalcDST, $CalcCurrent )
    Local $SolarCalculateResult[1] = [0]
    Local $jday = _getJD ( $CalcMonth, $CalcDay, $CalcYear )
    Local $tl = _getTimeLocal ( $CalcHour, $CalcMinute, $CalcSecond, $CalcDST )
    Local $tz
    If $CalcTz = "" Then
        $tz = Int ( $CalcLong / 15 )
    Else
        $tz = $CalcTz
    EndIf
    Local $total = $jday + $tl / 1440.0 - $tz / 24.0
    Local $T = _calcTimeJulianCent ( $total )
    Local $CalcSolNoon = _calcSolNoon ( $jday, $CalcLong, $tz, $CalcDST ) ; Returns String: "HH:MM:SS"
    ;MsgBox ( 0, "Solar Noon:", $CalcSolNoon )
    Local $rise = _calcSunriseSet ( 1, $jday, $CalcLat, $CalcLong, $tz, $CalcDST ) ; Returns String: "HH:MM" or "HH:MM DD Mon"
    ;MsgBox ( 0, "Sun Rise:", $rise )
    Local $set  = _calcSunriseSet ( 0, $jday, $CalcLat, $CalcLong, $tz, $CalcDST ) ; Returns String: "HH:MM" or "HH:MM DD Mon"
    ;MsgBox ( 0, "Sun Set:", $set )
    _ArrayAdd ( $SolarCalculateResult, $rise )
    _ArrayAdd ( $SolarCalculateResult, $CalcSolNoon )
    _ArrayAdd ( $SolarCalculateResult, $set )
    $SolarCalculateResult[0] += 3
    If $CalcCurrent Then
        Local $CalcAzEl = _calcAzEl ( 1, $T, $tl, $CalcLat, $CalcLong, $tz ) ; Returns Array: 0 - Azimuth, 1 - Elevation, 2 - Illumination-Disposition
        _ArrayConcatenate ( $SolarCalculateResult, $CalcAzEl )
        $SolarCalculateResult[0] += 3
    EndIf
    Return $SolarCalculateResult
EndFunc
#cs
</SCRIPT>
#ce

Here is an example usage script:

#include <SunPosition.au3>
#include <GUIConstantsEx.au3>
#include <WindowsConstants.au3>
Global $MyLat = "37.818599"
Global $MyLong = "-122.478418"

$SunPosition_Window = GUICreate ( "Sun Position", 220, 180 )

GUICtrlCreateLabel ( "Sun Rise:", 10, 12 )
$Rise_Time = GUICtrlCreateInput("", 100, 10 )
GUICtrlCreateLabel ( "Solar Noon:", 10, 32 )
$Noon_Time = GUICtrlCreateInput("", 100, 30 )
GUICtrlCreateLabel ( "Sun Set:", 10, 52 )
$Set_Time = GUICtrlCreateInput("", 100, 50 )
GUICtrlCreateLabel ( "Azimuth:", 10, 72 )
$Sun_Azimuth = GUICtrlCreateInput("", 100, 70 )
GUICtrlCreateLabel ( "Elevation:", 10, 92 )
$Sun_Elevation = GUICtrlCreateInput("", 100, 90 )
GUICtrlCreateLabel ( "Illum Disp.:", 10, 112 )
$Light_Disposition = GUICtrlCreateInput("", 100, 110, 110 )
GUICtrlCreateLabel ( "Dimmer Lvl:", 10, 135 )
$Light_Level_Control = GUICtrlCreateLabel ( "", 70, 135, 23, 17 )
GUICtrlCreateLabel ( "Switched Lts:", 100, 135 )
$Light_Switch_Control = GUICtrlCreateLabel ( "", 170, 135, 23, 17 )
$Dummy = GUICtrlCreateInput("", 100, 210 )
GUICtrlCreateLabel ( "Light Blocking Curtians:", 10, 155 )
$Privacy_Curtain_Control = GUICtrlCreateLabel ( "", 125, 155, 23, 17 )

GUISetState(@SW_SHOW)

Global $RefreshCounter = 0
Global $Light_Level_Register = 0
Global $Curtain_Opening_Register = 0
Global $CalcLightMagnatude = 0
While 1
    Sleep ( 2 )
    $msg = GUIGetMsg()
    Switch $msg
        Case $GUI_EVENT_CLOSE
            ExitLoop
        Case ""
            $RefreshCounter += 1
            If $RefreshCounter >=25 Then
                $TestVector = _GetSolarAstronimicalData ( $MyLat, $MyLong, @MON, @MDAY, @YEAR, Int ( $MyLong / 15 ), @HOUR, @MIN, @SEC, True )
                GUICtrlSetData ( $Rise_Time, $TestVector[1] )
                GUICtrlSetData ( $Noon_Time, $TestVector[2] )
                GUICtrlSetData ( $Set_Time, $TestVector[3] )
                GUICtrlSetData ( $Sun_Azimuth, $TestVector[4] )
                GUICtrlSetData ( $Sun_Elevation, $TestVector[5] )
                GUICtrlSetData ( $Light_Disposition, $TestVector[6] )

                ; 0.8333 add to actual elevation for lights
                ; 1.722233333 subtract from actual elevation for lights
                ; 1.722233333 x 3 = 5.166699999 Max for Light control
                $Light_Level_Register = $TestVector[5] + 0.8333
                $Light_Level_Register -= 3.444466666
                $Light_Level_Register *= -1
                If $Light_Level_Register < 0 Then
                    $Light_Level_Register = 0
                ElseIf $Light_Level_Register > 5.166699999  Then
                    $Light_Level_Register = 255
                Else
                    $Light_Level_Register = $Light_Level_Register / 5.166699999 
                    $Light_Level_Register = int ( 255 * $Light_Level_Register )
                EndIf
                $CalcLightMagnatude = 0
                $CalcLightMagnatude = ( $Light_Level_Register * 65536 ) + ( $Light_Level_Register * 256 ) + Round ( ( $Light_Level_Register * 0.75 ), 0 )
                GUICtrlSetBkColor ( $Light_Level_Control, $CalcLightMagnatude )
                GUICtrlSetData ( $Light_Level_Control, $Light_Level_Register )
                If $Light_Level_Register < 191 Then
                    GUICtrlSetColor ( $Light_Level_Control, 0xffffff )
                Else
                    GUICtrlSetColor ( $Light_Level_Control, 0x000000 )
                EndIf
                GUICtrlSetState( $Dummy, $GUI_FOCUS)
                If $Light_Level_Register > 95 Then
                    GUICtrlSetData ( $Light_Switch_Control, " On" )
                    GUICtrlSetBkColor ( $Light_Switch_Control, 0x00ff00 )
                Else
                    GUICtrlSetData ( $Light_Switch_Control, " Off" )
                    GUICtrlSetBkColor ( $Light_Switch_Control, 0xff0000 )
                EndIf
                
                ; 0.8333 + 3.444466666 = 4.277766666 add to actual elevation for curtains
                ; 1.722233333 x 2 = 3.444466666 Max for Curtain control
                $Curtain_Opening_Register = $TestVector[5] + 4.277766666
                $Curtain_Opening_Register *= -1
                If $Curtain_Opening_Register < 0 Then
                    $Curtain_Opening_Register = 0
                ElseIf $Curtain_Opening_Register > 3.444466666  Then
                    $Curtain_Opening_Register = 255
                Else
                    $Curtain_Opening_Register = $Curtain_Opening_Register / 3.444466666 
                    $Curtain_Opening_Register = int ( 255 * $Curtain_Opening_Register )
                EndIf
                $Curtain_Opening_Register = 255 - $Curtain_Opening_Register
                GUICtrlSetData ( $Privacy_Curtain_Control, $Curtain_Opening_Register )
                
                $RefreshCounter = 0
            EndIf
    EndSwitch
WEnd
GUIDelete()

I personally believe this data can be used for pointing solar panels or, in home automation, opening and closing blinds/curtains to avoid over-exposure of a room. It could also be used for scheduling lights, etc…

Have fun with it!

-Tim

Edited April 8, 2011 by tim292stro

ldub · March 22, 2011

Interesting, Thanks,

Edited March 22, 2011 by ldub

atticdweller · January 23, 2014

Hi Tim,

I started using your Solar Position UDF. I'm getting some errors in the Azimuth values. Some days give repeat Azimuth values for consecutive days. I plotted out the change in the Azimuth each day at noon between the Equinoxes of this year. The zero values indicate that the azimuth value at noon was the same for 2 days in a row. Kind of a strange error in continuity for the data.

If anyone else has experienced this, let me know. Thanks!

tim292stro · July 21, 2014

Been a while since I have used AutoIt (and this site, almost forgot what I used as a password) - I have personnaly moved to Linux and thus Perl.

I'll take a look at the code and see what caused this - do you know the date range for the above graph? I'm guessing there is some kind of overflow - there is also an sublte shift in what looks like month blocks. Both effects are probably another bug yet to be caught, and from looking at the data, I'm guessing it's somewhere in the conversion from calendar days to Julian days.

jaberwacky · July 21, 2014

I think I do spot one error. Forgive me if I'm wrong. (Probably not related to the issue above.)

function isLeapYear(yr)
{
  return ((yr % 4 == 0 && yr % 100 != 0) || yr % 400 == 0);
}

Func _isLeapYear ( $yr = 0 )
    Return ( ( Mod ( $yr, 4 ) And Mod ($yr, 100 ) ) Or ( Mod ( $yr, 400 ) = 0 ) ) ; Mod($yr, 4) does not test for 0
    Return ( ( Mod ( $yr, 4 ) = 0 And Mod ($yr, 100 ) <> 0) Or ( Mod ( $yr, 400 ) = 0 ) ) ; Should this be the new function?
EndFunc

Or better yet, I believe there is a native AutoIt function for this.

Edited July 21, 2014 by jaberwacky

tim292stro · July 21, 2014

Yes, in the Date.au3 there is a function called _DateIsLeapYear(), and would probably have been half the code.

The logical test of "modulo 4 current year" will either provide a zero (false) or a positive non-zero (true). What I did was short hand since it was providing a value to a boolean test (the AND), but yes you could also check if it's numerically equal to zero (this is an extra "test" step). Note that both the modulo math values must be true (positive non-zero) for the AND to return true as well, I may have saved 6 characters here . I'll agree what I did there was "ugly" since it wasn't consistent in methodology on the same line - since I did the "= 0" version of the check after the OR, but it's still functional as intended.

If nothing else this should illustrate the point that there are many ways to do things in AutoIt.

Love the user name by the way, makes me feel like a should do a presentation about positive forward looking statements .

Edited July 21, 2014 by tim292stro

jaberwacky · July 21, 2014

presentation about positive forward looking statements

Sometimes I can be kind of numb. Would you care to explain like I'm five?

Sign In

Solar Position UDF

Recommended Posts

tim292stro

Link to comment

Share on other sites

ldub

Link to comment

Share on other sites

atticdweller

Link to comment

Share on other sites

tim292stro

Link to comment

Share on other sites

jaberwacky

Link to comment

Share on other sites

tim292stro

Link to comment

Share on other sites

jaberwacky

Link to comment

Share on other sites

Create an account or sign in to comment

Create an account

Sign in

Recently Browsing 0 members

Browse

AutoIt Resources

Release

Beta