From: "Glen F. Pankow" Date: 2010-01-20T03:33:41+09:00 Subject: Re: [QUIZ] Solar System (#227) --Boundary_(ID_GkqsRnC8mcPrCjllCkjopw) Content-type: text/plain; charset=ISO-8859-1; format=flowed Content-transfer-encoding: 7BIT Daniel Moore wrote: > > ## Solar System (#227) > > This week's quiz is to create a solar system simulator. It will a date > as input then return a list of the distances to each planet from Earth > at that date. Not a whole lot of Ruby-specific stuff here, FWIW. --Boundary_(ID_GkqsRnC8mcPrCjllCkjopw) Content-type: text/plain; name=quiz_227 Content-transfer-encoding: 7BIT Content-disposition: inline; filename=quiz_227 #! /usr/bin/ruby # # Quiz 227 -- Solar System # # Given a date, computes the (approximate) distances from the Earth to each # planet in the solar system. We'll throw in the Sun and Pluto as well. # # The formulae for this script were taken from the 'Keplerian Elements for # Approximate Positions of the Major Planets', # (http://ssd.jpl.nasa.gov/?planet_pos). We consider only dates in the range # 1800 AD - 2050 AD. # # Glen Pankow 01/16/2009 Original version. # # Note -- the nomenclature is ugly in order to match the formulae. For # example, the local variable _T has a leading underscore to keep it from # being interpreted as a Ruby global... # Planets = %w{ Sun Mercury Venus Earth Mars Jupiter Saturn Uranus Neptune Pluto } Deg_to_rad = Math::PI / 180.0 Rad_to_deg = 180.0 / Math::PI # # Table 1: Keplerian elements and their rates, with respect to the mean # ecliptic and equinox of J2000, valid for time-interval 1800 AD - 2050 AD. # keplerians = { # x-naught, x-dot 'Mercury' => [ 0.38709927, 0.00000037, 0.20563593, 0.00001906, 7.00497902, -0.00594749, 252.25032350, 149472.67411175, 77.45779628, 0.16047689, 48.33076593, -0.12534081 ], 'Venus' => [ 0.72333566, 0.00000390, 0.00677672, -0.00004107, 3.39467605, -0.00078890, 181.97909950, 58517.81538729, 131.60246718, 0.00268329, 76.67984255, -0.27769418 ], 'Earth' => [ 1.00000261, 0.00000562, 0.01671123, -0.00004392, -0.00001531, -0.01294668, 100.46457116, 35999.37244981, 102.93768193, 0.32327364, 0.0, 0.0 ], 'Mars' => [ 1.52371034, 0.00001847, # a [au, au/cty] 0.09339410, 0.00007882, # e [ , /cty] 1.84969142, -0.00813131, # I [deg, deg/cty] -4.55343205, 19140.30268499, # L [deg, deg/cty] -23.94362959, 0.44441088, # omega-bar [deg, deg/cty] 49.55953891, -0.29257343 ], # big-omega [deg, deg/cty] 'Jupiter' => [ 5.20288700, -0.00011607, 0.04838624, -0.00013253, 1.30439695, -0.00183714, 34.39644051, 3034.74612775, 14.72847983, 0.21252668, 100.47390909, 0.20469106 ], 'Saturn' => [ 9.53667594, -0.00125060, 0.05386179, -0.00050991, 2.48599187, 0.00193609, 49.95424423, 1222.49362201, 92.59887831, -0.41897216, 113.66242448, -0.28867794 ], 'Uranus' => [ 19.18916464, -0.00196176, 0.04725744, -0.00004397, 0.77263783, -0.00242939, 313.23810451, 428.48202785, 170.95427630, 0.40805281, 74.01692503, 0.04240589 ], 'Neptune' => [ 30.06992276, 0.00026291, 0.00859048, 0.00005105, 1.77004347, 0.00035372, -55.12002969, 218.45945325, 44.96476227, -0.32241464, 131.78422574, -0.00508664 ], 'Pluto' => [ 39.48211675, -0.00031596, 0.24882730, 0.00005170, 17.14001206, 0.00004818, 238.92903833, 145.20780515, 224.06891629, -0.04062942, 110.30393684, -0.01183482 ] } # 'Earth' is actually the Earth-Moon Barycenter. # # _T_eph = date_str_to_ephemeris_date(date_str) # # Return the Julian Ephemeris Date for the date string of the form # 'mm/dd/yyyy' (or simple variants thereof), assuming UTC. ArgumentError is # raised if the year field is before 1800 AD or after 2050 AD or is # not of the expected form (using very simplistic checking). # # (Well, technically, the conversion is consistent with the Julian Date as # defined by the U.S. Naval Meteorology and Oceanography Command's (NMOC) # Astronomical Applications Department's Julian Date Converter # (http://aa.usno.navy.mil/data/docs/JulianDate.php), which I'm assuming is # consistent with the NASA JPL meaning of the term Julian Ephemeris Date.) # def date_str_to_ephemeris_date(date_str) # # Parse the date string; convert it into a Time object. Raise exceptions # on any error. # raise ArgumentError, "The date string '#{date_str}' is not of the expected form" \ " ('mm/dd/yyyy')." \ unless (date_str =~ %r{^(\d+)/(\d+)/(\d\d+)$}) month, day, year = $1.to_i, $2.to_i, $3.to_i raise ArgumentError, "Invalid month field '#{$1}'." \ unless ((1 <= month) && (month <= 12)) raise ArgumentError, "Invalid day field '#{$2}'." \ unless ((1 <= day) && (day <= 31)) if (year <= 70) year += 2000 elsif (year < 100) year += 1900 end raise ArgumentError, "Dates earlier than 1800 AD are not currently supported." \ unless (year >= 1800) raise ArgumentError, "Dates later than 2050 AD are not currently supported." \ unless (year <= 2050) date = Time.gm(year, month, day) # # Convert it into the Julian Ephemeris Date. # 2451545.0 + (date - Time.gm(2000)) / 86400.0 end # # _E_deg = keplers_equation(_M_deg, e_rad) # # Solve Kepler's equation, _M_deg = _E_deg - e_rad * sin(_E_deg), for _E_deg. # def keplers_equation(_M_deg, e_rad) _M_rad = _M_deg * Deg_to_rad e_deg = e_rad * Rad_to_deg n = 0 e_n_deg = _M_deg + e_deg * Math.sin(_M_rad) loop do e_n_rad = e_n_deg * Deg_to_rad delta_M_deg = _M_deg - (e_n_deg - e_deg * Math.sin(e_n_rad)) delta_E_deg = delta_M_deg / (1.0 - e_rad * Math.cos(e_n_rad)) e_n_deg += delta_E_deg return e_n_deg if (delta_E_deg.abs <= 1e-6) n += 1 end end # # x, y, z = position_point(keplerians, _T) # # Compute the coordinates, in the J2000 ecliptic plane, x-axis aligned toward # the equinox, in AU units, of the planetary body defined by the Keplerian # fields for the number of centuries past J2000 <_T>. # def position_point(keplerians, _T) # # 1. Compute the value of each of the planet's six elements... # a_au = keplerians[0] + keplerians[1] * _T # semimajor axis e_rad = keplerians[2] + keplerians[3] * _T # eccentricity _I_deg = keplerians[4] + keplerians[5] * _T # inclination _L_deg = keplerians[6] + keplerians[7] * _T # mean longitude obar_deg = keplerians[8] + keplerians[9] * _T # lon of perihelieon _O_deg = keplerians[10] + keplerians[11] * _T # lon of ascending node _I_rad = _I_deg * Deg_to_rad _O_rad = _O_deg * Deg_to_rad # # 2. Compute the argument of perihelion and the mean anomaly... # o_deg = obar_deg - _O_deg o_rad = o_deg * Deg_to_rad _M_deg = _L_deg - obar_deg # _M_deg += b * _T * _T + c * cos(f * _T) + s * sin(f * _T) # ^^^ Jupiter-Pluto; 3000BC-3000AD # # 3. Modulus the mean anomoly... then obtain the eccentric anomoly... # _M_deg -= 360.0 while (_M_deg > 180.0) _M_deg += 360.0 while (_M_deg < -180.0) _M_rad = _M_deg * Deg_to_rad _E_deg = keplers_equation(_M_deg, e_rad) _E_rad = _E_deg * Deg_to_rad # # 4. Compute the planet's heliocentric coordinates in its orbital plane... # x_prime = a_au * (Math.cos(_E_rad) - e_rad) y_prime = a_au * Math.sqrt(1.0 - e_rad * e_rad) * Math.sin(_E_rad) # z_prime = 0.0 # # 5. Compute the coordinates in the J2000 ecliptic plane... # cos_o = Math.cos(o_rad) ; sin_o = Math.sin(o_rad) cos_O = Math.cos(_O_rad) ; sin_O = Math.sin(_O_rad) cos_I = Math.cos(_I_rad) ; sin_I = Math.sin(_I_rad) cos_o_cos_O = cos_o * cos_O cos_o_sin_O = cos_o * sin_O sin_o_cos_O = sin_o * cos_O sin_o_sin_O = sin_o * sin_O x_ecl = ( cos_o * cos_O - sin_o_sin_O * cos_I) * x_prime \ + (-sin_o * cos_O - cos_o_sin_O * cos_I) * y_prime y_ecl = ( cos_o * sin_O + sin_o_cos_O * cos_I) * x_prime \ + (-sin_o * sin_O + cos_o_cos_O * cos_I) * y_prime z_ecl = sin_o * sin_I * x_prime \ + cos_o * sin_I * y_prime return x_ecl, y_ecl, z_ecl end #=========================================================================== # # Process the command-line arguments; convert the user-specified date into # the number of centuries past J2000. # raise ArgumentError, "Usage: #{$0} " unless (ARGV.size > 0) _T_eph = date_str_to_ephemeris_date(ARGV[0]) _T = (_T_eph - 2451545.0) / 36525.0 # # Determine the distances to each planet. # print " distance distance\n" print "planet (AU) (10^6 km)\n" print "------- -------- --------\n" Planets.each do |planet| next if (planet == 'Earth') # # Get the coordinates of the target planet and that of the Earth. # if (planet == 'Sun') tx = ty = tz = 0.0 else tx, ty, tz = position_point(keplerians[planet], _T) end ex, ey, ez = position_point(keplerians['Earth'], _T) # # And compute and print the distances! # x_diff = ex - tx ; y_diff = ey - ty ; z_diff = ez - tz dist_au = Math.sqrt(x_diff * x_diff + y_diff * y_diff + z_diff * z_diff) dist_Mkm = dist_au * 149.597_870_700 printf "%-7s %8.5f %8.3f\n", planet, dist_au, dist_Mkm end # # From the NASA Jet Propulsion Laboratory's HORIZONS application # (http://ssd.jpl.nasa.gov/horizons.cgi) for 01/16/2009: # # planet HORIZONS this code error # ---------- ---------- ----------- ----- # Sun 0.98373 au 0.98368 au 0.01% # Mercury 0.79262 au 0.80194 au 1.18% # Venus 1.71118 au 1.71107 au 0.01% # Mars 0.67940 au 0.67826 au 0.17% # Jupiter 5.78717 au 5.79397 au 0.12% # Saturn 9.08126 au 9.06789 au 0.15% # Uranus 20.60344 au 20.60612 au 0.01% # Neptune 30.87790 au 30.87685 au 0.00% # Pluto 32.67755 au 32.67345 au 0.01% # --Boundary_(ID_GkqsRnC8mcPrCjllCkjopw)--