Trajectory and Orbits - Computational Physics - Lecture Notes | PHYS 4007, Study notes of Physics

Download Trajectory and Orbits - Computational Physics - Lecture Notes | PHYS 4007 and more Study notes Physics in PDF only on Docsity! PHYS-4007/5007: Computational Physics Course Lecture Notes Section X Dr. Donald G. Luttermoser East Tennessee State University Version 4.1 Abstract These class notes are designed for use of the instructor and students of the course PHYS-4007/5007: Computational Physics taught by Dr. Donald Luttermoser at East Tennessee State University. Donald G. Luttermoser, ETSU X–3 Figure X–2: The spherical-polar coordinate system. y z x P r r̂ θ̂ φ̂ θ φ Figure X–3: The circular-cylindrical coordinate system. y z x P ρ ρ z z ρ̂ ẑ φ̂ φ X–4 PHYS-4007/5007: Computational Physics c) Circular-Cylindrical (Cylindrical) Coordinates (ρ, φ, z): ~∇ = ρ̂ ∂ ∂ρ + φ̂ 1 ρ ∂ ∂φ + ẑ ∂ ∂z . (X-5) 4. In this section, we will be using the gravitational potential field to describe the potential energy, then from Newton’s Universal Law of Gravity, ~Fg(r) = − Gm1m2 r2 r̂ , (X-6) we see that this force only depends upon r. Here the nega- tive sign indicates that the gravitational force is attractive, with the m’s representing the masses of bodies 1 and 2 which are separated by a distance r. G = 6.673 × 10−11 N m2/kg2 = 6.673× 10−8 dyne cm2/g2 is the Universal Gravitation Constant. a) Using Eq. (X-6) in conjunction with Eq. (X-4), we imme- diately see that U = U(r) since F = F (r). i) To solve the gradient equation for the potential en- ergy (Eq. X-2), we can take the dot product (i.e., inner product) of both sides of the equation of Eq. (X-2) with the differential displacement d~r and in- tegrate from the initial (1) to the final position (2). ii) Since the force here is conservative, the work done by this force is independent of the path or trajec- tory taken. Hence, we can use the standard definite integral instead of a path integral to solve for U : ∫ 2 1 ~F · d~r = − ∫ 2 1 ~∇U · d~r (X-7) − ∫ 2 1 ( Gm1m2 r2 r̂ ) · (dr r̂) = − ∫ 2 1 ( ∂U ∂r r̂ ) · (dr r̂) ∫ 2 1 Gm1m2 r2 dr = ∫ 2 1 dU dr dr (X-8) where the partials become full derivatives since U = U(r). Donald G. Luttermoser, ETSU X–5 iii) Integrating Eq. (X-8), we get − Gm1m2 r 2 1 = ∫ 2 1 dU = U2 − U1 , or U1 − U2 = Gm1m2 ( 1 r2 − 1 r1 ) . (X-9) b) From Eq. (X-9), we can immediately see that the potential energy for a gravitational field takes on the form U = −Gm1 m2 r . (X-10) i) From this equation we see that U → 0 as r → ∞. ii) Also we can see that gravitational potential en- ergy is a negative energy. 5. Besides potential energy, any body in motion has a kinetic energy associated with it. a) The left-hand side of Eq. (X-7) is the definition of the work done on a particle by a force field: W ≡ ∫ 2 1 ~F · d~r = U1 − U2 . (X-11) b) Using one of the forms of Newton’s Second Law of Motion in Eq. (X-1) and rewriting d~r as (d~r/dt)dt, we can write ~F · d~r = ( m d~v dt ) · ( d~r dt dt ) = m d~v dt · ~v dt . Note that d dt (~v · ~v) = ~v · d~v dt + d~v dt · ~v = 2d~v dt · ~v , so d~v dt · ~v = 1 2 d dt (~v · ~v) X–8 PHYS-4007/5007: Computational Physics 8. When dealing with motion in gravitational fields, there are two regimes that are typically encountered: (1) trajectories (near a surface of a large mass, e.g., Earth) and (2) orbits (where 2 masses can be considered as point-like). a) For orbits, we use the general form of the gravitational potential as described in Eqs. (X-10,16): Ug = − GMm r , (X-21) where we are now using M (the larger mass) for m1 and m (the smaller mass) for m2 and Φ = − GM r . (X-22) b) For trajectories, typically the maximum height (ymax = h) reached is small with respect to R⊕ and hence g ≈ constant. As such, we can write Eqs. (X-9 and X-11) as W = U1 − U2 = GM⊕m ( 1 r2 − 1 r1 ) . i) If we take point ‘2’ to be the Earth’s surface and ‘1’ to be the position of the projectile, then ∆U = GM⊕m ( 1 R⊕ − 1 R⊕ + h ) = GM⊕m   R⊕ + h R⊕ (R⊕ + h) − R⊕ R⊕ (R⊕ + h)   = GM⊕m  R⊕ + h− R⊕ R⊕ (R⊕ + h)   = GM⊕m   h R⊕ (R⊕ + h)   . ii) If R⊕ is much greater than h (which it will be for experiments near the Earth’s surface), h  R⊕. As Donald G. Luttermoser, ETSU X–9 such, R⊕+h ≈ R⊕ and the equation above becomes ∆U = GM⊕m   h R⊕ (R⊕)   = GM⊕m   h R2⊕   = GM⊕mh R2⊕ = m GM⊕ R2⊕ h . (X-23) iii) Using Eqs. (X-15 and X-16), we see that ~g = −~∇Φ = −~∇ ( GM⊕ r ) = − d dr ( GM⊕ r ) r̂ = GM⊕ r2 r̂ (X-24) and at the Earth’s surface, g = GM⊕ R2⊕ , (X-25) where ‘g’ is referred to as the Earth’s surface grav- ity. iv) Using Eq. (X-25) in Eq. (X-23), we finally get ∆U = mgh = mg∆y , where ∆y = y− y◦ = h is just the change in height from our initial position y◦ (typically the ground) and y is an arbitrary position in the trajectory above y◦. v) If we arbitrarily set y◦ = 0, then the potential at that position is zero, and y represents the posi- tion above the ground (y◦). As such, the potential energy becomes U = mgy . (X-26) X–10 PHYS-4007/5007: Computational Physics 9. Trajectory calculations can be difficult due to non-gravitational forces that enter the calculations. a) A drag force due to air friction which usually takes the form ~Fr = −mkvn ~v v , (X-27) where ~Fr represents the retarding (i.e., drag) force, v is the magnitude of the velocity, ~v is the velocity vector (hence ~Fr is in the opposite direction of ~v due to the negative sign, note that the ratio ~v/v is essentially just a unit vector in the direction of ~v), and k is the drag coefficient. As such, the total force on the object now becomes ~F = ~Fg + ~Fr . (X-28) b) If the downrange distance of the projectile is large enough such that the Earth’s surface can no longer be represented as a flat plane, the Earth’s rotation has to be taken into account. To an observer in the rotating coordinate sys- tem, the effective force (ignoring air friction) is ~Feff = m~af −m~ω × (~ω × ~r) − 2m~ω × ~vr . (X-29) i) ~Ff = m~af is the force in the fixed coordinate sys- tem (which is just Newton’s 2nd law). This force is said to be an inertial force since it applies only to a static coordinate system. ii) ~Fcf = −m~ω × (~ω × ~r) is the centrifugal force, which results from trying to write an inertial force law for a noninertial (i.e., accelerating) reference frame (note that ~ω is called the angular velocity). The minus sign in this term implies that this pseudo- force (see below) is directed outwards from the cen- Donald G. Luttermoser, ETSU X–13 i) Since v2 = v2x + v 2 y , we can write Fr,x = Fr cos θ = Fr vx v , (X-36) where θ is the angle between ~vx and ~v. ii) Likewise Fr,y = Fr sin θ = Fr vy v . (X-37) iii) Hence, the drag force has components of Fr,x = −mkvvx Fr,y = −mkvvy . (X-38) d) The x-component of Newton’s 2nd law gives dvx dt = −kvvx (X-39) with dx dt = vx . (X-40) e) The y-component gives dvy dt = −g − kvvy (X-41) with dy dt = vy . (X-42) f) The solution to these first-order DEs can be found with simple forward-difference equations: xi+1 = xi + vx,i ∆t vx,i+1 = vx,i − kvvx,i ∆t yi+1 = yi + vy,i ∆t vy,i+1 = vy,i − g∆t− kvvy,i ∆t . (X-43) X–14 PHYS-4007/5007: Computational Physics g) Eqs. (X-43) can be then solved as an initial value problem as described in §IX with x◦, y◦, vx,◦, and vy,◦ supplied by the user. i) The ∆t steps are chosen to give the errors that follow Eq. (VII-15). ii) Calculations are carried out until a certain τ = total time is reached or some condition of xi, yi, vx,i, or vy,i is satisfied. h) But what is the drag coefficient k? i) Since the projectile is trying to push air of mass dmair out of the way, where dmair ≈ ρA v dt , (X-44) where ρ is the density of the air and A is the frontal area, we can guess that k is a function of ρ [and possibly t if the object is rotating → then A = A(t)]. ii) k(y = 0) = k◦ is usually given (based on air tun- nel measurements), and the following expression is used for k(y): k(y) = ρ(y) ρ◦ k◦ , (X-45) where ρ◦ is the density of air when the k◦ measure- ment was made. i) Now we need a description of ρ(y)! i) One could supply a data table of ρ as a func- tion of height (see Appendix B.1 in Fundamentals of Atmospheric Modeling by Mark Jacobson, 1999, Cambridge University Press). Donald G. Luttermoser, ETSU X–15 ii) One could solve the following set of differential equations to determine ρ(y): P = NkBT (ideal gas law) (X-46) ρ = ∑ Nimi (mass conservation) (X-47) dP dy = −ρ g (hydrostatic equilibrium) (X-48) dT T =  R cp   dP P (Poisson’s equation) (X-49) dQ = cp dT − αdP (energy conserv.) (X-50) where P = pressure, T = temperature, N = total particle density, ρ total mass density, Ni = num- ber density of species i, mi = mass of species i, kB = Boltzmann’s constant, g = surface gravity, R = universal gas constant, cp = specific heat of air at constant pressure, Q = total heat (determined from solar radiation incident on the Earth’s atmo- sphere), and α = specific volume of air. iii) Note that the solution to these equation would still only be an approximation since we have left out condensation, evaporation, sublimation, chemical reactions, and wind from these equations. 2. Trajectories with h ∼ R⊕ and x ∼ R⊕. a) Rotating Coordinate Systems. i) Let’s consider 2 sets of coordinate axes: =⇒ one ‘fixed’ = inertial frame (the primed [′] coordinates), =⇒ one ‘rotating’ with respect to the fixed system and possibly in linear motion with respect to the fixed frame = noninertial frame. X–18 PHYS-4007/5007: Computational Physics  d~r ′ dt   fixed =   d~R dt   fixed + ( d~r dt ) rot + ~ω × ~r . If we define ~vf ≡ ~̇rf ≡  d~r ′ dt   fixed (X-58) ~V ≡ ~̇Rf ≡   d~R dt   fixed (X-59) ~vr ≡ ~̇rr ≡ ( d~r dt ) rot (X-60) we may write ~vf = ~V + ~vr + ~ω × ~r , (X-61) where ~vf = velocity relative to the fixed axes ~V = linear velocity of the moving origin ~vr = velocity relative to the rotating axes ~ω = angular velocity of the rotating axes ~ω × ~r = velocity due to the rotation of the moving axes. b) The Coriolis Force. i) Newton’s 2nd law ~F = m~a is only valid in an inertial frame, therefore ~F = m~af = m ( d~vf dt ) fixed , (X-62) where the differentiation must be carried out with respect to the fixed system. ii) If we limit ourselves to cases of constant angular acceleration (ω̇ = 0), using Eq. (X-61) we can write ~F = m~̈Rf +m ( d~vr dt ) fixed +m~ω× ( d~r dt ) fixed . (X-63) Donald G. Luttermoser, ETSU X–19 iii) The second term can be evaluated by substitut- ing ~vr for ~Q in Eq. (X-56): ( d~vr dt ) fixed = ( d~vr dt ) rot + ~ω × ~vr = ~ar + ~ω × ~vr , (X-64) where ~ar is the acceleration in the rotating coordi- nate system. iv) The last term in Eq. (X-63) can be obtained di- rectly from Eq. (X-55): ~ω × ( d~r dt ) fixed = ~ω × ( d~r dt ) rot + ~ω × (~ω × ~r) = ~ω × ~vr + ~ω × (~ω × ~r) . (X-65) v) Combining Eqs. (X-63)-(X-65), we obtain ~F = m~̈Rf+m~ar+2m~ω×~vr+m~ω×(~ω×~r) , ~̇ω = 0 (X-66) where ~̈Rf is the acceleration of the origin of the moving coordinate system relative to the fixed sys- tem. vi) In the case of trajectories on the Earth’s surface, the origin of the rotating coordinate system is sta- tionary (in the ~r direction) with respect to the fixed coordinate system. As such, ~̈Rf = 0 we can finally write ~F = m~af = m~ar+m~ω×(~ω×~r)+2m~ω×~vr . (X-67) vii) From this equation, the effective force on a par- ticle measured by an observer in the rotating frame X–20 PHYS-4007/5007: Computational Physics is then ~Feff = m~ar = m~af −m~ω × (~ω × ~r) − 2m~ω × ~vr . (X-68) which is what we wrote down in the Introduction of this section of the notes in Eq. (X-29). c) Coding Problems with h ∼ R⊕ and x ∼ R⊕ — Part 1: The Fixed Coordinate System. i) Unlike the preceding case where h  R⊕ and x  R⊕, for our current problem we must use a three- dimensional coordinate system. ii) In reality, the Earth is constantly being subjected to a variety of motions: =⇒ The Earth’s own rotation. =⇒ The Earth’s orbital velocity around the Sun. =⇒ The Sun’s orbital velocity about the center of the Milky Way Galaxy. =⇒ The Milky Way’s and Andromeda (M31) galaxy’s motion towards each other. =⇒ The Local Group of galaxy’s motion in the Virgo Supercluster of galaxies. =⇒ The Hubble flow of all of the galaxies, clusters, and superclusters in the Universe. As can be seen, quite a lot of velocities to worry about. However, for our problem here, the Earth’s rotational velocity will dominate these other veloc- ities. iii) Choose the fixed coordinate system’s origin to lie at the Earth’s center and its ẑ′ direction to lie from a line from the center to the North Pole (follow- Donald G. Luttermoser, ETSU X–23 x) The Earth is subjected to extremely small, but unpredictable, variations in its rotation rate (mainly due to gravitational perturbations from other bod- ies in the Solar System and a non-smooth slow- ing down due to the tides raised from the Sun and Moon). To precisely predict the positions of bod- ies in the Solar System, we require a steady time standard. =⇒ Universal Time is then replaced by Ephemeris Time (E.T.) in celestial mechanics. =⇒ At the beginning of 1900 A.D., an ephemeris second was defined as 1/31,556,925.97474 the length of the tropical year 1900 and both U.T. and E.T. were in agreement. =⇒ Today, these times differ by about 56 seconds. xi) As such, to calculate φ′, use the following for- mula: φ′ = 2π tU.T. 24 hr , (X-70) where tU.T. is the current Universal Time in decimal sidereal hours. xii) As such, your code should have the ability for the user to input: =⇒ A time of launch (in U.T. sidereal hours). =⇒ The latitude of launch (in decimal degrees). =⇒ The altitude of the ground with respect to sea level. X–24 PHYS-4007/5007: Computational Physics d) Coding Problems with h ∼ R⊕ and x ∼ R⊕ — Part 2: The Rotating Coordinate System. i) We now have all of the fixed (primed ) frame coor- dinates (i.e., r′, λ′, and φ′), we now need the rotat- ing frame coordinates. =⇒ x is defined in the eastern direction. =⇒ y is defined in the northern direction. =⇒ z is the altitude (which follows from the right- hand rule). ii) Rotational coordinate frame transformations are then made with r = ( x2 + y2 + z2 )1/2 (X-71) θ = cos−1   z (x2 + y2 + z2)1/2   (X-72) φ = tan−1 y x , (X-73) where r is measured from the launch point (the rotating origin), θ is the angle from the z (i.e., al- titude) axis (note that this makes θ different from what it was in the h  R⊕ case when it was mea- sured from the ground), and φ is the angle sub- tended from the east point on the horizon moving towards north. Donald G. Luttermoser, ETSU X–25 iii) Now your code should allow the user to input: =⇒ A launch velocity (typically in m/s) = ṙ. =⇒ The projection angle of launch (in decimal degrees) = γ. =⇒ The direction angle with respect to the east direction rotating towards north = φ. iv) Then from these in input, r will be determined in the next step and θ = 90◦ − γ. e) Coding Problems with h ∼ R⊕ and x ∼ R⊕ — Part 3: Solving the System of Equations. i) Calculate the Coriolis and centrifugal cross prod- ucts in Eq. (X-68). ii) Solve for vr in Eq. (X-61) by using the forward- difference technique shown in Eqs. (X-43). How- ever, now, replace y with z and the two x equations become four — 2 for x and 2 for y which will now have the same functional form. Also, instead of g, now use ar from Eq. (X-68) in the z direction where af = g in this equation. iii) Solve for the unknowns: downrange distance, max- imum height, impact velocity, etc. X–28 PHYS-4007/5007: Computational Physics Figure X–5: Geometry setup of Kepler’s 2nd Law of Motion. e) Referring to Fig. (X-5), we see that describing the path ~r(t), the radius vector sweeps out an area 12r 2 dθ in a time interval dt: dA = 1 2 r2 dθ , (X-84) and dividing by the time interval, the areal velocity is dA dt = 1 2 r2 dθ dt = 1 2 r2θ̇ = ` 2µ = constant, (X-85) which shows that the areal velocity is constant in time =⇒ Kepler’s Second Law of Planetary Motion. f) The conservation of energy gives us E = T + U = constant (X-86) = 1 2 µ ( ṙ2 + r2θ̇2 ) + U(r) , (X-87) or E = 1 2 µṙ2 + 1 2 `2 µr2 + U(r) . (X-88) Donald G. Luttermoser, ETSU X–29 3. Equation of Motion. a) Using the total energy equation (Eq. X-88), we can solve for the velocity: ṙ = dr dt = √√√√2 µ (E − U) − ` 2 µ2r2 . (X-89) b) Note that dθ = dθ dt dt dr dr = θ̇ ṙ dr = `/µr2 ṙ dr , (X-90) or ṙ = ` µr2 dr dθ , (X-91) hence ` µr2 dr dθ = √√√√2 µ (E − U) − `2 µ2r2 . (X-92) c) This leads to an equation of θ(r): θ(r) = ∫ (`/r2) dr √ 2µ(E − U − `2/2µr2) + constant. (X-93) d) If F (r) ∝ rn, Eq. (X-93) becomes an elliptical integral. i) For n = 1,−2, and −3 → solution are circular functions. ii) n = 1 gives harmonic motion. iii) n = −2 gives a central-force law (e.g., gravity, see later). iv) n = −3 gives an equation that is not important in physics. X–30 PHYS-4007/5007: Computational Physics e) We can also solve this problem with Lagrangians: ∂L ∂r − d dt ∂L ∂ṙ = 0 , (X-94) where the Lagrangian is given by Eq. (X-80). i) Term one of Eq. (X-94) is ∂L ∂r = µrθ̇2 − ∂U ∂r . ii) Term two of Eq. (X-94) is ∂L ∂ṙ = µṙ , d dt ∂L ∂ṙ = µr̈ . iii) Putting these two equations together gives µrθ̇2 − µr̈ − ∂U ∂r = 0 or µ(r̈ − rθ̇2) = − ∂U ∂r = F (r) . (X-95) f) Now let u ≡ 1 r , (X-96) then du dθ = − 1 r2 dr dθ = − 1 r2 dr dt dt dθ = − 1 r2 ṙ θ̇ . (X-97) g) Since θ̇ = `/µr2, we get du dθ = −µ ` ṙ . (X-98) h) Next, d dθ du dθ = d2u dθ2 = d dθ ( −µ ` ṙ ) = dt dθ d dt ( −µ ` ṙ ) = − µ `θ̇ r̈ = −µ 2 `2 r2 r̈ . (X-99) Donald G. Luttermoser, ETSU X–33 d) For gravity, F (r) = − Gm1m2 r2 , (X-110) so U(r) = − ∫ F (r) dr = − Gm1m2 r (X-111) and V (r) = − Gm1m2 r + `2 2µr2 , (X-112) where V → 0 as r → ∞. See Fig. (X-6) on the next page for plots of these potentials. i) For the bottom plot in Fig. (X-6), if E = E1 > 0, the motion is unbounded. ii) If 0 > E(= E2) > E3, the motion has turning points (i.e., bounded between r2 and r4 → the ap- sidal distances). iii) If E = E3, only one solution exists → motion is circular. iv) Note that E < E3 is not possible =⇒ ṙ would be imaginary! 6. Planetary Motion — Kepler’s Problem. a) Let us re-examine Eq. (X-93). i) For the potential energy given by Eq. (X-111) let’s define k = Gm1m2 (typically we will let m1 be the bigger mass such that m1  m2) giving the poten- tial energy as U(r) = − Gm1m2 r = − k r . (X-113) X–34 PHYS-4007/5007: Computational Physics Figure X–6: The various potentials in a gravitational field. Donald G. Luttermoser, ETSU X–35 ii) If we choose the origin of θ so that the integration constant in Eq. (X-93) is zero, we can show (with a little algebra) that cos θ = `2 µk · 1 r − 1 √ 1 + 2E` 2 µk2 . (X-114) iii) Let us now define the following constants: α ≡ ` 2 µk (X-115) ε ≡ √√√√1 + 2E`2 µk2 . (X-116) b) Then we can rewrite Eq. (X-114) as α r = 1 + ε cos θ . (X-117) i) This is the equation of a conic section with one focus at the origin. ii) The quantity ε is called the eccentricity of the orbit. iii) The quantity 2α is termed the latus rectum of the orbit. c) The minimum value for r occurs when cos θ is a maxi- mum (i.e., for θ = 0). Thus the choice of the integration constant in Eq. (X-93) be zero corresponds to measuring θ from rmin. i) For an arbitrary orbit, rmin is called the pericen- ter. ii) For solar orbits it is called the perihelion. X–38 PHYS-4007/5007: Computational Physics ii) The semiminor axis of an elliptical orbit is given by b = α√ 1 − ε2 = ` 2µ|E| , (X-120) as such, the semiminor axis depends both on E and `. Comparing Eq. (X-118) with Eq. (X-120), we can write the semiminor axis in terms of the semimajor as b = a √ 1 − ε2 = √ αa . (X-121) g) The extrema (i.e., turning points) of the orbit is given by rmin = a(1 − ε) = α 1 + ε (X-122) rmax = a(1 + ε) = α 1 − ε (X-123) h) To find the period of elliptical motion, take the areal ve- locity equation (Eq. X-85) and integrate: dt = 2µ ` dA ∫ τ 0 dt = 2µ ` ∫ A 0 dA τ = 2µ ` A , (X-124) where the entire area of the ellipse A is swept out in one orbital period τ . i) Since A = πab for an ellipse, we can write the orbital period in terms of the energy of a particle in orbit a τ = 2µ ` πab = 2µ ` k 2|E| ` √ 2µ|E| = πk √ µ 2 |E|−3/2 . (X-125) Donald G. Luttermoser, ETSU X–39 j) Likewise, we can use Eqs. (X-115) and (X-121) to write τ = 2µ ` πab = 2µ ` πa √ αa = 2πµ ` √√√√ ` 2 µk a3/2 , (X-126) or τ 2 = 4π2µ k a3 = 4π2 G(m1 +m2) a3 , (X-127) which is the analytic proof to Kepler’s Third Law of Planetary Motion. Note that for planets in the solar system, m1 = M  m2 (where m2 is the mass of a planet), giving τ 2 ∼= 4π2 GM a3 = K a 3 , where K is a constant (as long as m2 M ) =⇒ this is the 3rd law form as Kepler stated it. 7. Kepler’s Equation. a) We have written an equation of r(θ) (see Eq. X-117). We will now develop an equation for θ, called the true anomaly, as a function of time. b) We will start by setting up the ratio A τ = πab τ = dA dt . (X-128) i) Let θ = 0 at t = 0, where dA = 12r 2dθ, then πab τ ∫ t 0 dt′ = ∫ A 0 dA′ = 1 2 ∫ θ 0 r2 dθ . (X-129) ii) Using Eq. (X-117), we can write r in terms of θ as r = α 1 + ε cos θ , (X-130) X–40 PHYS-4007/5007: Computational Physics and hence πab τ t = α2 2 ∫ θ 0 dθ (1 + ε cos θ)2 = α2 2(1 − ε2)   2√ 1 − ε2 tan−1  (1 − ε) tan(θ/2)√ 1 − ε2   − ε sin θ 1 + ε cos θ ] . iii) Note that ab = α2 ( 1 − ε2 )−3/2 , (X-131) so we can write 2πt τ = 2 tan−1   √√√√1 − ε 1 + ε tan θ 2   −ε √ 1 − ε2 sin θ 1 + ε cos θ . (X-132) c) As can be seen, retrieving θ(t) from Eq. (X-132) will not be easy! Due to this difficulty, we will approach the prob- lem through analytic geometry. i) We will start by circumscribing the ellipse (repre- senting an orbit) with a circle and define a coordi- nate system with the origin set at one of the foci of the ellipse as shown in Figure X-9. ii) The equation for the ellipse is (x+ aε)2 a2 + y2 b2 = 1 . (X-133) • The circle has radius of ‘a.’ • Point ‘P ’ lies on the ellipse, ‘Q’ on the circle. Donald G. Luttermoser, ETSU X–43 j) If we make use of the trigonometric identity tan α 2 = √√√√1 − cosα 1 + cosα , we can take the square root of Eq. (X-143) to write tan θ 2 = √√√√1 + ε 1 − ε tan ψ 2 . (X-144) k) Now, θ(t) can be found directly from ψ(t). i) Differentiating Eq. (X-144) gives dθ = √√√√1 + ε 1 − ε cos2(θ/2) cos2(ψ/2) dψ . (X-145) ii) Rewriting Eq. (X-141) gives r = a(1 − ε)1 + cosψ 1 + cos θ = a(1 − ε) cos2(ψ/2) cos2(θ/2) . (X-146) l) We will now make use of this expression for r in Eq. (X- 129): πab τ t = 1 2 ∫ θ 0 r2 dθ . i) However for this integral, let r2 = (r1)(r2), where r1 is given by Eq. (X-138) and r2 is given by Eq. (X-146). ii) Then, using Eq. (X-145) for dθ we get r2 dθ = [a(1 − ε cosψ)]  a(1 − ε)cos 2(ψ/2) cos2(θ/2)   ·   √√√√1 + ε 1 − ε · cos 2(θ/2) cos2(ψ/2) dψ   = a2 √ 1 − ε2 (1 − ε cosψ) dψ . X–44 PHYS-4007/5007: Computational Physics iii) So πab τ t = a2 √ 1 − ε2 2 ∫ ψ 0 (1 − ε cosψ) dψ . iv) Since b = a √ 1 − ε2 , ab = a2 √ 1 − ε2 , and 2πt τ = ψ − ε sinψ . (X-147) v) Let M = 2πt/τ = mean anomaly =⇒ measures angular deviation of a body moving in a circular orbit of period τ (note that M = 2π when the mass completes one orbit, t = τ ), then M = ψ − ε sinψ , (X-148) which is referred to as Kepler’s Equation. m) In order to find ψ(t), Kepler’s Equation must be inverted by some approximation procedure which usually involves expanding the sine term in a Taylor’s series (see E.W. Brown 1931, Monthly Notices of the Royal Astronomical Society, 92, 104). Then, Eq. (X-144) relates ψ and θ and hence the time dependence of the true anomaly (i.e., the anomaly for an elliptical orbit) can be found. 8. Following the same geometric arguments that we made to de- rive Kepler’s Equation, we can develop a velocity equation as a function of r. The orbital velocity can be found with v2 = ẋ2 + ẏ2 , (X-149) or v2 = a2ψ̇2 sin2ψ + a2(1 − ε2)ψ̇2 cos2ψ = a2ψ̇2(1 − ε2 cos2 ψ) . (X-150) Donald G. Luttermoser, ETSU X–45 a) Differentiating Eq. (X-147) gives 2π τ = ψ̇(1 − ε cosψ) . (X-151) b) Solving this equation for ψ̇ and inserting this value in Eq. (X-150) gives v2 = ( 2π τ )2 a2 1 − ε2 cos2ψ (1 − ε cosψ)2 = ( 2π τ )2 a2 1 + ε cosψ 1 − ε cosψ = ( 2π τ )2 a2 2 − (1 − ε cosψ) 1 − ε cosψ . (X-152) c) Note that r a = 1 − ε cosψ , so we can write v2 = ( 2π τ )2 a3 ( 2 r − 1 a ) . (X-153) d) Using Kepler’s 3rd Law: v2 = Gm1m2 µ ( 2 r − 1 a ) or v2 = G (m1 +m2) ( 2 r − 1 a ) . (X-154) D. Programming Orbital Motion in the Solar System. 1. As we have seen, the analytic solution of planetary motion is extremely complicated. This problem is much easier to solve numerically. 2. Here, we will limit ourselves to motion of a mass (m2 = mp = m) in orbit about the Sun (m1 = M ). Using these masses in Eq. X–48 PHYS-4007/5007: Computational Physics xi+1 = xi + vx,i+1 ∆t (X-165) vy,i+1 = vy,i − GM yi r3i ∆t (X-166) yi+1 = yi + vy,i+1 ∆t . (X-167) 5. To carry out an orbit calculation, one must first define the semi- major axis size ‘a’ and the eccentricity ε of the orbit as input parameters. In addition, the masses of the primary (e.g., the Sun for solar orbits) and of the secondary must be supplied. 6. Finally, we need to set the initial values for x, y, vx, and vy, and the initial time interval ∆t. This is usually done by starting at the pericenter point with θ = θ◦ = 0. a) Then Eq. (X-122) gives us the value of r: r◦ = rmin = a(1 − ε) . b) From this, the starting Cartesian displacements are x◦ = r◦ cos θ◦ = r◦ = rmin (X-168) y◦ = r◦ sin θ◦ = 0 . (X-169) c) The velocity of an object in an elliptical orbit is given by Eq. (X-154). For r = rmin, we have v◦ = √√√√G (M +m) ( 2 rmin − 1 a ) . (X-170) As such, the initial velocities in Cartesian coordinates are given by vx,◦ = 0 (X-171) vy,◦ = v◦ (X-172) since at θ = 0 (i.e., r = rmin), the orbital velocity vec- tor points entirely in the +y direction (see Figure X-10). Donald G. Luttermoser, ETSU X–49 Note that this will cause the mass to move in the counter- clockwise direction as drawn in Figure X-10. 7. One then numerically solves Eqs. (X-164) through (X-167) us- ing the Euler-Cromer or 4th order Runge-Kutta (RK4) methods. The best results are usually obtained with RK4 with an adaptive grid (i.e., time interval). See the Numerical Recipes books for details.