Orbit-1 ******************************************************************************* Modified 5/26/2006 Prepared by Graham Kendall grahamkendall74135@yahoo.com ======= http://www.grahamkendall.net ======= http://www.aerospaceweb.org/question/spacecraft/q0164.shtml http://scienceworld.wolfram.com/physics/Two-BodyProblem.html http://mathworld.wolfram.com/Ellipse.html http://scienceworld.wolfram.com/physics/EllipticalOrbit.html http://www.coastalbend.edu/acdem/math/campbel7/stars/index.htm http://mathworld.wolfram.com/KeplersEquation.html http://www.geocities.com/SiliconValley/2902/kepler.htm http://www.akiti.ca/KeplerEquation.html http://www.willbell.com/math/mc12.htm http://www.projectpluto.com/kepler.htm http://www.xylem.f2s.com/kepler/kepler.html http://www.astro.uu.nl/~strous/AA/en/reken/kepler.html http://www.math.ubc.ca/~cass/courses/m309-8a/java/m309gfx/kepler.html http://en.wikipedia.org/wiki/Kepler's_laws_of_planetary_motion http://www.nezumi.demon.co.uk/astro/kepler/kepler.htm http://en.wikipedia.org/wiki/Space_mathematics Fundamentals of Astrodynamics, Bate, Mueller, White 2. Methods of Orbit Determination for the Microcomputer, Boulet 3. Astronomical Algorithms, Meeus 4. Methods of Orbit Determination, Escobal 5. Celestial Mechanics, Moulton 6. Modern Astrodynamics: fundamentals and perturbation methods, Bond, Allman 7. Adventures in Celestial Mechanics, Szebehely Astronomical Formulae For Calculators Meeus Celestial Basic by ??? ====== ----- http://scienceworld.wolfram.com/physics/topics/CelestialMechanics.html Much orbit math data http://lheawww.gsfc.nasa.gov/users/ddavis/courses/phys315/week3/week3.pdf advanced orbit http://cassfos02.ucsd.edu/physics/ph7/ch4.html good orbit lesson == the AU is the unit of distance, which is the distance earth to sun. q = perihelion= 0.4255 AU for this example. e = eccentricity = 0.2 for this example. a = semimajor axis = a = q/(1-e) = 0.531875 b = a*SQRT[1-e^2] = .521128942777 Area of the ellipse = Pi*a*b The area is 0.870772377707 c = a*e = 0.106375 = Distance center to focus. k = 0.01720209895 Gravity constant. r= 0.572141626031 c^2=a^2-b^2 1/e^2=1+(b/c)^2 New we prepare to solve for position when we know these quantities and elapsed time. There is a quantity called mean anomaly which goes from 0 to 2 * Pi in one orbit. M=k*t/a^1.5 assume t = time in days since passing perihelion point = 40 days for this example. Compute M M= 1.77389155705 The following is called Kepler's Equation. Solve for EE EE is always solved in radians in Kepler's Equation. M=EE-e*Sin[EE] where EE is eccentric anomaly in radians. e= eccentricity and m= mean anomaly. Calculators and computers have root finders. This method will work in any case. EE=M first approximation EE = EE + ( EE - e * Sin[EE] - M )/(e * Cos[EE] - 1 ) Repeat this until the value of EE reaches a final value.. Eccentric anomaly is an angle around the center of the ellipse. EE goes from 0 to 2Pi in one orbit. Draw an ellipse with long axis being horizontal. Draw a circle around it, touching at the left and right side, with the same center as the ellipse. Assume the planet is at a point in the ellipse. Mark this point on the ellipse with a 0. Draw a vertical line upward from this 0 until it reaches the circle at X. Draw a line from X to the center of the circle C. The true anomaly is the angle clockwise around F from the perihelion to the planet. * X | * o| * Ellipse | * Circle *** o | * 0 Planet * * * o # * * * * # * * o True Anomaly * * # ** o EE Focus ** C--------------F-------------| Perihelion Center |<-- q ------>| <------ c ---->| <--------------- a --------->| See D10-Kepler Jprg When we solve Kepler's Equation, we get EE= 1.95900897924 radians Now you can find the radius and true anomaly. Switch to degrees after this step and use th in degrees for the rest of the problem, if you choose. We use a, e, and EE to find radius and true anomaly angle th around the focus. Use radian mode s1 = r * Cos[th] = a * ( Cos[EE] - e ) = -.3077708130153 s2 = r * Sin[th] = a * (SQRT[ 1 - e^2 ]) * Sin[EE] = .482350232585 Note that EE is in radians and th is in degrees. Compute the right side in radians and switch to degrees when computing th. Remember Tan[th]= s1/s2 If you use ArcTan function, add 180 degrees to th if s2 is negative. If your calculator has rectangular to polar function, it computes the proper quadrant. r = SQRT[s1^2+s2^2] r = .5721416260 AU and th = 122.535231561 degrees. Go to degree mode th= ArcTan[s2/s1] + 180 since s1 is negative The goal is to find Right Ascension (RA) and Declination (Dec) of the object as seen from earth. RA is the angle in degrees around the equator west to east starting from the equinox point, where the sun crosses the equator in the spring. Ra is in hours, minutes and seconds in the star position tables. 360 degrees = 24 hours. Declination is degrees from the equator, + north and - south. After we get r and th, we go through a series of coordinate rotations. Computers normally compute totally in radians. Some can directly use degrees. Calculators can handle both. r= q * (1+e)/(1 + e * Cos[th] ) = p / (1 + e * Cos[th] ) I will list the six orbit elements. q=perihelion. e= eccentricity DATE = Date of perihelion so we can get t, days after perihelion. The other three are angles. I = inclination of the orbit to the plane of the earth's orbit in degrees. It is 0 to 90 degrees if the planet goes the same direction as the other planets, otherwise it is 90 to 180 degrees. Another one called W, the argument of the perigee. The orbit as seen from the north comes up through the plane of the earth's orbit at a point which I call PP, the puncture point. Now consider the plane of the planet orbit. There is the perihelion point and the puncture point PP. The angle between them is W. It is positive if the planet reaches PP before it meets q. There is a reference point in space called the Equinox point. It is the place on the equator where the sun crosses it in the spring. It changes slowly over time due to the precession of the earth, one cycle each 26,000 years. There is an angle in the plane of the earth orbit between this Equinox point and PP. This is the Longitude of the Ascending Node, called capital omega in the literature and O in my equations, in degrees. These are the six orbit elements. We will now rotate the orbit through I,W, O and IE. In my example q = 0.4255 AU e = 0.2 t = 40 days I = 72 degrees = inclination W = 105 degrees = argument of perigee O = 293 degrees = longitude of the ascending node IE = tilt of earth = 23.441028 degrees Rotate the ellipse around W. U = th + W = 227.535231561 degrees x1 = r * Cos[U] = -0.386273822557 AU y = r * Sin[U] = -0.422064656473 AU Now rotate ellipse around I. y1 = y * Cos[I] = -0.130425151575 AU z2 = y * Sin[I] = -0.401407341836 AU Now rotate ellipse around O. x3 = x1 * Cos[O] - y1 * Sin[O] = -0.27098619163 AU y2 = y1 * Cos[O] + x1 * Sin[O] = 0.304605761767 AU Now rotate ellipse around IE. y3 = y2 * Cos[IE] - z2 * Sin[IE] = 0.439148484296 AU z3 = z2 * Cos[IE] + y2 * Sin[IE] = -0.247105509699 AU The position of the sun in equatorial coordinate system is XS = -0.931108260968 AU YS = 0.371439715781 AU ZS = 0.161052202235 AU This is an earth centered system. This changes continually. We change origin from the sun to the earth. This computes earth to planet dustance. s1 = x3 + XS = -1.2020944526 AU s2 = y3 + YS = 0.810587642107 AU s3 = z3 + ZS = -0.086053307464 AU s1 = ro * Cos[RA] * Cos[Dec] = -1.2020944526 s2 = ro * Sin[RA] * Cos[Dec] = 0.810588200078 s3 = ro * Sin[Dec] = -0.086053307464 ro = earth to planet distance in AU. ro= SQRT[s1^2 + s2^2 + s3^2] ******************************************************************************* Tan[RA] =s2/s1 If s1 negative add 180 degrees to computed value of RA, if ArcTan[s2/s1] is computed with a calculator or computer. If rectangular to polar computation is used, then RA is in the proper quadrant. Sin[Dec] = s3/ro ro = 1.45240816398 AU RA = 146.007690781 degrees Dec = -3.3966901959 degrees This is now completed. ------ Find area of ellipse a=q/(1-e) b=a*sqrt(1-e^2) area=pi*a*b = 0.870772377706 in polar coordinates area = 0.5 * Integral { 0 to 2pi} r^2 dtheta area = 0.5 * Integral {0 to 2pi} q^2*(1+e)^2/(1+e*Cos(theta))^2 d theta =0.870772377706