class SunRiseSet

Calculates the sun rise and sunset times, with civil, naval and astronomical twilight values. Not sure of the origin of the code. I have seen a fortran version www.srrb.noaa.gov/highlights/sunrise/program.txt a .pl www.mso.anu.edu.au/~brian/grbs/astrosubs.pl and .vb versions too. All had the same comments, so are of a common origin.

Constants

ASTRO_TWILIGHT

In degrees from the Zenith. Represents the Time when the sun is not interfering with viewing distant stars.

CIVIL_TWILIGHT

In degrees from the Zenith. Represents the Time when we turn car lights on and off

LATITUDE_DEFAULT

because I live here

LONGITUDE_DEFAULT

In degrees from the Zenith. Represents the Time when we can see the first light (dawn)

SUN_RISE_SET

0.833 is allowing for the bending in the atmosphere.

VERSION

Attributes

astroTwilightEnd[R]

@return [DateTime] end of astronomical twilight (sky is now fully dark)

astroTwilightStart[R]

@return [DateTime] Naval Twilight begins (Sun is begining to lighten the sky )

civilTwilightEnd[R]

@return [DateTime] End of Civil Twilight

civilTwilightStart[R]

@return [DateTime] Civil Twilight begins

navalTwilightEnd[R]

@return [DateTime] end of naval twilight

navalTwilightStart[R]

@return [DateTime] Naval Twilight begins

solNoon[R]

@return [DateTime] Sun is at the midpoint for today (varies throughout the year)

sunrise[R]

@return [DateTime] Sun rise

sunset[R]

@return [DateTime] Sun set

Public Class Methods

new(datetime, latitude=LATITUDE_DEFAULT, longitude=LONGITUDE_DEFAULT) click to toggle source

@return [SunRiseSet] Constructor for any datetime and location @param [DateTime, jd, offset] datetime @param [Float] latitude @param [Float] longitude

# File lib/sunriseset.rb, line 51
def initialize(datetime, latitude=LATITUDE_DEFAULT, longitude=LONGITUDE_DEFAULT)
  @latitude, @longitude = latitude, longitude
  @julian_date = DateTime.jd(datetime.jd.to_f)
  @julian_day = @julian_date.jd.to_f #Shorthand for later use, where we want this value as a float.
  @zone = datetime.offset #datetime.utc_offset/86400 #time zone offset + daylight savings as a fraction of a day
  calcSun
end
now(latitude=LATITUDE_DEFAULT, longitude=LONGITUDE_DEFAULT) click to toggle source

@return [SunRiseSet] Constructor for date == today, at location specified @param [Float] latitude @param [Float] longitude

# File lib/sunriseset.rb, line 69
def self.now(latitude=LATITUDE_DEFAULT, longitude=LONGITUDE_DEFAULT)
  self.new(DateTime.now, latitude, longitude)
end
today(latitude=LATITUDE_DEFAULT, longitude=LONGITUDE_DEFAULT) click to toggle source

@return [SunRiseSet] Constructor for date == today, at location specified @param [Float] latitude @param [Float] longitude

# File lib/sunriseset.rb, line 62
def self.today(latitude=LATITUDE_DEFAULT, longitude=LONGITUDE_DEFAULT)
  self.new(DateTime.now, latitude, longitude)
end

Public Instance Methods

time_format(t) click to toggle source
# File lib/sunriseset.rb, line 73
def time_format(t)
  if t == nil
    "Not Found"
  else
    t.strftime('%H:%M:%S %d-%m')
  end
end
to_s() click to toggle source

@return [String] dumps key attributes as a multiline string

# File lib/sunriseset.rb, line 82
def to_s
  "Astro Twilight #{time_format(@astroTwilightStart)}\n" +
  "Naval Twilight #{time_format(@navalTwilightStart)}\n" +
  "Civil Twilight #{time_format(@civilTwilightStart)}\n" +
  "Sun Rises #{time_format(@sunrise)}\n" +
  "Solar noon  #{time_format(@solNoon)}\n" +
  "Sun Sets #{time_format(@sunset)}\n" +
  "End of Civil Twilight  #{time_format(@civilTwilightEnd)}\n" +
  "Naval Twilight #{time_format(@navalTwilightEnd)}\n" +
  "Astro Twilight #{time_format(@astroTwilightEnd)}\n" 
end
version() click to toggle source

@return [String] the constant VERSION

# File lib/sunriseset.rb, line 95
def version
  VERSION
end

Private Instance Methods

calcEccentricityEarthOrbit(t) click to toggle source

calculate the eccentricity of earth's orbit @param [Float] t number of Julian centuries since J2000.0 @return [Float] the unitless eccentricity

# File lib/sunriseset.rb, line 224
def calcEccentricityEarthOrbit(t)
  0.016708634 - t * (0.000042037 + 0.0000001267 * t)  # unitless
end
calcEquationOfTime(t) click to toggle source

calculate the difference between true solar time and mean solar time @param [Float] t number of Julian centuries since J2000.0 @return [Float] equation of time in minutes of time

# File lib/sunriseset.rb, line 316
def calcEquationOfTime(t)
  epsilon = calcObliquityCorrection(t)
  l0 = calcGeomMeanLongSun(t)
  e = calcEccentricityEarthOrbit(t)
  m = calcGeomMeanAnomalySun(t)

  y = Math.tan(degToRad(epsilon)/2.0)
  y *= y

  sin2l0 = Math.sin(2.0 * degToRad(l0))
  sinm   = Math.sin(degToRad(m))
  cos2l0 = Math.cos(2.0 * degToRad(l0))
  sin4l0 = Math.sin(4.0 * degToRad(l0))
  sin2m  = Math.sin(2.0 * degToRad(m))

  radToDeg(y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0 -
              0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m)*4.0  # in minutes of time
end
calcGeomMeanAnomalySun(t) click to toggle source

Calculate the Geometric Mean Anomaly of the Sun @param [Float] t number of Julian centuries since J2000.0 @return [Float] the Geometric Mean Anomaly of the Sun in degrees

# File lib/sunriseset.rb, line 217
def calcGeomMeanAnomalySun(t)
  357.52911 + t * (35999.05029 - 0.0001537 * t) #in degrees
end
calcGeomMeanLongSun(t) click to toggle source

calculate the Geometric Mean Longitude of the Sun @param [Float] t number of Julian centuries since J2000.0 @return [Float] the Geometric Mean Longitude of the Sun in degrees

# File lib/sunriseset.rb, line 203
def calcGeomMeanLongSun(t)
  l0 = 280.46646 + t * (36000.76983 + 0.0003032 * t)
  while(l0 > 360.0)
    l0 -= 360.0
  end
  while(l0 < 0.0)
    l0 += 360.0
  end
  return l0;    # in degrees
end
calcHourAngleSunrise(solarDec, angle=SUN_RISE_SET) click to toggle source

calculate the hour angle of the sun at sunrise for the latitude @param [Float] solarDec declination angle of sun in degrees @param [SUN_RISE_SET,CIVIL_TWILIGHT,NAVAL_TWILIGHT,ASTRO_TWILIGHT] angle @return [Float] hour angle of sunrise in radians 0.833 is an approximation of the reflaction caused by the atmosphere

# File lib/sunriseset.rb, line 340
def calcHourAngleSunrise(solarDec, angle=SUN_RISE_SET)
  latRad = degToRad(@latitude) 
  sdRad  = degToRad(solarDec)
  #puts "latRad = #{radToDeg(latRad)}, sdRad = #{radToDeg(sdRad)}, angle = #{angle}"

  #ha_arg = Math.cos(degToRad(angle + 0.833))/(Math.cos(latRad)*Math.cos(sdRad))-Math.tan(latRad) * Math.tan(sdRad)
  ha_arg = Math.cos(degToRad(angle))/
          (Math.cos(latRad)*Math.cos(sdRad)) - Math.tan(latRad) * Math.tan(sdRad)
  Math.acos(ha_arg)  # in radians
end
calcHourAngleSunset( solarDec, angle=SUN_RISE_SET) click to toggle source

calculate the hour angle of the sun at sunset for the @param [Float] solarDec declination angle of sun in degrees @param [SUN_RISE_SET,CIVIL_TWILIGHT,NAVAL_TWILIGHT,ASTRO_TWILIGHT] angle @return [Float] hour angle of sunset in radians

# File lib/sunriseset.rb, line 355
def calcHourAngleSunset( solarDec, angle=SUN_RISE_SET)
  -calcHourAngleSunrise(solarDec,angle)   # in radians
end
calcJDFromJulianCent(t) click to toggle source

convert centuries since J2000.0 to Julian Day. @param [Float] t number of Julian centuries since J2000.0 @return [Float] the Julian Day corresponding to the t value

# File lib/sunriseset.rb, line 196
def calcJDFromJulianCent(t)
  t * 36525.0 + 2451545.0
end
calcMeanObliquityOfEcliptic(t) click to toggle source

calculate the mean obliquity of the ecliptic @param [Float] t number of Julian centuries since J2000.0 @return [Float] mean obliquity in degrees

# File lib/sunriseset.rb, line 278
def calcMeanObliquityOfEcliptic(t)
  seconds = 21.448 - t*(46.8150 + t*(0.00059 - t*(0.001813)))
  23.0 + (26.0 + (seconds/60.0))/60.0   # in degrees
end
calcObliquityCorrection(t) click to toggle source

calculate the corrected obliquity of the ecliptic @param [Float] t number of Julian centuries since J2000.0 @return [Float] corrected obliquity in degrees

# File lib/sunriseset.rb, line 286
def calcObliquityCorrection(t)
  e0 = calcMeanObliquityOfEcliptic(t)
  omega = 125.04 - 1934.136 * t
  e0 + 0.00256 * Math.cos(degToRad(omega))    # in degrees
end
calcSolNoonUTC(t) click to toggle source

calculate the Universal Coordinated Time (UTC) of solar noon for the given day at the given location on earth @param [Float] t number of Julian centuries since J2000.0 @return [Float] time in minutes from zero Z

# File lib/sunriseset.rb, line 404
def calcSolNoonUTC(t)
  # First pass uses approximate solar noon to calculate eqtime
  tnoon = calcTimeJulianCent(calcJDFromJulianCent(t) - @longitude/360.0)
  eqTime = calcEquationOfTime(tnoon)
  solNoonUTC = 720 - (@longitude * 4) - eqTime; # min

  newt = calcTimeJulianCent(calcJDFromJulianCent(t) -0.5 + solNoonUTC/1440.0)

  eqTime = calcEquationOfTime(newt)
  solNoonUTC = 720 - (@longitude * 4) - eqTime; # min

  return solNoonUTC
end
calcSun() click to toggle source

calculate time of sunrise and sunset for the entered date and location. In the special cases near earth's poles, the date of nearest sunrise and set are reported. @return [nil] fills in class DateTime attributes :sunrise, :civilTwilightStart, :navalTwilightStart, :astroTwilightStart :sunset, :civilTwilightEnd, :navalTwilightEnd and :astroTwilightEnd

# File lib/sunriseset.rb, line 108
def calcSun

    # Calculate sunrise for this date
    # if no sunrise is found, set flag nosunrise

    @sunrise = calcSunriseUTC(@julian_day)
    @civilTwilightStart = calcSunriseUTC( @julian_day, CIVIL_TWILIGHT)
    @navalTwilightStart = calcSunriseUTC( @julian_day, NAVAL_TWILIGHT)
    @astroTwilightStart = calcSunriseUTC( @julian_day, ASTRO_TWILIGHT)
    
    # Calculate sunset for this date
    # if no sunrise is found, set flag nosunset
    @sunset = calcSunsetUTC(@julian_day)
    @civilTwilightEnd = calcSunsetUTC( @julian_day, CIVIL_TWILIGHT)
    @navalTwilightEnd = calcSunsetUTC( @julian_day, NAVAL_TWILIGHT)
    @astroTwilightEnd = calcSunsetUTC( @julian_day, ASTRO_TWILIGHT)

    # Calculate solar noon for this date
    t = calcTimeJulianCent( @julian_day )
    
    @solNoon = to_datetime(@julian_day, calcSolNoonUTC(t))

    # No sunrise or sunset found for today
    doy = @julian_date.yday
    if(@sunrise == nil)
      if ( ((@latitude > 66.4) && (doy > 79) && (doy < 267)) ||
         ((@latitude < -66.4) && ((doy < 83) || (doy > 263))) )
         # if Northern hemisphere and spring or summer, OR
         # if Southern hemisphere and fall or winter, use
         # previous sunrise and next sunset
        newjd = findRecentSunrise(@julian_day, @latitude, @longitude)
        @sunrise = calcSunriseUTC(newjd) + @zone
      elsif ( ((@latitude > 66.4) && ((doy < 83) || (doy > 263))) ||
        ((@latitude < -66.4) && (doy > 79) && (doy < 267)) )
        # if Northern hemisphere and fall or winter, OR
        # if Southern hemisphere and spring or summer, use
        # next sunrise and previous sunset
        newjd = findNextSunrise(@julian_day, @latitude, @longitude)
        @sunrise = calcSunriseUTC(newjd) + @zone
      else
        raise "Cannot Find Sunrise!"
      end

    end

    if(@sunset == nil)
      if ( ((@latitude > 66.4) && (doy > 79) && (doy < 267)) ||
        ((@latitude < -66.4) && ((doy < 83) || (doy > 263))) )
        # if Northern hemisphere and spring or summer, OR
        # if Southern hemisphere and fall or winter, use
        # previous sunrise and next sunset
        newjd = findNextSunset(@julian_day, @latitude, @longitude)
        @sunset = calcSunsetUTC(newjd) + @zone
      elsif ( ((@latitude > 66.4) && ((doy < 83) || (doy > 263))) ||
        ((@latitude < -66.4) && (doy > 79) && (doy < 267)) )
        # if Northern hemisphere and fall or winter, OR
        # if Southern hemisphere and spring or summer, use
        # next sunrise and last sunset
        newjd = findRecentSunset(@julian_day, @latitude, @longitude)
        @sunset = calcSunsetUTC(newjd) + @zone
      else
        raise "Cannot Find Sunset!"
      end
    end
end
calcSunApparentLong(t) click to toggle source

calculate the apparent longitude of the sun @param [Float] t number of Julian centuries since J2000.0 @return [Float] sun's apparent longitude in degrees

# File lib/sunriseset.rb, line 269
def calcSunApparentLong(t)
  o = calcSunTrueLong(t)
  omega = 125.04 - 1934.136 * t
  o - 0.00569 - 0.00478 * Math.sin(degToRad(omega)) # in degrees
end
calcSunDeclination(t) click to toggle source

calculate the declination of the sun @param [Float] t number of Julian centuries since J2000.0 @return [Float] sun's declination in degrees

# File lib/sunriseset.rb, line 306
def calcSunDeclination(t)
  e = calcObliquityCorrection(t)
  lambda = calcSunApparentLong(t)
  sint = Math.sin(degToRad(e)) * Math.sin(degToRad(lambda))
  radToDeg(Math.asin(sint)) # in degrees
end
calcSunEqOfCenter(t) click to toggle source

calculate the equation of center for the sun @param [Float] t number of Julian centuries since J2000.0 @return [Float] in degrees

# File lib/sunriseset.rb, line 231
def calcSunEqOfCenter(t)

  m = calcGeomMeanAnomalySun(t)

  mrad = degToRad(m)
  sinm = Math.sin(mrad)
  sin2m = Math.sin(mrad+mrad)
  sin3m = Math.sin(mrad+mrad+mrad)

  sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) + sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289 # in degrees
end
calcSunRadVector(t) click to toggle source

calculate the distance to the sun in AU @param [Float] t number of Julian centuries since J2000.0 @return [Float] sun radius vector in AUs

# File lib/sunriseset.rb, line 260
def calcSunRadVector(t)
  v = calcSunTrueAnomaly(t)
  e = calcEccentricityEarthOrbit(t)
  (1.000001018 * (1.0 - e * e)) / (1.0 + e * Math.cos(degToRad(v))) # in AUs
end
calcSunRtAscension(t) click to toggle source

calculate the right ascension of the sun @param [Float] t number of Julian centuries since J2000.0 @return [Float] sun's right ascension in degrees

# File lib/sunriseset.rb, line 295
def calcSunRtAscension(t)
  e = calcObliquityCorrection(t)
  lambda = calcSunApparentLong(t)
  tananum = (Math.cos(degToRad(e)) * Math.sin(degToRad(lambda)))
  tanadenom = (Math.cos(degToRad(lambda)))
  radToDeg(Math.atan2(tananum, tanadenom))  # in degrees
end
calcSunTrueAnomaly(t) click to toggle source

calculate the true anamoly of the sun @param [Float] t number of Julian centuries since J2000.0 @return [Float] sun's true anamoly in degrees

# File lib/sunriseset.rb, line 253
def calcSunTrueAnomaly(t)
  calcGeomMeanAnomalySun(t) + calcSunEqOfCenter(t)  # in degrees
end
calcSunTrueLong(t) click to toggle source

calculate the true longitude of the sun @param [Float] t number of Julian centuries since J2000.0 @return [Float] sun's true longitude in degrees

# File lib/sunriseset.rb, line 246
def calcSunTrueLong(t)
  calcGeomMeanLongSun(t) + calcSunEqOfCenter(t) # in degrees
end
calcSunriseUTC( julian_day, angle=SUN_RISE_SET) click to toggle source

calculate the Universal Coordinated Time (UTC) of sunrise for the given day at the given location on earth @param [Float] julian_day @param [SUN_RISE_SET,CIVIL_TWILIGHT,NAVAL_TWILIGHT,ASTRO_TWILIGHT] angle @return [DateTime] Date and Time of event

# File lib/sunriseset.rb, line 364
def calcSunriseUTC( julian_day, angle=SUN_RISE_SET)
  begin
    t = calcTimeJulianCent( julian_day )

    # *** Find the time of solar noon at the location, and use
        #     that declination. This is better than start of the
        #     Julian day

    noonmin = calcSolNoonUTC(t)
    tnoon = calcTimeJulianCent( julian_day+noonmin/1440.0)

    # *** First pass to approximate sunrise (using solar noon)

    eqTime = calcEquationOfTime(tnoon)
    solarDec = calcSunDeclination(tnoon)
    hourAngle = calcHourAngleSunrise(solarDec,angle)
    delta =  -@longitude - radToDeg(hourAngle)
    timeDiff = 4 * delta; # in minutes of time
    timeUTC = 720 + timeDiff - eqTime;  # in minutes

    # *** Second pass includes fractional jday in gamma calc

    newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC/1440.0)
    eqTime = calcEquationOfTime(newt)
    solarDec = calcSunDeclination(newt)
    hourAngle = calcHourAngleSunrise(solarDec,angle)
    delta = -@longitude - radToDeg(hourAngle)
    timeDiff = 4 * delta
    timeUTC = 720 + timeDiff - eqTime; # in minutes

    to_datetime(julian_day,timeUTC)
  rescue Math::DomainError => error
    return nil #didn't find a Sunrise today. Will be raised by
  end
end
calcSunsetUTC(julian_day, angle=SUN_RISE_SET) click to toggle source

calculate the Universal Coordinated Time (UTC) of sunset for the given day at the given location on earth @param [Float] julian_day @param [SUN_RISE_SET,CIVIL_TWILIGHT,NAVAL_TWILIGHT,ASTRO_TWILIGHT] angle @return [DateTime] Date and Time of event

# File lib/sunriseset.rb, line 423
def calcSunsetUTC(julian_day, angle=SUN_RISE_SET)
  begin
    t = calcTimeJulianCent(julian_day)

    # *** Find the time of solar noon at the location, and use
        #     that declination. This is better than start of the
        #     Julian day

    noonmin = calcSolNoonUTC(t)
    tnoon = calcTimeJulianCent(julian_day+noonmin/1440.0)

    # First calculates sunrise and approx length of day

    eqTime = calcEquationOfTime(tnoon)
    solarDec = calcSunDeclination(tnoon)
    hourAngle = calcHourAngleSunset(solarDec, angle)

    delta = -@longitude - radToDeg(hourAngle)
    timeDiff = 4 * delta
    timeUTC = 720 + timeDiff - eqTime

    # first pass used to include fractional day in gamma calc

    newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC/1440.0)
    eqTime = calcEquationOfTime(newt)
    solarDec = calcSunDeclination(newt)
    hourAngle = calcHourAngleSunset(solarDec, angle)

    delta = -@longitude - radToDeg(hourAngle)
    timeDiff = 4 * delta
    timeUTC = 720 + timeDiff - eqTime; # in minutes

    to_datetime(julian_day,timeUTC)
  rescue Math::DomainError => error
    return nil # no Sunset
  end
end
calcTimeJulianCent(julian_day) click to toggle source

convert Julian Day to centuries since J2000.0. @param [Float] julian_day Julian Day @return [Float] T value corresponding to the Julian Day

# File lib/sunriseset.rb, line 189
def calcTimeJulianCent(julian_day)
  (julian_day - 2451545.0)/36525.0
end
degToRad(angleDeg) click to toggle source

@param [Float] degrees @return [Float] angle in radians

# File lib/sunriseset.rb, line 182
def degToRad(angleDeg)
  return (Math::PI * angleDeg / 180.0)
end
findNextSunrise(julianday) click to toggle source

calculate the julian day of the next sunrise starting from the given day at the given location on earth @param [Float] julianday @return [Float] julian day of the next sunrise

# File lib/sunriseset.rb, line 493
def findNextSunrise(julianday)
  time = calcSunriseUTC(julianday)
  while(!isNumber(time))
    julianday += 1.0
    time = calcSunriseUTC(julianday)
  end

  return julianday
end
findNextSunset(julianday) click to toggle source

calculate the julian day of the next sunset starting from the given day at the given location on earth @param [Float] julianday @return [Float] julian day of the next sunset

# File lib/sunriseset.rb, line 507
def findNextSunset(julianday)
  time = calcSunsetUTC(julianday)
  while(!isNumber(time))
    julianday += 1.0
    time = calcSunsetUTC(julianday)
  end

  return julianday
end
findRecentSunrise(julianday) click to toggle source

calculate the julian day of the most recent sunrise starting from the given day at the given location on earth @param [Float] julianday @return [Float] julian day of the most recent sunrise

# File lib/sunriseset.rb, line 465
def findRecentSunrise(julianday)
  time = calcSunriseUTC(julianday)
  while(!isNumber(time))
    julianday -= 1.0
    time = calcSunriseUTC(julianday)
  end

  return julianday
end
findRecentSunset(julianday) click to toggle source

calculate the julian day of the most recent sunset starting from the given day at the given location on earth @param [Float] julianday @return [Float] julian day of the most recent sunset

# File lib/sunriseset.rb, line 479
def findRecentSunset(julianday)
  time = calcSunsetUTC(julianday)
  while(!isNumber(time))
    julianday -= 1.0
    time = calcSunsetUTC(julianday)
  end

  return julianday
end
radToDeg(angleRad) click to toggle source

@param [Float] radian @return [Float] angle in degrees

# File lib/sunriseset.rb, line 176
def radToDeg(angleRad)
  return (180.0 * angleRad / Math::PI)
end
to_datetime(x,minutes) click to toggle source

convert julian day and minutes to datetime @param [Float] minutes @return [DateTime]

# File lib/sunriseset.rb, line 520
def to_datetime(x,minutes)
  jd = DateTime.jd(@julian_day)  + (minutes/1440.0)  + @zone 
end