HELPFUL DESCRIPTIONS. a : semimajor axis For elliptical orbits, the semimajor axis is half the length of the longest line segment that can be drawn inside the ellipse. For hyperbolic orbits, the semimajor axis is the distance from the vertex to the nearest point on either branch. e : eccentricity The orbit is a circle when e=0. The orbit is an ellipse when 0<e<1. The orbit is a parabola when e=1. The orbit is a hyperbola when e>1. i : inclination to ecliptic I prefer to use the principal value of the arc-cosine for my range of inclination [0 , pi radians]. I take the ecliptic to be the current plane of Earth's orbit. (This is the hillbilly tutorial, remember?) L : longitude of ascending node The angle, subtended at the sun, measured in the ecliptic, from the Vernal Equinox to the ascending node. The ascending node is the point at which the orbit crosses the ecliptic having the Z component of its velocity in the direction of the North Ecliptic Pole. (Most people use a capital omega for this quantity.) w : argument of perihelion The angle, subtended at the sun, measured in the plane of the orbit in the direction of motion, from the ascending node to the orbit's perihelion. (Most people use a lower-case omega for this quantity.) T : time of perihelion passage The time (e.g. a decimal calendar or Julian date) at which the orbiting object passes through the perihelion of its orbit. Q : true anomaly The angle, subtended at the sun, measured in the plane of the orbit in the direction of motion, from the perihelion to the current position of the orbiting object. No Greek letters can be rendered here, sorry. M : mean anomaly The angle, subtended at the sun, measured in the plane of the orbit in the direction of motion, from the perihelion to the point at which the orbiting object would be found, if it moved constantly by its average angular velocity. u : eccentric anomaly The angle, subtended at the geometric center of the orbit, measured in the plane of the orbit in the direction of motion, from the perihelion to a point on the circle circumscribing the orbit, which point is found by extending a line perpendicular to the major axis through the current position of the orbiting object until it intersects the circumscribing circle. No pictures, sorry. WHERE TO FIND ORBITAL ELEMENTS. Planetary orbital elements (NASA/JPL): http://ssd.jpl.nasa.gov/elem_planets.html Except don't believe the values for Earth's inclination and longitude of ascending node. Set them both equal to zero. The other elements in that table are OK. Why they didn't just go ahead and give the times of perihelion passage (or a mean anomaly at epoch), I've no idea. Planetary orbital elements (some other guy): http://www.astro.uu.nl/~strous/AA/en/reken/hemelpositie.html Don't believe his value for Earth's longitude of ascending node, either. Set it to zero. Asteroid orbital elements: http://ssd.jpl.nasa.gov/sb_elem.html Convert a calendar date to a Julian date: http://wwwmacho.mcmaster.ca/JAVA/JD.html Convert a Julian date to a calendar date: http://wwwmacho.mcmaster.ca/JAVA/CD.html Jerry Abbott
Restrictions inherent in this procedure. This method for finding transfer orbits will treat only ellipses and hyperbolas. Circular and parabolic transfer orbits will not be considered. The default coordinate system will be heliocentric ecliptic coordinates in which the XY plane (the ecliptic) shall be the current plane of Earth's orbit, with positions being relative to the center of the sun, and with velocities being sun-relative, unless otherwise noted. All transfer orbits sought by this method will have an apside at one, but not at both, endpoints of the intended trajectory. That is, either departure or arrival will occur at either aphelion or perihelion of the transfer orbit (and both of those instances of "or" are really "xor" - exclusive or). Jerry Abbott
Initial data. The unfortunate fact to someone wanting to figure a transfer orbit is being unable to know, ahead of time, where the arrival position is. What you can actually see is where the spaceship and the rendezvous object are at the moment of departure. And without knowing in advance where the arrival will be, you must do some trial and error work in order to ensure that the spaceship and the rendezvous object will reach the arrival position in the same amount of transit time. But while developing the theory of transfer orbit calculation, it is getting ahead of things to worry about transit time matching. We will do the trial-and-error work later, once the theory is in place for determining the transfer orbit's elements and the transit time in the transfer orbit. For now, we will assume that we know where the arrival will occur. Let these be the orbital elements of Earth. a = 1.00000011 AU e = 0.01671022 i = 0 L = 0 w = 102.94719 degrees T = JD 2453009.3 Let these be the orbital elements of Vesta. a = 2.3626478 AU e = 0.08887781 i = 7.13485 degrees L = 103.94712 degrees w = 149.67895 degrees T = JD 2452941.1 We assert that a spaceship will depart from Vesta at t1 = 4 February 2004, 19h 12m UT t1 = JD 2453040.3 We assert that the spaceship will arrive at Earth at t2 = 16 September 2004, 21h 36m UT t2 = JD 2453265.4 We thus have established a required transit time for the spaceship, going from Vesta to Earth, of t2-t1, or 225.1 days. Our job is to find out whether these assertions can possibly be true, astrodynamically speaking. That is, does there exist a valid transfer orbit from the departure position at t1 to the arrival position at t2? Jerry Abbott
Preburn and Rendezvous State Vectors. Given its orbital elements and a specified moment of time, you can find the heliocentric position and sun-relative velocity at that moment. Find the orbital period. P = (365.256898326 days) a^1.5 Find the mean motion. m = 2 pi radians / P Find the mean anomaly. M = m (t-T) Find the eccentric anomaly. (Danby's Method used.) u(0) = M Repeat over subscript j, f0 = u(j) - e sin u(j) - M f1 = 1 - e cos u(j) f2 = e sin u(j) f3 = e cos u(j) d1 = -f0/f1 d2 = -f0/(f1 + d1 f2 / 2) d3 = -f0/(f1 + d1 f2 / 2 + d2^2 f3 / 6) u(j+1) = u(j) + d3 Until | u(j+1) - u(j) | < 1E-12 The eccentric anomaly, u, is the converged value from the loop. Finding the canonical position vector. x''' = a [ (cos u) - e ] y''' = a (1 - e^2)^0.5 sin u z''' = 0 Rotation by the argument of the perihelion. x'' = x''' cos w - y''' sin w y'' = x''' sin w + y''' cos w z'' = z''' = 0 Rotation by the inclination. x' = x'' y' = y'' cos i z' = y'' sin i Rotation by the longitude of ascending node. x = x' cos L - y' sin L y = x' sin L + y' cos L z = z' This vector [x,y,z] is the heliocentric position vector at time [t] of a planet having orbital elements [a , e , i , L , w , T]. Finding the true anomaly. Q = arctan2( y''' , x''' ) The arctan2 function is the two-argument arctangent. It is related to the principal value arctan function as follows: Function arctan2( y , x ); begin if x=0 then begin if y>0 then Q:=+pi/2 if y=0 then Q:=0 if y<0 then Q:=-pi/2 end else begin Q:=arctan(y/x) if x<0 then Q:=Q+pi if x>0 and y<0 then Q:=Q+2 pi end arctan2:=Q end Velocity in an elliptical orbit. You will want to convert the semimajor axis from astronomical units to meters. 1 astronomical unit = 1.49597870691E+11 meters GMsun = 1.32712440018E+20 m^3 sec^-2 VX’’’ = -sin Q { GMsun / [ a (1-e^2) ] }^0.5 VY’’’ = (e + cos Q) { GMsun / [ a (1-e^2) ] }^0.5 VZ’’’ = 0 The coordinate rotations by w, i, and L proceed as with the position vector. Putting in the numbers... I will from now on usually present angles in radians in the interval [0, 2 pi.) Finding the HEC state vector for Vesta on JD 2453040.3 a = 2.3626478 AU e = 0.08887781 i = 0.1245266 radians L = 1.81422 radians w = 2.612391 radians T = JD 2452941.1 P = 1326.468 days m = 0.004736778 radians per day t-T = 99.2 days M = 0.4698881 radians u = 0.5135515 radians x''' = +1.847896 y''' = +1.156106 z''' = 0 x'' = -2.178777 y'' = -0.06506175 z'' = 0 x' = -2.178777 y' = -0.06455795 z' = -0.008080998 x = +0.5877971 y = -2.098983 z = -0.008080998 Q = 0.5590546 radians Vx''' = -10318.3 m/s Vy''' = +18221.6 m/s Vz''' = 0 ...intermediate steps skipped... Vx = +20234.0 m/s Vy = +4724.0 m/s Vz = -2600.6 m/s Finding the HEC state vector for Earth on JD 2453265.4 a = 1.00000011 AU e = 0.01671022 i = 0 L = 0 w = 1.796767 radians T = JD 2453009.3 P = 365.2569 days m = 0.0172021 radians per day t-T = 256.1 days M = 4.405455 radians u = 4.389608 radians x''' = -0.3339517 y''' = -0.9482125 z''' = 0 x'' = +0.9989289 y'' = -0.1130118 z'' = 0 x' = +0.9989289 y' = -0.1130118 z' = 0 x = +0.9989289 y = -0.1130118 z = 0 Q = 4.373764 radians Vx''' = +28097.2 m/s Vy''' = -9397.8 m/s Vz''' = 0 ...intermediate steps skipped... Vx = +2863.6 m/s Vy = +29488.5 m/s Vz = 0 Summary: Preburn and rendezvous state vectors. The HEC position and velocity of Vesta at the moment of departure are: x1 = +0.5877971 y1 = -2.098983 z1 = -0.008080998 Vx(preburn) = +20234.0 m/s Vy(preburn) = +4724.0 m/s Vz(preburn) = -2600.6 m/s The HEC position and velocity of Earth at the moment of arrival are: x2 = +0.9989289 y2 = -0.1130118 z2 = 0 Vx(rendezvous) = +2863.6 m/s Vy(rendezvous) = +29488.5 m/s Vz(rendezvous) = 0 I could have designated the positions with (preburn) or with (rendezvous), too, but they are also equal to the positions of the transit object at departure and at arrival, respectively. For quantities pertaining to a transfer orbit, departure will usually be denoted by the subscript 1, while arrival will usually be denoted by the subscript 2. Jerry Abbott
Inclination of the transfer orbit. We have two HEC position vectors, R1 and R2, which we suppose to be in the plane of the transfer orbit. The vectors must be distinct and not exactly opposite in direction. We find the cross product, R1xR2. Xn' = y1 z2 - z1 y2 Yn' = z1 x2 - x1 z2 Zn' = x1 y2 - y1 x2 We find the magnitude of the primed normal vector, Rn'. Rn' = [ (Xn')^2 + (Yn')^2 + (Zn')^2 ]^0.5 We divide the primed vector's components each by the magnitude of the primed vector. Xn = Xn' / Rn' Yn = Yn' / Rn' Zn = Zn' / Rn' The vector [Xn, Yn, Zn] is a unit normal vector with respect to the transfer orbit. i = pi/2 - arcsin(Zn) The variable i is the inclination of the transfer orbit plane with respect to the ecliptic. Plugging in the numbers... We suppose that these two HEC position vectors are in the transfer orbit plane. x1 = +0.5877971 y1 = -2.098983 z1 = -0.008080998 x2 = +0.9989289 y2 = -0.1130118 z2 = 0 The cross product gives Xn' = -9.132484E-4 Yn' = -8.072343E-3 Zn' = +2.030307 The magnitude, Rn', is 2.030324 AU The unit normal to the transfer orbit is Xn = -4.498044E-4 Yn = -3.975890E-3 Zn = +0.9999920 The inclination of the transfer orbit, i = 0.003996730 radians Jerry Abbott
Eccentricity and semimajor axis of the transfer orbit. We find the lengths of the three sides of a triangle having the sun, the point of departure, and the point of arrival as their vertices. r1 = [ (x1)^2 + (y1)^2 + (z1)^2 ]^0.5 r1 is the distance from the sun to the point of departure. r2 = [ (x2)^2 + (y2)^2 + (z2)^2 ]^0.5 r2 is the distance from the sun to the point of arrival. d = [ (x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2 ]^0.5 d is the line-of-light separation between the points of departure and arrival. In order to get the eccentricity and semimajor axis of the transfer orbit from these three values, we must restrict our allowed transfer orbits to those having an apside at either endpoint of the intended trajectory. That is, perihelion or aphelion must occur at departure or at arrival. (Both instances of "or" in the previous sentence are exclusive. If one endpoint of the intended trajectory has an apside, the other endpoint may not have the other one, or this procedure will fail.) We proceed by considering, first, the Law of Cosines, understanding that the angle between r1 and r2 is the difference in true anomaly Q2-Q1. cos(Q2-Q1) = ( r1^2 + r2^2 - d^2 ) / (2 r1 r2) Secondly, we solve the polar equation of the orbit, with the pole at the solar focus, for the cosine of the true anomaly. cos Q = [ a (1 - e^2) - r ] / (e r) There are four general cases for a subsequent solution: 1. Departure occurs from the perihelion of the transfer orbit. Outward bound trajectory, Q1 = 0. 2. Arrival occurs at the aphelion of the transfer orbit. Outward bound trajectory, Q2 = pi radians. 3. Arrival occurs at the perihelion of the transfer orbit. Inward bound trajectory, Q2 = 0. 4. Departure occurs from the aphelion of the transfer orbit. Inward bound trajectory, Q1 = pi radians. Now, I'm not about to solve all four cases. I'm just going to solve case #1 explicitly, to provide proof of principle. You can solve the other three cases yourself, if you want to check my work. I'm only going to provide the results for cases #2, #3, and #4. INSIST Q1=0. Therefore, r1 = a(1-e). Why? Because Q1=0 means that departure occurs at the transfer orbit's perihelion, and the heliocentric distance at perihelion equals the product of the semimajor axis and of one minus the eccentricity. We thus have two equations for cos Q2, namely, cos Q2 = ( r1^2 + r2^2 - d^2 ) / (2 r1 r2) cos Q2 = [ r1 (1-e^2) / (1-e) - r2 ] / (e r2) Meaning that... e r2 ( r1^2 + r2^2 - d^2 ) = 2 r1^2 r2 (1+e) - 2 r1 r2^2 ...solve, solve, solve... e = 2 r1 (r1 - r2) / (r2^2 - r1^2 - d^2) a = r1 / (1-e) How easy that was. Case 1. Perihelion at Departure: Q1 = 0. e = 2 r1 (r1 - r2) / ( r2^2 - r1^2 - d^2 ) If e is in [0,1) then an elliptical transfer orbit exists, and a = r1 / (1-e) If e>1 then a hyperbolic transfer orbit exists, and a = r1 / (e-1) If e<0 then there is no transfer orbit having perihelion at departure. Case 2. Aphelion at Arrival: Q2 = pi radians. e = 2 r2 (r1 - r2) / ( r1^2 - r2^2 - d^2 ) If e is in [0,1) then an elliptical transfer orbit exists, and a = r2 / (1+e) If e is not in [0,1) then there is no transfer orbit having aphelion at arrival. Case 3. Perihelion at Arrival: Q2 = 0. e = 2 r2 (r2 - r1) / ( r1^2 - r2^2 - d^2 ) If e is in [0,1) then an elliptical transfer orbit exists, and a = r2 / (1-e) If e>1 then a hyperbolic transfer orbit exists, and a = r2 / (e-1) If e<0 then there is no transfer orbit having perihelion at arrival. Case 4. Aphelion at Departure: Q1 = pi radians. e = 2 r1 (r2 - r1) / ( r2^2 - r1^2 - d^2 ) If e is in [0,1) then an elliptical transfer orbit exists, and a = r1 / (1+e) If e is not in [0,1) then no transfer orbit exists having aphelion at departure. Since our example problem involves a transfer from Vesta to Earth, the trajectory will be inward bound and the transfer orbit we seek will be found in either Case #3 or Case #4. Plugging in the numbers... r1 = 2.179748 r2 = 1.005301 d = 2.028097 Case 1: e = -0.6519092 (unacceptable) Case 2: e = -6.339075 (unacceptable) Case 3: e = 6.339075 (hyperbolic orbit), a = 0.1882913 Case 4: e = 0.6519092 (elliptical orbit), a = 1.319533 So we shall propose that two orbits be investigated as possible transfer orbits between r1 at the time of departure and r2 at the time of arrival. Proposed hyperbolic transfer orbit. a = 0.1882913 e = 6.339075 i = 0.003996730 radians L = ? w = ? T = ? Proposed elliptical transfer orbit. a = 1.319533 e = 0.6519092 i = 0.003996730 radians L = ? w = ? T = ? Jerry Abbott
Velocity in transfer orbit at apsidal end of intended trajectory. I define a subscript variable K to designate the apsidal endpoint of the trajectory. When K=1, the apside is at departure. When K=2, the apside is at arrival. We must determine the sun-relative velocity of the spaceship at the apsidal end of its intended trajectory. Although it would be possible to use the general method (shown for ellipses in an earlier post), there is another way to proceed in this special circumstance. We've found the heliocentric position vector to the apsidal endpoint of the intended trajectory, rK = [xK , yK, zK] and a unit vector normal to the transfer orbit Rn = [Xn , Yn, Zn]. We will obtain the cross product Rn x rK. VxK'' = Yn zK - Zn yK VyK'' = Zn xK - Xn zK VzK'' = Xn yK - Yn xK We find the magnitude of this double-primed vector. VK'' = [ (VxK'')^2 + (VyK'')^2 + (VzK'')^2 ]^0.5 We divide the components of the double-primed vector by the magnitude of the double-primed vector, obtaining a unit vector in the direction of the velocity in the transfer orbit at the apsidal endpoint. VxK' = VxK'' / VK'' VyK' = VyK'' / VK'' VzK' = VzK'' / VK'' We refer to the Vis Viva equation to get the sun-relative speed in the transfer orbit at the apsidal endpoint. To keep the units consistent, you may want to convert [a] and [rK] from astronomical units to meters. 1 astronomical unit = 1.49597870691E+11 meters GMsun = 1.32712440018E+20 m^3 sec^-2 Elliptical orbits: VK = [ GMsun ( 2/rK - 1/a ) ]^0.5 Hyperbolic orbits: VK = [ GMsun ( 2/rK + 1/a ) ]^0.5 VxK = VK VxK' VyK = VK VyK' VzK = VK VzK' The vector VK = [VxK , VyK , VzK] is the velocity in the transfer orbit at the apsidal endpoint of the intended trajectory, referred to heliocentric ecliptic coordinates. This method only works at apsides, not in general. It works at apsides because the heliocentric position vector and the sun-relative velocity are perpendicular to each other there. Plugging in the numbers... Xn = -4.498044E-4 Yn = -3.975890E-3 Zn = +0.9999920 The proposed hyperbolic transfer orbit. x2 = +0.9989326 y2 = -0.1129745 z2 = 0 Vx2'' = +0.1129736 Vy2'' = +0.9989246 Vz2'' = +4.022925E-3 V2'' = 1.005301 Vx2' = +0.1123779 Vy2' = +0.9936575 Vz2' = +4.001713E-3 a = 0.1883722 AU r2 = 1.005301 AU V2 = 80463.3 m/s Vx2 = +9042.3 m/s Vy2 = +79953.0 m/s Vz2 = +322.0 m/s The proposed elliptical transfer orbit. x1 = +0.5877971 y1 = -2.098983 z1 = -0.008080998 Vx1'' = +2.098999 Vy1'' = +0.5877888 Vz1'' = +3.281149E-3 V1'' = 2.179748 Vx1' = +0.9629546 Vy1' = +0.2696590 Vz1' = +1.505288E-3 a = 1.319533 AU r1 = 2.179748 AU V1 = 11902.5 m/s Vx1 = +11461.5 m/s Vy1 = +3209.6 m/s Vz1 = +17.9 m/s The significance of the work so far. We have solved for three of the orbital elements of the transfer orbit (i, e, a), and we have determined a state vector for the transfer orbit at the apsidal endpoint of its trajectory. State vector for the proposed hyperbolic transfer orbit. x2 = +0.9989289 y2 = -0.1130118 z2 = 0 Vx2 = +9042.3 m/s Vy2 = +79953.0 m/s Vz2 = +322.0 m/s State vector for the proposed elliptical transfer orbit. x1 = +0.5877971 y1 = -2.098983 z1 = -0.008080998 Vx1 = +11461.5 m/s Vy1 = +3209.6 m/s Vz1 = +17.9 m/s As yet, we know nothing about the velocity at the non-apsidal endpoint of the intended trajectory. But we shall soon fix that! Jerry Abbott
Longitude of the ascending mode of the transfer orbit. The angular momentum per unit mass, h, is equal to the cross product of the heliocentric radius vector and the velocity vector at the apsidal endpoint: rK x VK. Angular momentum is a conserved quantity, and the components of h are constant for the entire transfer orbit. You will want to convert xK, yK, zK from astronomical units to meters. hX = yK VzK - zK VyK hY = zK VxK - xK VzK hZ = xK VyK - yK VxK h = [ (hX)^2 + (hY)^2 + (hZ)^2 ]^0.5 The longitude of the ascending node, L, of the transfer orbit can be calculated as follows: L = arctan2(hX , -hY) Plugging in the numbers... The proposed hyperbolic orbit. hX = -5.441887E+12 m^2 sec^-1 hY = -4.811776E+13 m^2 sec^-1 hZ = +1.210085E+16 m^2 sec^-1 h = 1.2100947E+16 m^2 sec^-1 L = 6.170569 radians The proposed elliptical orbit. hX = -1.745789+12 m^2 sec^-1 hY = -1.543129E+13 m^2 sec^-1 hZ = +3.881186E+15 m^2 sec^-1 h = 3.881217E+15 m^2 sec^-1 L = 6.170532 radians The longitude of the ascending node had better be the same for both the proposed hyperbolic orbit and the proposed elliptical orbit! (That's because they share the heliocentric position vectors for arrival and departure and therefore occur in the same plane.) Jerry Abbott
Argument of the perihelion of the transfer orbit. cos(w+QK) = (xK cos L + yK sin L) / rK if (i = 0) or (i = pi radians), then sin(w+QK) = (yK cos L - xK sin L) / rK if (i<>0) and (i<>pi radians), then sin(w+QK) = zK / (rK sin i) w = arctan2[ sin(w+QK) , cos(w+QK) ] - QK If necessary, adjust w to the interval [0, 2 pi). Plugging in the numbers... The proposed hyperbolic transfer orbit. x2 = +0.9989326 y2 = -0.1129745 z2 = 0 r2 = 1.005301 AU i = 0.003996730 radians L = 6.170569 radians Q2 = 0 (This orbit's perihelion occurs at arrival at a point on Earth's orbit.) w = 0 Yes, the argument of the perihelion of the proposed hyperbolic transfer orbit is zero. Why? Because the arrival occurs at a point which is both perihelion and on the ecliptic plane, where the Earth is to be found. This node is the ascending node because Vz2 is positive. Hence, the argument of the perihelion, for this orbit, is zero. The proposed elliptical transfer orbit. x1 = +0.5877971 y1 = -2.098983 z1 = -0.008080998 r1 = 2.179748 AU i = 0.003996730 radians L = 6.170532 radians Q1 = pi radians (This orbit's aphelion is at departure from Vesta.) w = 1.9560970 radians Jerry Abbott
True anomaly in the transfer orbit of the non-apsidal end of the intended trajectory. I define a subscript variable J to designate the non-apsidal endpoint of the intended trajectory. If the apside is at departure (K=1), then J=2. If the apside is at arrival (K=2), then J=1. We want to find the heliocentric vector to the non-apsidal endpoint in the coordinate system in which the XY plane contains the transfer orbit, with the +X axis extended through the transfer orbit's perihelion. Once we have this canonical vector, we can use it to get the true anomaly at the non-apsidal endpoint of the intended trajectory. xJ’ = yJ sin L + xJ cos L yJ’ = yJ cos L - xJ sin L zJ’ = zJ xJ’’ = xJ’ yJ’’ = zJ’ sin i + yJ’ cos i zJ’’ = zJ’ cos i - yJ’ sin i xJ’’’ = yJ’’ sin w + xJ’’ cos w yJ’’’ = yJ’’ cos w - xJ’’ sin w zJ’’’ = zJ’’ This triple-primed position vector is the one we were looking for. Important check: Within a reasonable allowance for roundoff error, the value of zJ’’’ should be zero. QJ = arctan2 (yJ’’’ , xJ’’’) Finding the true anomaly of the non-apsidal endpoint of the intended trajectory is a milestone in solving the transfer orbit problem. Plugging in the numbers... The proposed hyperbolic transfer orbit. x1 = +0.5878052 y1 = -2.098982 z1 = -8.082041E-3 L = 6.170569 radians i = 3.99673E-3 radians w = 0 ...intermediate steps skipped... x1''' = +0.8199621 y1''' = -2.019646 z1''' = -1.008436E-5 Q1 = 5.098051 radians The proposed elliptical transfer orbit. x2 = +0.9989289 y2 = -0.1130118 z2 = 0 L = 6.170532 radians i = 3.99673E-3 radians w = 1.9560970 radians ...intermediate steps skipped... x2''' = -0.3778306 y2''' = -0.9315980 z2''' = +8.933364E-11 Q2 = 4.327089 radians Jerry Abbott
Eccentric and mean anomalies in the transfer orbit at the non-apsidal endpoint. For hyperbolic orbits, we proceed as follows: B = cosh uJ = (1/e) (1 + rJ /a) uJ’ = ln [ B + (B^2 - 1)^0.5 ] = arc-cosh(B) (Depending on whether your calculator has an arc-cosh button on it.) If sin QJ < 0 then uJ = -uJ’ else uJ = uJ’ MJ = e sinh uJ - uJ For elliptical orbits, we proceed as follows: sin uJ = (rJ /a) sin QJ / (1-e^2)^0.5 cos uJ = e + (rJ /a) cos QJ uJ = arctan2(sin uJ , cos uJ) MJ = uJ - e sin uJ The proposed hyperbolic transfer orbit. e = 6.33678 a = 0.1883722 AU r1 = 2.179749 AU Q1 = 5.098051 radians u1' = 1.307609 radians sin Q1 is negative. u1 = -1.307609 radians M1 = -9.550008 radians The proposed elliptical transfer orbit. e = 0.6519092 a = 1.319533 r2 = 1.005301 Q2 = 4.327089 radians sin u2 = -0.9310417 cos u2 = +0.3655728 u2 = 5.086543 radians M2 = 5.693351 radians Jerry Abbott
The transit time in the transfer orbit. The proposed hyperbolic transfer orbit. GMsun = 1.32712440018E+20 m^3 sec^-2 1 astronomical unit = 1.49597870691E+11 meters You will want to convert the semimajor axis [a] from astronomical units to meters. m = (86400 sec/day) ( GMsun / a^3 )^0.5 Where m is the hyperbolic mean motion in radians per day. Perihelion at Departure: K=1, J=2. dt = M2 / m, if Q1=0 (thus M1=0) Perihelion at Arrival: K=2, J=1. dt = -M1 / m, if Q2=0 (thus M2=0) In our example, we were investigating a hyperbolic orbit connecting the positions of Vesta at the time of departure and of Earth at the time of arrival. Since the orbit is inward bound (arrival is closer to the sun than departure), the perihelion was made to coincide with arrival. Thus Q2=0. M1 = -9.550008 radians m = 0.2104051 radians per day dt = -M1/m = 45.4 days The proposed elliptical transfer orbit. P = (365.256898326 days) a^1.5 m = 2 pi radians / P Where P is the orbital period in days, and m is the mean motion in radians per day. Apside at Departure: K=1, J=2. dt = M2 / m, if Q1=0 (thus M1=0) dt = (M2-pi) / m, if Q1=pi (thus M1=pi) Apside at Arrival: K=2, J=1. dt = (2 pi - M1) / m, if Q2=0 (thus M2=0) dt = (pi - M1) / m , if Q2=pi (thus M2=pi) In our example, we were investigating an elliptical orbit connecting the positions of Vesta at the time of departure and of Earth at the time of arrival. The orbit had its aphelion at the point of departure. Thus Q1=pi radians. P = 553.6415 days m = 0.01134884 radians per day M2 = 5.693351 radians dt = (M2-pi) / m = 224.85 days Significance of our work so far. In an earlier post, I wrote: Obviously the proposed hyperbolic orbit fails the transit time test and is therefore not a valid transfer orbit. A spaceship leaving Vesta at the time of departure will reach the point of arrival long before Earth does, which would be rather foolish. The spaceship pilot would lose his job. However, the proposed elliptical orbit is a near-match (six hours difference). Some adjustment in the time of arrival will probably result in a perfect match. The only thing left to do is calculate the delta-vees. Jerry Abbott
The delta-vees for departure and for arrival. The canonical velocity in a hyperbolic orbit, in general, is found from Vx’’’ = -(a/r) { GMsun / a }^0.5 sinh u Vy’’’ = +(a/r) { GMsun / a }^0.5 (e^2 - 1)^0.5 cosh u Vz’’’ = 0 The canonical velocity in an elliptical orbit, in general, is found from Vx’’’ = -sin Q { GMsun / [ a (1-e^2) ] }^0.5 Vy’’’ = (e + cos Q) { GMsun / [ a (1-e^2) ] }^0.5 Vz’’’ = 0 These triple-primed vectors would be rotated (negatively) by the angular elements of the orbit (w,i,L) to heliocentric ecliptic coordinates. Rotation by the argument of the perihelion, w. Vx'' = Vx''' cos w - Vy''' sin w Vy'' = Vx''' sin w + Vy''' cos w Vz'' = Vz''' = 0 Rotation by the inclination, i. Vx' = Vx'' Vy' = Vy'' cos i Vz' = Vy'' sin i Rotation by the longitude of ascending node, L. Vx = Vx' cos L - Vy' sin L Vy = Vx' sin L + Vy' cos L Vz = Vz' By comparing transit times for Earth and for the spaceship from their respective starting points (at the time of departure) to the arrival location, we determined that the elliptical orbit in our example problem was a valid (or almost so) transfer orbit from Vesta to Earth. That orbit had aphelion at the departure from Vesta. Using a special method (involving a cross product and the Vis Viva equation) that works only at apsides, we determined the HEC velocity in the transfer orbit at departure: [Vx1, Vy1, Vz1]. Using the general method for finding velocity in an elliptical orbit, we obtained Vesta's own HEC velocity at departure: [ Vx(preburn) , Vy(preburn) , Vz(preburn) ] The delta-vee for departure is dV1x = Vx1 - Vx(preburn) dV1y = Vy1 - Vy(preburn) dV1z = Vz1 - Vz(preburn) We may now use the general method for finding velocity in an elliptical orbit to determine the velocity in the transfer orbit at the non-apsidal endpoint of the intended trajectory (i.e., at Earth). The result will be: [Vx2, Vy2, Vz2]. We may use the general method, also, to determine the HEC velocity of Earth in its own orbit at the moment of arrival, namely: [ Vx(rendezvous) , Vy(rendezvous) , Vz(rendezvous) ] The delta-vee required of the spaceship to match speed with Earth is dV2x = Vx(rendezvous) - Vx2 dV2y = Vy(rendezvous) - Vy2 dV2z = Vz(rendezvous) - Vz2 Since the spaceship's transit time wasn't exactly a match for the Earth's, I won't plug in the numbers on this occasion. These equations do, however, show the way it is done. UPDATE: See a later post for the delta-vees for departure and arrival for the elliptical transfer orbit. Jerry Abbott
Long-path and superperiodic elliptical transfer orbits. The elliptical transfer orbits treated in the preceding posts are what I've come to call "elliptical transfer orbits of the short path," meaning that the transit object moves through an arc of true anomaly less than pi radians. For each elliptical transfer orbit of the short path, there exists an elliptical orbit of the long path, in which the transit object moves along the same orbit in the opposite direction. Long path elliptical orbits differ from short path elliptical orbits in the following ways. At each point where the HEC position vectors are identical, the HEC velocity vectors differ (only) in sign. The sum of the transit times for the long and short paths is equal to the period of the transfer orbit. The inclinations sum to pi radians. The angular momentum vector of the long path orbit is the reverse of that of the short path. The longitude of the ascending node on the long path will be that of the short path plus or minus pi radians, whichever will keep it inside the interval [0 , 2 pi). The sum of the short-path and long-path arguments of the perihelion should add up to pi radians +/- some whole multiple of 2 pi radians. The sum of the non-apsidal endpoint true anomalies for the short path elliptical orbit and for the long-path elliptical orbit is 2 pi. And the like. There is yet another class of intercept orbit, which might be called a super-periodic transfer orbit. They may be either prograde or retrograde, as we normally reckon such, and involve transit times exceeding the transfer orbit's period. Such extra possible orbits gives a spaceship pilot more opportunities to find a valid transfer orbit, however at the cost of increasing his transit time, often greatly. Jerry Abbott
Comments? Applause? Honorary master's degree? We have now reached the end of this dissertation. I will entertain questions and engage in a certain amount of pedagogical banter. Jerry Abbott
Matching up transit times. Once you have the above procedure coded, you can use trial-and-error to find valid transfer orbits. You'd choose a time for departure (t1) and a time for arrival (t2), thereby fixing the positions of the rendezvous planet at those times, as well as fixing the transit time for the rendezvous planet (t2-t1). Then you crank through the procedure to determine whether (or not) the proposed transfer orbit connecting the location of the spaceship at t1 with the rendesvous planet's location at t2 actually does occur in the same time difference. If the transit time along the proposed transfer orbit isn't even close, discard that proposed transfer orbit: it's no good. If the transit time along the proposed transfer orbit is close, you'd play around with different departure or arrival times until you have an exact transit time match with both the rendezvous planet and the spaceship. Now, there's another way to proceed, if you don't mind the work. For any two objects - Vesta and Earth in our example - it would be possible to construct a "Big Book of Object 1 to Object 2 Transfer Orbits" comprised mostly of a table indexed by the heliocentric longitudes of the two objects. The index interval might be tenths of a degree, or milliradians, or something convenient to spaceship pilots. This table is the result of some reverse interpolation. To construct it to begin with, you need to use the mean anomalies of Object 1 and Object 2 as indices going forward, with the transfer orbit(s) elements, transit times, and the heliocentric longitudes of both planets at the time of departure being outputs. Finding the heliocentric longitude of the preburn position is relatively simple: one use of the arctan2 function on y1,x1. The rendezvous planet must be backtracked in its own orbit by the transit time in the transfer orbit. To do this, a reduction of the rendezvous planet's orbital elements is conducted with the "time of interest" being t2-dt. The resulting HEC position vector might be labeled [xbt, ybt, zbt]. The heliocentric longitude of the rendezvous planet at the time of departure is arctan2(ybt,xbt). The output is only made in the event that t2-dt = t1, to within some appropriately chosen error tolerance. Otherwise, the proposed transfer orbit in that case would be found to be invalid and accordingly of no interest. Having found the heliocentric longitudes of both planets at the moment of departure, one now has the indices needed to construct the reverse interpolation table from which the sampled set of valid transfer orbits may be read off as a first rough approximation for spaceship pilots wanting to go from a point in one orbit to rendezvous with a planet in a different orbit. Jerry Abbott
Call me impressed! One topic that came up here some time ago was whether the advantage to burning as much of a rocket's propellant while low in a gravitational field applies only when the field is considered as varying in strength (as when a probe flies by a planet), or if it also applies to modest rockets going straight up to an apogee and then falling straight back down in a field that is considered constant, e.g. a flight of a couple of miles up from Earth's surface. A key assumption was made that the beefing-up of the structure to allow for higher acceleration would be ignored. So you've got two rockets to compare, both with same structure weight and same propellant weight and same total impulse. Your thoughts?
It sounds like the first rocket is taking advantage of the gravity well of a planet by holding on to its fuel's mass until most of its potential energy has become kinetic energy, then blowing the fuel away contrary to the spaceship's velocity, which transfers some (how much depends on the exhaust speed) of the fuel mass's kinetic energy into the spaceship. The ship then rises out of the gravity well still having all of the energy from the fall, the energy from the thrust, plus the gravitational potential energy stolen from the fuel. For the second spaceship, there's also an advantage to losing the weight of the propellant as fast as possible. The least efficient (zero efficiency) launching rocket is one that can do no more than hover by making thrust equal only to its own weight. When launching from a planet with an atmosphere, though, there may be a tradeoff between thrust and aerodynamic pressure, so that you'd want to get above the soupy part of the atmosphere before flooring the pedal. I'm really not much of an expert in astrodynamics. I really am a hillbilly. I live in the Alleghenies near Hillsboro, West Virginia, where I raise goats and grow potatoes and apple and walnut trees. Now and then I like to figure things out, like I did with transfer orbits a few years ago. Jerry Abbott
Depends on what you want I guess. If you want to gain speed using planetary a flyby, it's more advantageous to fire at pericentre. This is because speed is largest there. The speed that you add by burning fuel will always be the same (it depends on the fuel, engine, etc). So, because kinetic energy goes with the square of the velocity, adding that same velocity to a large speed will give you a larger increase in kinetic energy (and because potential energy stays constant during the burn, the total energy increase is the same as the kinetic energy increase for both cases.)
Thanks to both of you for your responses. For the case of the modest rocket going straight up and straight down in a constant-g field... If we ignore the effect of an atmosphere, we ought to be able to graph apogee vs. thrust and get some sort of curve, with fixed parameters being structure mass, propellant mass, and total impulse, and the constraint that thrust times burn time equals the total impulse. As Jenab noted, "The least efficient (zero efficiency) launching rocket is one that can do no more than hover by making thrust equal only to its own weight." To my intuition, that means the curve will not have its peak over toward the left side of the graph where the thrust is low. In fact, my intuition suggests the curve is monotonically increasing with increasing thrust. (In the real world, structure mass would have to increase when thrust gets bigger to keep the rocket from collapsing or blowing up, so that its non-constancy would have to be taken into account, so the real-world version of the graph would actually have some maximum apogee for some finite value of thrust, beyond which increasing thrust lowers the apogee; "the point of diminishing returns" you might call it.) Remcook, when you say, "because potential energy stays constant during the burn," are you making the approximation that the burn time is short compared to the time it takes for the gravitational potential to change appreciably when the probe is near its minimum point in the potential field? And another question for Jenab: can a West Virginian still grab onto a slow-moving coal hopper and go places in your state? Or are those days long gone? As it happens, I recently read an article on M&K Junction in West Virginia, across the Cheat River from Rowlesburg.