Appendix B. Set of BASIC subroutines to calculate night length Subroutine NightLevelInit is first called to initialize the public variables. Function NightLength is then called and returns the number of hours for the sun to travel from zenith to zenith on date d at location Longitude and Latitude. Functions GetUTRise and GetUTSet can then be called to get the times of sunrise and sunset in UT (formerly GMT) hours. Source of equations are marked, for example, as [PDS 25] for Duffet-Smith (1988) page 25, or [JM 10] for Meeus (1991) page 10. Public DTR As Double 'constant degrees-to-radians Public RTD As Double 'constant radians-to-degrees Public Pi As Double 'constant 3.14159... Public TwoPi As Double 'constant 6.28318... Public DaysPerYear As Double 'constant 365.25... Public Epoch As Date 'epoch used for sun motion equations Public EclipticLongitudeAtEpoch As Double Public EclipticLongitudeOfPerigreeAtEpoch As Double Public EccentricityAtEpoch As Double Public ObliquityOfTheElliptic As Double Public UTrise As Single 'Time of sunrise in UT Public UTset As Single 'time of sunset in UT Sub NightLevelInit() DTR = 3.14159265358979 / 180 RTD = 180 / 3.14159265358979 Pi = 3.14159265358979 TwoPi = 2 * Pi DaysPerYear = 365.242191 Epoch = #12/31/1989# EclipticLongitudeAtEpoch = 279.403303 EclipticLongitudeOfPerigreeAtEpoch = 282.768422 EccentricityAtEpoch = 0.016713 ObliquityOfTheElliptic = 23.441884 End Sub Function NightLength(Latitude As Single, Longitude As Single, d As Date, _ Zenith As Single) As Single ' Calculation is based [PDS 93] with enhancements to improve accuracy Dim RA1 As Double 'Right ascension of sun on day d at 00UT Dim HA1 As Date 'Hour angle of sun on day d at 00UT Dim Declination1 As Double 'Declination of sun on day d at 00UT Dim d2 As Date 'Date of day afer d Dim RA2 As Double 'Right ascension of sun on day d+1 at 00UT Dim HA2 As Date 'Hour angle of sun on day d+1 at 00UT Dim Declination2 As Double 'Declination of sun on day d+1 at 00UT Dim N As Double 'Julian Date of day d at 00UT Dim R As Double 'Atmospheric refraction in degrees Dim LST1r As Double 'Local siderial time of sunrise calculated from sun's' position on midnitght before Dim LST1s As Double 'Local siderial time of sunset calculated from sun's position on midnitght before Dim GST1r As Double 'Grwich siderial time of sunrise calculated from sun's position on midnitght before Dim GST1s As Double 'Grwich siderial time of sunset calculated from sun's position on midnitght before Dim LST2r As Double 'Local siderial time of sunrise calculated from sun's position on midnitght after Dim LST2s As Double 'Local siderial time of sunset calculated from sun's position on midnitght after Dim GST2r As Double 'Grwich siderial time of sunrise calculated from sun's position on midnitght after Dim GST2s As Double 'Grwich siderial time of sunset calculated from sun's position on midnitght after Dim GSTr As Double 'Greenwich siderial time of sunrise calculated from the sun's interpolated position Dim GSTs As Double 'Greenwich siderial time of sunset calculated from the sun's interpolated position Dim T00 As Double 'Greenwich siderial time of 00UT on day d at Greenwich Dim T00p As Double 'Greenwich siderial time of 00UT on day d at observer's longitude Dim v As Double N = Int(d) + 2415018.5 'Calculate Julian Date R = GetRefraction(Zenith) 'When sun is actually at Zenith, it appears to be R higher Call getRAD(d, RA1, Declination1) 'For day d, get Right Ascension and Declination of sun at 00UT LST1r = GetLSTr(Latitude, Declination1, RA1, Zenith + R) 'Calculate local siderial time of sunrise at angle Zenith If LST1r < 0 Then NightLength = LST1r 'special case Exit Function End If LST1s = GetLSTs(Latitude, Declination1, RA1, Zenith + R) 'Calculate local siderial time of sunset at angle Zenith If LST1s < 0 Then NightLength = LST1s 'special case Exit Function End If GST1r = Range0_24(LST1r + Longitude / 15) 'Convert to Greenwich Siderial Time GST1s = Range0_24(LST1s + Longitude / 15) 'Convert to Greenwich Siderial Time d2 = d + 1 'For the next day, calculate the same parameters Call getRAD(d2, RA2, Declination2) LST2r = GetLSTr(Latitude, Declination2, RA2, Zenith + R) If LST2r < 0 Then NightLength = LST2r Exit Function End If LST2s = GetLSTs(Latitude, Declination2, RA2, Zenith + R) If LST2s < 0 Then NightLength = LST2s Exit Function End If GST2r = Range0_24(LST2r + Longitude / 15) GST2s = Range0_24(LST2s + Longitude / 15) If GST1r > GST2r Then GST2r = GST2r + 24 If GST1s > GST2s Then GST2s = GST2s + 24 T00 = UT_to_GST(d) 'Calculate the Greenwich siderial time at Geenwich at 00UT T00p = T00 + 1.002737909 * Longitude / 15 'Greenwich siderial time at observer at 00UT If T00p < 0 Then T00p = T00p + 24 If GST1r < T00p Then 'Use this time to update Greenich siderial times of rising, if necessary GST1r = GST1r + 24 GST2r = GST2r + 24 End If If GST1s < T00p Then 'Use this time to update Greenich siderial times of setting, if necessary GST1s = GST1s + 24 GST2s = GST2s + 24 End If GSTr = (24.065709816*GST1r - T00*(GST2r - GST1r)) / (24.065709816 + GST1r - GST2r) GSTs = (24.065709816*GST1s - T00*(GST2s - GST1s)) / (24.065709816 + GST1s - GST2s) 'interpolate times UTrise = GST_to_UT(d, GSTr) UTset = GST_to_UT(d, GSTs) NightLength = Range0_24((GSTr - GSTs) * 0.9972695663) 'night length is time difference converted from siderial time End Function Function GetUTRise() As Single GetUTRise = Utrise 'Return already calculated UT of rise (hours) End Function Function GetUTSet() As Single GetUTSet = UTset 'Return already calculated UT of set (hours) End Function Supporting Astronomical Subroutines Function GetLSTr(Latitude As Single, Declination As Double, RAscension As Double, _ Zenith As Single) As Double ' ' Get local siderial time of sun rising at angle Zenith with Declination and RightAscension ' ' [PDS 52, modified for arbitrary Zenith from JM 98] ' ' Returns special cases of -1 for sun never gets high enough in the sky ' and -2 for sun never gets low enough in the sky Dim H As Double H = (Cos(Zenith * DTR) - Sin(Latitude * DTR) * Sin(Declination * DTR)) / _ (Cos(Latitude * DTR) * Cos(Declination * DTR)) If H >= -1 And H <= 1 Then GetLSTr = 24 - RTD * ACos(H) / 15 + 24 * RAscension / 6.28318530717959 While GetLSTr > 24 GetLSTr = GetLSTr - 24 Wend ElseIf H > 0 Then GetLSTr = -1 Else GetLSTr = -2 End If End Function Function GetLSTs(Latitude As Single, Declination As Double, RAscension As Double, _ Zenith As Single) As Double ' ' Get local siderial time of sun setting at angle Zenith with Declination and RightAscension ' [PDS 52, modified for arbitrary Zenith from JM 98] ' ' Returns special cases of -1 for sun never gets high enough in the sky ' and -2 for sun never gets low enough in the sky Dim H As Double H = (Cos(Zenith * DTR) - Sin(Latitude * DTR) * Sin(Declination * DTR)) / _ (Cos(Latitude * DTR) * Cos(Declination * DTR)) If H >= -1 And H <= 1 Then GetLSTs = RTD * ACos(H) / 15 + 24 * RAscension / 6.28318530717959 While GetLSTs >= 24 GetLSTs = GetLSTs - 24 Wend ElseIf H > 0 Then GetLSTs = -1 Else GetLSTs = -2 End If End Function Sub getRAD(d1 As Date, RightAscension As Double, Declination As Double) ' ' Calculation of the position of the sun at 00UT on date d1 - returns RightAscension and Declination ' ' First we determine the Longitude of the sun (lambda) using standard values for date Epoch from [PDS 86] ' This value is calculated in the ecliptic plane's coordinate set [PDS 30] (the plane of the orbit of the ' earth around the sun), thus the Latitude of the sun (in these coordinates) is zero. ' These are converted to equitorial units [PDS 27] to facilitate determination of position of sun in earth-based ' units. Dim d As Long, N As Double, Mdot As Double Dim Ec As Double, lambda As Double Dim x As Double, y As Double d = Int(d1 - Epoch) N = d * 360 / DaysPerYear While N < 0 N = N + 360 Wend While N > 360 N = N - 360 Wend Mdot = N + EclipticLongitudeAtEpoch - EclipticLongitudeOfPerigreeAtEpoch If Mdot < 0 Then Mdot = Mdot + 360 Ec = Sin(Mdot * DTR) * EccentricityAtEpoch * 360 / Pi lambda = Mdot + Ec + EclipticLongitudeOfPerigreeAtEpoch If lambda > 360 Then lambda = lambda - 360 lambda = lambda * DTR Declination = ASin(Sin(DTR * ObliquityOfTheElliptic) * Sin(lambda)) * RTD y = Sin(lambda) * Cos(DTR * ObliquityOfTheElliptic) x = Cos(lambda) RightAscension = ATan2(x, y) End Sub Function GetRefraction(Zenith As Single) As Double ' ' Atmospheric refraction [JM 102] to be added to Zenith ' when calculating LSTr's and LSTs's. Refraction varies with Zenith angle up to a max of .553 degrees at horizon ' Dim R As Double If Zenith <= 90 Then R = 1 / (60 * Tan(DTR * (90 - Zenith + 7.31 / (94.4 - Zenith)))) GetRefraction = R - 0.06 * Sin(DTR * (14.7 * R + 13)) Else GetRefraction = 0.552687141853696 'use value for zenith = 90 End If End Function Function GST_to_UT(d As Date, GST As Double) As Double ' ' Convert Greewich Siderial Time GST on date d to Universal Time [PDS 18] ' Dim S As Double, T As Double, T0 As Double, N As Double, H As Single N = Int(d) + 2415018.5 S = N - 2451545 T = S / 36525 T0 = 6.697374558 + 2400.051336 * T + 0.000025862 * T * T T0 = Range0_24(GST - T0) GST_to_UT = T0 * 0.9972695663 End Function Function UT_to_GST(d As Date) As Double ' ' Convert Universal Time d to Greewich Siderial Time [PDS 17] ' Dim S As Double, T As Double, T0 As Double, N As Double, H As Single H = 24 * (d - Int(d)) N = Int(d) + 2415018.5 S = N - 2451545 T = S / 36525 T0 = 6.697374558 + 2400.051336 * T + 0.000025862 * T * T T0 = T0 + H * 1.002737909 T0 = Range0_24(T0) UT_to_GST = T0 End Function Trigonometric Subroutines Function ACos(x As Double) As Double ' ' returns arccosine in radians of x ' ACos = Atn(-x / Sqr(1 - x * x)) + 2 * Atn(1) End Function Function ASin(x As Double) As Double ' ' returns arcsine in radians of x ' ASin = Atn(x / Sqr(1 - x * x)) End Function Function ATan2(x As Double, y As Double) As Double ' ' Returns arctangent (y/x) with results from 0 to 360; also tolerates x = 0 ' Dim Pi As Double Pi = 3.14159265359 If x = 0 Then If y > 0 Then ATan2 = Pi / 2 Else ATan2 = -Pi / 2 End If Else ATan2 = Atn(y / x) If x < 0 Then If y >= 0 Then ATan2 = ATan2 + Pi Else ATan2 = ATan2 - Pi End If End If End If If ATan2 < 0 Then ATan2 = ATan2 + 2 * Pi End Function Function Range0_24(a As Double) As Double Range0_24 = a While Range0_24 > 24 Range0_24 = Range0_24 - 24 Wend While Range0_24 < 0 Range0_24 = Range0_24 + 24 Wend End Function 16 Hill and Braun Geolocation by Light Level 7 1