Planetary Science

Planetary Science

The study of planets, moons, and systems

This obviously covers a lot of topics (planet and solar system formation, moon formation, meteors, geology...), but we are only focused on one: clestial orbits, i.e. the orbits of the moons and planets. To code it properly you need a pretty good understanding of how everything works as far as orbits go. Now to the point of this article:




Rainmeter

I accidentally got heavily involved in Rainmeter, a desktop customization tool that lets you display custom skins like weather forecasts, CPU monitors, RSS feeds and so on. You can choose from thousands of skins made by other users, or you can make your own (here are some examples of ones I've made, including the one this page is refering to). After using it for a couple years I thought, if it can display weather forecasts and more specifically sunset/rise times, then why not planet set and rise times. I figured there were probably a few equations I could put in a Rainmeter skin, and have it output the rise and set time of a particular planet. As it turns out, it's much more complicated than that.


Overview

Each planet has several characteristics that determine its orbit. These are called orbital elements:

  • N = longitude of the ascending node
  • i = inclination to the ecliptic
  • w = argument of perihelion
  • a = semi-major axis
  • e = eccentricity
  • M = mean anomaly

    • The Sun is a somewhat trivial example so we'll use Mars:

    • N = 49.5574 + 2.11081E-5 * d
    • i = 1.8497 - 1.78E-8 * d
    • w = 286.5016 + 2.92961E-5 * d
    • a = 1.523688
    • e = 0.093405 + 2.516E-9 * d
    • M = 18.6021 + 0.5240207766 * d

    The 'd' in these equations is the decimal day number since 2000 Jan 0.0 UT. This is our time scale. The equation for 'd' follows (y=year, m=month, D=date, UT=UT in hours+decimals):

    • d = 367*y - 7 * ( y + (m+9)/12 ) / 4 + 275*m/9 + D - 730530

    Note: these are ALL integer divisions.

    • d = d + UT/24.0

    Now that you've calculated d you can get the orbital elements and the related orbital elements:

    • w1 = N + w = longitude of perihelion
    • L = M + w1 = mean longitude
    • q = a*(1-e) = perihelion distance
    • Q = a*(1+e) = aphelion distance
    • P = a ^ 1.5 = orbital period (years if a is in AU, astronomical units)
    • T = Epoch_of_M - (M(deg)/360_deg) / P = time of perihelion
    • v = true anomaly (angle between position and perihelion)
    • E = eccentric anomaly

    Most of these calculations are straight forward except for true anomaly. Orbits aren't quite perfect circles so you can't use just cosine and sine to calculate the position along the circle. True anomaly will get you the location along a Keplerian orbit using the eccentric and mean anomaly. Once you have the true anomaly you can calculate the position in 3-dimensional space.
    These first calculations will get you the rise and set times as if you were standing on the Sun (heliocentric). You have to shift the point of reference to Earth (geocentric). You do this simply by multiplying the y and z components by cos(ecliptic) and sin(ecliptic) respectively. From this you can easily compute the Right Ascension and Declination using the atan2 function. These are known as the equatorial coordinates (much prefered over azimuthal coordinates because of tracking objects).

    Note: these equations are written out and explained in the code section.

    Code

    Summary of LUA

    Input:

  • Planet
  • The value you want returned (right ascension, declination, rise time, set time)
  • Your latitude
  • Your longitude
    • The rest of the inputs are for the date you want the information for (past, present, future)
  • Year
  • Month
  • Day
  • Hour
  • Minute

    • Initializations and Time Scale

    We start out with the header. This section is specific to Rainmeter; it's how Rainmeter sends the input: 

    PROPERTIES =
{
  planet=0;
  returnVal=0;
  latitude=0;
  longitude=0;
  Year=0;
  Month=0;
  Day=0;
  Hour=0;
  Minute=0;
}

function Initialize()
 pl = tonumber(PROPERTIES.planet)
 rtn = tonumber(PROPERTIES.returnVal)
 lat = tonumber(PROPERTIES.latitude)
 lon = tonumber(PROPERTIES.longitude)
 pY = tonumber(PROPERTIES.Year)
 pM = tonumber(PROPERTIES.Month)
 pD = tonumber(PROPERTIES.Day)
 pH = tonumber(PROPERTIES.Hour)
 pM2 = tonumber(PROPERTIES.Minute)
 pS = 30  --seconds
end
    The beginning of function Update() converts the input to a date in seconds because that's what LUA uses to get the current time. If there is no time input then we assume you want to use the current time. Then we convert the date in secs back to a normal date. We calculate the decimal time, and then the Julian Date and century for J2000. 
    function Update()
 --Date
 dateinsec = os.time{year=pY, month=pM, day=pD, hour=pH-4, min=pM2, sec=pS}
 if pY>0 then
  datetable=os.date("!*t", dateinsec)
 else
  datetable=os.date("!*t")
 end
 y=datetable.year
 m=datetable.month
 d=datetable.day
 h=datetable.hour
 m2=datetable.min
 s=datetable.sec
 --Calculate Julian Date and century
 h=h+m2/60+s/3600
 JD=367*y-math.floor((7*math.floor((y+math.floor((m+9)/12))))/4)+math.floor(275*m/9)+d-730531.5+h/24
 cy=JD/36525

    Orbital Elements

    Now we will intialize the orbital elements for the specific planet we want. I'm not going to show all the planets (repetitive), just the moon since it's slightly different and the first planet.
    The first if statement just calls another function to get the moon elements, ignore it for now (it's a work in progress, too many variables to get an accurate position). The second if statement initializes the orbital elements for Mercury, passes them to the RA/DEC function to get the planets current position around the Sun. It then calls a function to get the rise and set time for the planet. The numbers will vary for each planet. 

     if pl==0 then   --Moon
  RA, DEC = moon(y, m, d, h, m2, s) --Moon orbital elements
  RISE, SET = riseSet(RA+90, DEC, lat, lon, y, m, d, h, m2, s)
 elseif pl==1 then     --Mercury orbital elements
  A = 0.38709893 + 0.00000066*cy
  e = 0.20563069 + 0.00002527*cy
  i = math.rad( 7.00487  -  23.51*cy/3600)
  O = math.rad(48.33167  - 446.30*cy/3600)
  w = math.rad(77.45645  + 573.57*cy/3600)
  L = math.rad(252.25084 + 538101628.29*cy/3600)%(2*math.pi)
  RA, DEC = compPos(cy, A, e, i, O, w, L)
  RISE, SET = riseSet(RA, DEC, lat, lon, y, m, d, h, m2, s)
    Now for the last part of this first section. This is for Rainmeter, it basically returns either the RA, DEC, RISE, or SET for the particular planet depending on what Rainmeter wants. Unfortunately this means the program has to calculate the RA, dec, and rise and set times every time even though Rainmeter can only use one piece of information at a time. This is a limitation with Rainmeter. 
     if rtn==0 then --this is for Rainmeter, this script returns RA, DEC, RISE, or SET based on what Rainmeter wants
  return RA
 elseif rtn==1 then
  return DEC
 elseif rtn==2 then
  return RISE
 elseif rtn==3 then
  return SET
 end
end
    Now to the heart of the program. Everything else before was just setting up constants.

    Compute Position Offset

    To get meaningful rise and set times we must first calculate a postion transformation from the Sun to Earth. Otherwise our rise and set times will be as if we are standing on the Sun.
    The first set of variables are the orbital elements of Earth (sun and earth can be interchanged here because if you think about it, from Earth's perspective the Sun orbits Earth, from the Sun's perspective Earth orbits the Sun).
    Then comes the anomaly. Orbits are not perfect circles, they are elliptic. The eccentric anomaly (E) is computed from the mean anomaly (M) and from the eccentricity (e). The while statement is for error correction. Then we can calculate the true anomaly (ve) and the distance from the Sun (re). We convert these values into rectangular coordinates. 

    function earthElements(cent)
 --Earth/Sun
 --Orbital Elements
 ae = 1.00000011 - 0.00000005*cent
 ee = 0.01671022 - 0.00003804*cent
 ie = math.rad((0.00005 -    46.94*cent/3600))
 oe = math.rad((-11.26064 - 18228.25*cent/3600))
 pe = math.rad((102.94719+1198.28*cent/3600))
 le = math.rad((100.46435 + 129597740.63*cent/3600))%(2*math.pi)

 me = (le-pe)%(2*math.pi)
 --True Anomaly
 E  = me + ee*math.sin(me)*(1.0 + ee*math.cos(me))
 E1=0
 while math.abs(E - E1)>0.0000000000001 do
  E1 = E
  E = E1 - (E1 - ee*math.sin(E1) - me)/(1 - ee*math.cos(E1))
 end
 ve = (2*math.atan(math.sqrt((1 + ee)/(1 - ee))*math.tan(0.5*E)))%(2*math.pi)
 re = ae*(1 - ee*ee)/(1 + ee*math.cos(ve))
 xe = re*math.cos(ve + pe)
 ye = re*math.sin(ve + pe)
 ze = 0.0;

 return xe, ye, ze
end

    Compute RA and DEC

    This is the true heart of the program. It computes the absolute position of a planet around the Sun. First we call the earthElements function which returns the postion offset for calculating rise and set times. We perform the same steps before to calculate the distance from the Sun and the true anomaly (rp, vp).
    Then we want to convert that into a position in 3-dimensional heliocentric space (xh, yh, zh).
    Next perform a spacial transformation from the Sun's position (heliocentric) to Earth's position (geocentric) (xg, yg, zg).
    xeq, yeq, and zeq takes into account the tilt of Earth's axis.
    Lastly we get to compute what we came for. RA and DEC using the atan2 function. 

    function compPos(cy, ap, ep, ip, op, pp, lp)
 xe, ye, ze = earthElements(cy)

 mp = (lp - pp)%(2*math.pi)
 --True Anomaly
 E  = mp + ep*math.sin(mp)*(1.0 + ep*math.cos(mp))
 E1=0
 while math.abs(E - E1)>0.0000000001 do
  E1 = E
  E = E1 - (E1 - ep*math.sin(E1) - mp)/(1 - ep*math.cos(E1))
 end
 vp = (2*math.atan(math.sqrt((1 + ep)/(1 - ep))*math.tan(0.5*E)))%(2*math.pi)

 rp = ap*(1 - ep*ep)/(1 + ep*math.cos(vp))

 xh = rp*(math.cos(op)*math.cos(vp + pp - op) - math.sin(op)*math.sin(vp + pp - op)*math.cos(ip))
 yh = rp*(math.sin(op)*math.cos(vp + pp - op) + math.cos(op)*math.sin(vp + pp - op)*math.cos(ip))
 zh = rp*(math.sin(vp + pp - op)*math.sin(ip))

 xg = xh - xe
 yg = yh - ye
 zg = zh - ze

 ecl = math.rad(23.439281)
 xeq = xg
 yeq = yg*math.cos(ecl) - zg*math.sin(ecl)
 zeq = yg*math.sin(ecl) + zg*math.cos(ecl)

 ra=math.deg((math.atan2(yeq, xeq))%(2*math.pi))
 dec=math.deg(math.atan(zeq/math.sqrt(xeq*xeq + yeq*yeq)))
 return ra, dec
end

    The Rise and Set Times

    Start with converting current to time UT and then recalculate the Julian Date for J2000. Then we use a couple of the orbital elements of the Sun (w, M) to calculate the mean longitude (L).
    Next we want the Sidereal time (GMST0) which is Greenwich Midnight. From this we use our longitude to calculate our local Sidereal time (LST) and from LST we can get the Local Hour Angle (LHA) or the angle the Earth has turned since the Sun was last in the south.
    LHA2 is the Sun's Local Hour Angle at sunrise/sunset. You will notice a "-0.833" in the equation. This is how many degrees the center of the Sun is below the horizon (with atmospheric refraction).

    • 0 degrees: Center of Sun's disk touches a mathematical horizon
    • -0.25 degrees: Sun's upper limb touches a mathematical horizon
    • -0.583 degrees: Center of Sun's disk touches the horizon; atmospheric refraction accounted for
    • -0.833 degrees: Sun's upper limb touches the horizon; atmospheric refraction accounted for
    • -6 degrees: Civil twilight (one can no longer read outside without artificial illumination)
    • -12 degrees: Nautical twilight (navigation using a sea horizon no longer possible)
    • -15 degrees: Amateur astronomical twilight (the sky is dark enough for most astronomical observations)
    • -18 degrees: Astronomical twilight (the sky is completely dark)

    The eqution of time will be explained in the next section. The offset section just accounts for you timezone and daylight savings time.
    To get sunrise you take UT_Sun_in_south - LHA2 and for sunset UT_Sun_in_south + LHA2 
    function riseSet(RA, DEC, lat, lon, y, m, d, h, m2, s)

 UT = h+m2/60
 JD = 367*y-math.floor((7*math.floor((y+math.floor((m+9)/12))))/4)+math.floor(275*m/9)+d-730530 + UT/24.0

 w = 282.9404 + 4.70935E-5 * JD
 M = 356.0470 + 0.9856002585 * JD

 L = (w + M)%360

 GMST0 = ((L + 180)%360)/15

 LST = GMST0 + UT + lon/15
 LHA = LST - RA/15

 UT_Sun_in_south = (RA/15 - GMST0 - lon/15)%24

 LHA2 = math.acos((math.sin(math.rad(-0.833)) - math.sin(math.rad(lat))*math.sin(math.rad(DEC)))/(math.cos(math.rad(lat))*math.cos(math.rad(DEC))))
 --equation of time
 JD2 = juliandate(y, m, d, h, m2, s)
 T = (JD2-2451545)/36525
 eqnTime = equationOfTime(T)

 --Offset
 currenttime = os.time()
 utc = os.date("!*t", currenttime)
 lcl = os.date("*t",  currenttime)
 offset = (utc.hour*60 + utc.min)  - (lcl.hour*60 + lcl.min)

 rise = -(math.deg(LHA2)/15)+UT_Sun_in_south+eqnTime/60-(offset/60)
 set = (math.deg(LHA2)/15)+UT_Sun_in_south+eqnTime/60-(offset/60)
 return rise, set
end

    Calculating the Equation of Time

    We start out with a slight variation of the Julian Date. This is because we need to calculate it using Greenwich time at 0 hours (normal Julian day begins at noon). 

    function juliandate(y, m, d, h, m2, s)
 YearDur = 365.25
 if (m<=2) then
  y=y-1
  m=m+12
 end
 A = math.floor(YearDur*(y+4716))
 B = math.floor(30.6001*(m+1))
 C = 2
 D = math.floor(y/100)
 E = math.floor(math.floor(y/100)*0.25)
 F = d-1524.5
 G = (h+(m2/60)+s/3600)/24
 jd = A+B+C-D+E+F+G
 return jd
end
    The Sun doesn't stay in the same position in the sky throughout the year. It's lower in the winter and higher in the summer. This is where the equation of time comes in (it's not so much of an equation as it is a correction). 
    function equationOfTime(t)
 epsilon = calcObliquityCorrection(t)
 LO = (280.46646 + t * (36000.76983 + t*(0.0003032)))%360
 e = 0.016708634 - t * (0.000042037 + 0.0000001267 * t)
 m = 357.52911 + t * (35999.05029 - 0.0001537 * t)
 y = (math.tan(math.rad(epsilon)/2.0))^2
 sin2l0 = math.sin(2*math.rad(LO))
 sinm = math.sin(math.rad(m))
 cos2l0 = math.cos(2*math.rad(LO))
 sin4l0 = math.sin(4*math.rad(LO))
 sin2m = math.sin(2*math.rad(m))
 ETime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m
 return math.deg(ETime)*4
end
    Now that we've calculate the RA, Dec, sunrise, and sunset we are ready to output it into something that makes a little more since (visually): Rainmeter

    Rainmeter Code

    Right Ascension is the angle measured eastward along the celestial equator from the vernal equinox. Declination is the angle measure from north or south of the celestial equator. These two direction coordinates are used in the equatorial coordinate system. The Azimuth coordinate system is what most people are used to: longitude and latitude to put it simply. We are going to ignore Azimuth because it doesn't make much sense to use in astronomy. Because the Earth is tilted, stars move in circular paths, RA and DEC take this into consideration so stars will "change" RA throughout the night but DEC will stay the same. This makes tracking objects easy because you only have to move your telescope along one axis.
    We will start with a basic Star and Planet Locator and ignore the planets and Sun for now. We want it in real time so we need to find out how many degrees the stars rotate in one day. You may believe the Earth rotates 360 degrees in one day and you would be wrong. It actually rotates 360.985608 degrees in one day (not getting in to this part). Here is a formula for calculating how much the stars rotate in one minute:

    Formula=-((((EarthDay+253)%366*360.985608)+(EarthHourFix*15.041067)+(EarthMinute*0.250687))%360)

    I arbitrarily chose the 253rd day of 2012 as a starting point. EarthDay is the number of days since then. 15.041067 is how much it rotates in an hour and 0.250687 degrees is how much it rotates in one minute. This is all you need if you just want to know what stars are out at a certain time of the day.
    Now we must convert RA and DEC into X and Y coordinates for Rainmeter to plot the positions of our planets on to our Star and Planet Locator. This is achieved via this formula:

    X = 3.9*(90-YposSun)*cos((XposSun+((-EarthHour)*15.041)+((-0-EarthMinute)*0.250687)+(-12-EarthDay)*0.98)*PI/180)+725
    Y = 3.9*(90-YposSun)*sin((XposSun+((-EarthHour)*15.041)+((-0-EarthMinute)*0.250687)+(-12-EarthDay)*0.98)*PI/180)+725

    I'm going to break this formula down, but first the variables: YposSun=DEC, XposSun=RA, EarthHour, EarthMinute=time since midnight, EarthDay=day of year

  • 3.9: First we need to know how many pixels are in a degree to scale the X and Y positions properly. For my resolution there are 3.9 pixels in one degree.
  • (90-YposSun): Next we notice that the ecliptic (path of the Sun) doesn't stay in one position in the sky throughout the year. It is offset by roughly 90 pixels (90/3.9=23degrees, the tilt of the Earth)
  • The rest of the formula is just how far the Sun has moved since 12:01am for the current day.
  • (-12-EarthDay)*0.98=takes in to account the equation of time (90-YposSun is just for the initial position, this calculates the position throughout the year)

    • And We're Done!

    That wasn't so hard was it.



    No comments:

    Post a Comment

    Subscribe to: Comments (Atom)