Jenab Posted May 21, 2004 Posted May 21, 2004 Okay fellows. Here's my challenge. There existed a transfer orbit from asteroid Vesta (IAU number 4) to Earth with a departure date of 4 February 2004. (Go look up the orbital elements of Earth and Vesta; that's part of the work!) Find: The classical orbital elements of the transfer orbit, the transit time, and the change-of-velocity vectors at departure and at arrival. The classical orbital elements are: semimajor axis, eccentricity, inclination to the ecliptic, longitude of the ascending node, argument of perihelion, and Julian date of perihelion passage. The delta-vee vectors should be referred to heliocentric ecliptic coordinates. Please show your enough of your work to prove that you know what you're doing. High praise will reward the poster who provides the correct answers first! Good luck to all who try. Jerry Abbott
Jenab Posted May 22, 2004 Author Posted May 22, 2004 Departure from Vesta: 4 February 2004, 19h 12m UT JD 2453040.3 Arrival at Earth. 16 September 2004, 21h 36m UT JD 2453265.4 The transit time is therefore 225.1 days. The point of departure in heliocentric ecliptic coordinates can be calculated from a reduction of Vesta's orbital elements and the time of departure. The point of arrival in heliocentric ecliptic coordinates can be calculated from a reduction of Earth's orbital elements and the time of arrival. With these points (the endpoints of the intended trajectory) both known, it becomes much simpler to solve for the elements of the transfer orbit. Had only the departure time and position been known, it would have been necessary to construct a large table of TRIAL transfer orbits from Vesta to various points on Earth's orbit, and checking each of them to see whether Earth and the spaceship had the same transit time. In general, they would not; i.e., most of the trial transfer orbits would not be actual transfer orbits between the two planets. Only when you discovered a match between transit times would you have found a valid transfer orbit. Now that you have the arrival time, you can straightforwardly proceed to solve for the elements of the transfer orbit. If nobody ventures a serious attempt at a solution in a few days, I'll begin solving the problem myself. Jerry Abbott
Jenab Posted May 24, 2004 Author Posted May 24, 2004 My guess is that not many people are interested in transfer orbits, and those who are intrigued by this challenge are daunted by what they think is a difficult task. But, really, if you've had vector algebra in high school and differential calculus and analytic geometry in college, you should be able to figure this out. And if you've had an astronomy course in celestial mechanics, the problem should be a cakewalk. But, okay, for those who need a bit of jumpstarting... The first thing to do is finding the position vector of departure and the position vector of arrival. That is, you want to calculate: 1. The position vector of Vesta, R1, at the time of departure, which is 4 February 2004, 19h 12m UT (or JD 2453040.3). 2. The position vector of Earth, R2, at the time of arrival, which is 16 September 2004, 21h 36m UT (or JD 2453265.4). Both position vectors should be referred to heliocentric ecliptic coordinates. Once these vectors are known, a cross product (R1xR2) will give you the vector normal to the transfer orbit, which should be divided by its own magnitude to result in a unit normal to the transfer orbit. The inclination of the transfer orbit is found from the Z component of this unit normal vector. The next elements to tackle are the eccentricity and the semimajor axis. Use the Pythagorean Theorem to get the line-of-light distance, d, between the point of departure and the point of arrival. Now comes a neat trick. To find the eccentricity of the transfer orbit, solve the polar equation of the orbit (with the pole at the solar focus) simultaneously with the law of cosines. After a bit of algebra, a set of equations for the eccentricity, e, pops right out. After which, solving for the semimajor axis, a, is a trivial one-step process. Jerry Abbott
Tesseract Posted May 24, 2004 Posted May 24, 2004 You posted three times in a row, i would think that nobody wants to answer the question.Try a more interesting thread next time.
jordan Posted May 24, 2004 Posted May 24, 2004 But, really, if you've had vector algebra in high school and differential calculus and analytic geometry in college, you should be able to figure this out. And if you've had an astronomy course in celestial mechanics, the problem should be a cakewalk. Although I enjoy a good challenge, you just about summed up my problem right there. Sorry, but I have no idea what to do on this one. And don't mind Tesserect's cordial way with words.
Jenab Posted May 24, 2004 Author Posted May 24, 2004 You posted three times in a row, i would think that nobody wants to answer the question.Try a more interesting thread next time. By the time I'm finished in this thread, it will contain a mini-lesson in finding transfer orbits. Somebody might think it worthwhile, even if that somebody has not yet found it. I've had a chance to look around this forum, and I've noticed that a great deal of the discussions are mostly cursory treatments of popular science. I see here many of the same questions about relativity and quantum mechanics that I've seen in many other places, with the same hand-waving and lack of mathematical rigor. It's gee-whiz to the max, but the educational load in the average thread is really fairly low. For once, let's change that. Jerry Abbott
Jenab Posted May 24, 2004 Author Posted May 24, 2004 The orbital elements that I'll use for Vesta, except the mean anomaly, came from http://arnold.usno.navy.mil/murison/Asteroids/OrbitalElements.html Vesta's mean anomaly at epoch was found at http://cfa-www.harvard.edu/iau/Ephemerides/Bright/2004/00004.html mean anomaly at epoch, M, 15.64391 degrees epoch is JD 2453000.5 (27 Dec 2003) Vesta's period is found by raising the semimajor axis to the power of 1.5 and multiplying the result by 365.256898326 days. period, P, 1325.593765 days. mean daily motion, m = 0.27157641 degrees per day Vesta was 57.6040821 days past its perihelion at epoch, therefore... VESTA ORBITAL ELEMENTS time of perihelion passage, T, JD 2452942.9 semimajor axis, a, 2.36160912 astronomical units eccentricity, e, 0.08871313 inclination, i, 7.133087 degrees longitude of ascending node, L, 103.937391 degrees argument of perihelion, w, 150.232154 degrees EARTH ORBITAL ELEMENTS time of perihelion passage, T, JD 2453009.3 semimajor axis, a, 1.00000011 astronomical units eccentricity, e, 0.01671022 inclination, i, defined zero (ecliptic taken as plane of Earth's orbit) longitude of ascending node, L, defined zero argument of perihelion, w, 102.94719 degrees If you want to convert from calendar date to Julian date, go to... http://wwwmacho.mcmaster.ca/JAVA/JD.html If you want to convert from Julian date to calendar date, go to... http://wwwmacho.mcmaster.ca/JAVA/CD.html Dates and times for Earth's perihelions, aphelions, solstices and equinoxes can be found at... http://aa.usno.navy.mil/data/docs/EarthSeasons.html Jerry Abbott
Sayonara Posted May 24, 2004 Posted May 24, 2004 You posted three times in a row, i would think that nobody wants to answer the question.Try a more interesting thread next time. SILENCE! It's an interesting read.
Tesseract Posted May 24, 2004 Posted May 24, 2004 SILENCE! It's an interesting read. The kind of thing that will save the world.If you dont know this, how can you expect to be a valuable member of society?
Sayonara Posted May 24, 2004 Posted May 24, 2004 What are you talking about? In fact no, don't tell me... just stop polluting this thread.
Jenab Posted May 24, 2004 Author Posted May 24, 2004 Given its orbital elements and a specified moment of time, you can find where a planet was at that moment. The procedure for doing that is called a reduction of the orbital elements, and I present two examples below on how to do that. These particular examples were chosen, of course, because they are both part of the larger problem of finding a transfer orbit between the two points in space that will be the results. THE PROCEDURE Find the orbital period. P = (365.256898326 days) a^1.5 Find the mean motion. m = 2 pi radians / P Find the time difference from perihelion. dt = t - T Find the mean anomaly. M = m dt Find the eccentric anomaly. The equation relating the mean anomaly, the orbital eccentricity and the eccentric anomaly is known as Kepler's equation. It is transcendental in the eccentric anomaly, which is why it is necessary to use radian measure for this step. I'll use Danby's method to solve for the eccentric anomaly. Most of the time, the simpler Newton's method will also work, but Danby's method is faster and has fewer convergence failures. 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]. Solving for VESTA (t = JD 2453040.3) a = 2.36160912 astronomical units P = 1325.594 days m = 0.004739903 radians per day dt = 2453040.3 - 2452942.9 = 97.4 days M = 0.4616672 radians e = 0.08871313 u = 0.5045526 radians x''' = +1.857826 y''' = +1.137138 z''' = 0 w = 2.622046 radians x'' = -2.177249 y'' = -0.06469965 z'' = 0 i = 0.1244959 radians x' = -2.177249 y' = -0.06419891 z' = -0.008034048 L = 1.81405 radians x = +0.5867243 y = -2.097687 z = -0.008034048 Solving for EARTH (t = JD 2453265.4) a = 1.00000011 astronomical units P = 365.2569 days m = 0.0172021 radians per day dt = 2453265.4 - 2453009.3 = 256.1 days M = 4.405455 radians e = 0.01671022 u = 4.389608 radians x''' = -0.3339160 y''' = -0.9482244 z''' = 0 w = 1.796767 radians x'' = +0.9989326 y'' = -0.1129745 z'' = 0 i = 0 radians x' = +0.9989326 y' = -0.1129745 z' = 0 L = 0 radians x = +0.9989326 y = -0.1129745 z = 0 Summary. These two points in spacetime exist on the transfer orbit. Departure from Vesta x1 = +0.5867243 y1 = -2.097687 z1 = -0.008034048 t1 = JD 2453040.3 Arrival at Earth x2 = +0.9989326 y2 = -0.1129745 z2 = 0 t2 = JD 2453265.4 The spatial components of these vectors are referred to heliocentric ecliptic coordinates, and the values are in astronomical units. The time component for each vector is referred to Julian Date (JD). The accuracy of the vectors depends on the accuracy with which the orbital elements were known. Possibly, improved values could be found somewhere, but we'll go with these. Jerry Abbott
Jenab Posted May 24, 2004 Author Posted May 24, 2004 Refer to my previous post for nomenclature. We will eventually want to know the velocity of Vesta at departure and the velocity of Earth at arrival, in order to calculate the necessary delta-vee at each end of the intended trajectory. Let Q be the true anomaly of an object in its orbit around the sun. (I don't have access to the greek letters here.) To determine its value: Q = arctan2(y''' , x''') Where x''' and y''' are two of the components of the canonical position vector (see previous post). The arctan2 function is the two-argument arctangent function. It is related to the single-argument 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 GMsun is the solar gravitational constant. GMsun = 1.32712440018E+20 m^3 sec^-2 When using the equations below, make sure to keep the units consistant. You will probably want to enter the semimajor axis in meters, rather than attempt to adjust GMsun. 1 astronomical unit = 1.49597870691E+11 meters. Vx''' = -sin Q { GMsun / [ a (1 - e^2) ] }^0.5 Vy''' = (e + cos Q) { GMsun / [ a (1 - e^2) ] }^0.5 Vz''' = 0 In order to adjust this canonical velocity vector to heliocentric ecliptic coordinates, you'd rotate it in the same way that a position vector would be rotated. Again, see my preceding post. Plugging in the numbers... Orbital velocity of Vesta at JD 2453040.3 x''' = +1.857826 y''' = +1.137138 Q = 0.5492546 radians a = 2.36160912 e = 0.08871313 Canonical velocity Vx''' = -10158.3 m/s Vy''' = +18322.5 m/s Vz''' = 0 ...intermediate steps skipped... HEC velocity Vx = +20241.4 m/s Vy = +4735.7 m/s Vz = -2601.2 m/s Orbital velocity of Earth at JD 2453265.4 x''' = -0.3339160 y''' = -0.9482244 Q = 4.37380123 radians a = 1.00000011 e = 0.01671022 Canonical velocity Vx''' = +28097.7 m/s Vy''' = -9396.8 m/s Vz''' = 0 ...intermediate steps skipped... HEC velocity Vx = +2862.5 m/s Vy = +29488.7 m/s Vz = 0 Jerry Abbott
Jenab Posted May 25, 2004 Author Posted May 25, 2004 Given two heliocentric position vectors to points on an orbit, excepting a pair differing in true anomaly by pi radians, it is possible to determine the inclination of the orbit to the ecliptic. Position vector 1, R1, has these components. x1 = +0.5867243 y1 = -2.097687 z1 = -0.008034048 Position vector 2, R2, has these components. x2 = +0.9989326 y2 = -0.1129745 z2 = 0 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) Where is the inclination of the transfer orbit to the ecliptic plane. Plugging in the numbers... Xn' = -0.0009076426 Yn' = -0.0080254725 Zn' = +2.0291630445 Rn' = 2.02917911803 Xn = -0.0004472955 Yn = -0.0039550340 Zn = +0.9999920788 i = 0.00398025345 radians Jerry Abbott
Tesseract Posted May 25, 2004 Posted May 25, 2004 You should be a professor Jenab you clearly show your talents of slowly explaining things.
Jenab Posted May 26, 2004 Author Posted May 26, 2004 You should be a professor Jenab you clearly show your talents of slowly explaining things. That's what the eMode test said after declaring me a Visionary Philosopher with a 142 IQ. 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.178210435 r2 = 1.005300740 d = 2.027082617 Case #3. e = 6.287121148 a = 0.190141423 A hyperbolic orbit, with perihelion at arrival point. Case #4. e = 0.651493744 a = 1.318933507 An elliptical orbit, with aphelion at departure point. It is not at all likely that both of these orbits are valid transfer orbits. One of them might be (ok, is...I've already checked), but the other will turn out to have an inappropriate transit time. A spaceship using the wrong orbit will reach the "arrival point" either much sooner or much later than the Earth does, which would be kind of foolish. The spaceship pilot will have to choose the right orbit, if he wants to remain employed. Jerry Abbott
Jenab Posted May 26, 2004 Author Posted May 26, 2004 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 in Post #12, there is another way to proceed in this special circumstance. We found a unit vector normal to the transfer orbit Rn = [Xn , Yn, Zn] in Post #13. We also have the heliocentric position vector to the apsidal endpoint, rK = [xK , yK, zK], which was calculated in Post #11. 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. Elliptical orbits: VK = [ GMsun ( 2/rK - 1/a ) ]^0.5 Hyperbolic orbits: VK = [ GMsun ( 2/rK + 1/a ) ]^0.5 Again, remember to keep the units consistent. You may want to convert [a] and [rK] from astronomical units to meters, for example. 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... Possible hyperbolic transfer orbit Perihelion at arrival at Earth on JD 2453265.4 K=2 Xn = -0.0004472955 Yn = -0.0039550340 Zn = +0.9999920788 xK = +0.9989326 yK = -0.1129745 zK = 0 VxK'' = +0.11297361 VyK'' = +0.99892469 VzK'' = +0.00400135 VK'' = 1.00530074 VxK' = +0.11237792 VyK' = +0.99365757 VzK' = +0.00398025 GMsun = 1.32712440018E+20 m^3 sec^-2 1 astronomical unit = 1.49597870691E+11 meters a = 2.8444752E+10 meters rK = 1.50390850E+11 meters VK = 80190.5 m/s VxK = +9011.6 m/s VyK = +79681.9 m/s VzK = +319.2 m/s If this hyperbolic orbit turns out to be an actual transfer orbit - which we don't yet know for certain - we would determine the necessary delta-vee for matching velocity with Earth. Possible elliptical transfer orbit Aphelion at Departure at Vesta on JD 2453040.3 K=1 Xn = -0.0004472955 Yn = -0.0039550340 Zn = +0.9999920788 xK = +0.5867243 yK = -2.097687 zK = -0.008034048 VxK'' = +2.09770216 VyK'' = +0.58671606 VzK'' = +0.00325880 VK'' = 2.17821044 VxK' = +0.96303926 VyK' = +0.26935692 VzK' = +0.00149609 GMsun = 1.32712440018E+20 m^3 sec^-2 1 astronomical unit = 1.49597870691E+11 meters a = 1.97309644E+11 meters rK = 3.25855643E+11 meters VK = 11913.7 m/s VxK = +11473.4 m/s VyK = +3209.1 m/s VzK = +17.8 m/s If this elliptical orbit turns out to be an actual transfer orbit - which we don't yet know for certain - we would determine the necessary delta-vee for entering the transfer orbit from initially being at rest with respect to Vesta. Significance of our progress so far. We now have complete state vectors for two proposed transfer orbits at whichever endpoint of their intended trajectories occurs at an apside: [xK, yK, zK, VxK, VyK, VzK] The proposed hyperbolic transfer orbit has this state vector: x2 = +0.9989326 AU y2 = -0.1129745 AU z2 = 0 Vx2 = +9011.6 m/s Vy2 = +79681.9 m/s Vz2 = +319.2 m/s The proposed elliptical transfer orbit has this state vector: x1 = +0.5867243 AU y1 = -2.097687 AU z1 = -0.008034048 AU Vx1 = +11473.4 m/s Vy1 = +3209.1 m/s Vz1 = +17.8 m/s Jerry Abbott
Jenab Posted May 27, 2004 Author Posted May 27, 2004 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... Proposed hyperbolic transfer orbit hX = -5.39471769E+12 m^2 sec^-2 hY = -4.77006702E+13 m^2 sec^-2 hZ = +1.19075189E+16 m^2 sec^-2 h = 1.19076157E+16 m^2 sec^-2 L = 6.17056860 radians Proposed elliptical transfer orbit hX = -1.72886746E+12 m^2 sec^-2 hY = -1.53519637E+13 m^2 sec^-2 hZ = +3.88213341E+15 m^2 sec^-2 h = 3.88216415E+15 m^2 sec^-2 L = 6.17104239 radians Comment. Strictly, the longitude of the ascending node for both of these proposed transfer orbits should be equal, as they both contain the same two heliocentric radius vectors and therefore occur in the same plane. The most likely reason that they are not quite equal is roundoff error. Their closeness provides a check on my work so far. Jerry Abbott
Jenab Posted May 27, 2004 Author Posted May 27, 2004 The inclination [symbolized as i] of an orbit is the angle between the plane of the orbit and the plane of the ecliptic. Herein the ecliptic plane is assumed to be the plane of Earth's orbit. These two planes intersect in a line called "the line of nodes." The orbit itself crosses the ecliptic at two points, called nodes. At one node, called the ascending node, the planet (or whatever the orbiting object happens to be) crosses the ecliptic with the Z component of its velocity in the direction of the North Ecliptic Pole. At the other node (the descending node), the planet crosses the ecliptic with the Z component of its velocity toward the South Ecliptic Pole. The heliocentric ecliptic coordinate system is defined such that the sun is at the origin, with +X axis extending toward the Vernal Equinox, and +Z axis normal to the Earth's orbit. It is a right-handed system, in which the +Y axis is 90 degrees counterclockwise from the +X axis to an observer positioned somewhere up the +Z axis looking back toward the origin. The longitude of the ascending node [L] is the angle, subtended at the sun (measured in the ecliptic counterclockwise by the aforementioned observer) from the Vernal Equinox to the ascending node. The argument of the perihelion [w] is 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. The true anomaly [Q] is the angle, subtended at the sun, measured in the plane of the orbit in the direction of motion, from the orbit's perihelion to the current position of the planet. 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.00530074 i = 0.00398025345 radians L = 6.17056860 radians cos(w+Q2) = 1.00000000 (to the limit of justified precision) sin(w+Q2) = 0 w + Q2 = 0 This is the Case #3, having perihelion at arrival, so... Q2 = 0 w = 0 The zero result for the argument of the perihelion should not surprise anyone. The arrival occurs at the perihelion of the transfer orbit and also in the plane of Earth's orbit - i.e., in the ecliptic. The ascending node and the perihelion are in the same spot, and thus the argument of the perihelion is zero. The proposed elliptical transfer orbit x1 = +0.5867243 y1 = -2.097687 z1 = -0.008034048 r1 = 2.17821044 i = 0.00398025345 radians L = 6.17104239 radians cos(w+Q1) = +0.37543977 sin(w+Q1) = -0.92666979 w + Q1 = 5.09732665 radians This is the Case #4, having aphelion at departure, so... Q1 = pi radians w = 1.95573400 radians Jerry Abbott
Jenab Posted May 27, 2004 Author Posted May 27, 2004 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 desire a canonical position vector for the non-apsidal endpoint of the intended trajectory. That is, 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. 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.5867243 y1 = -2.097687 z1 = -0.008034048 L = 6.17056860 radians i = 0.00398025345 radians w = 0 x1' = +0.81874324 y1' = -2.01846369 z1' = -0.008034048 x1'' = +0.81874324 y1'' = -2.01847968 z1'' = -8.50786288E-9 Notice how small z1'' is. Actually, it would be zero, but for roundoff error. x1''' = +0.81874324 y1''' = -2.01847968 z1''' = -8.50786288E-9 Since w=0 for this orbit, the triple-primed vector is equal to the double-primed vector. Q1 = 5.09773397 radians The proposed elliptical transfer orbit x2 = +0.9989326 y2 = -0.1129745 z2 = 0 L = 6.17104239 radians i = 0.00398025345 radians w = 1.95573400 radians x2' = +1.00530063 y2' = -4.76296468E-4 z2' = 0 x2'' = +1.00530063 y2'' = -4.76292695E-4 z2'' = +1.89577565E-6 x2''' = -0.37793320 y2''' = -0.93155573 z2''' = +1.89577565E-6 Q2 = 4.32697753 radians Jerry Abbott
Jenab Posted May 27, 2004 Author Posted May 27, 2004 We have reached a point where we must begin treating hyperbolic and elliptical orbits somewhat differently. The proposed hyperbolic transfer orbit B = cosh uJ = (1/e) (1 + rJ /a) uJ’ = ln [ B + (B^2 - 1)^0.5 ] If sin QJ < 0 then uJ = -uJ’ else uJ = uJ’ MJ = e sinh uJ - uJ Here, the uJ is the hyperbolic eccentric anomaly in the transfer orbit of the non-apsidal endpoint of the intended trajectory, and MJ is the corresponding hyperbolic mean anomaly. When calculating MJ it is necessary that uJ be in radians. You will want to convert [a] from astronomical units to meters. GMsun = 1.32712440018E+20 m^3 sec^-2 1 astronomical unit = 1.49597870691E+11 meters m = (86400 sec/day) ( GMsun / a^3 )^0.5 The hyperbolic mean motion will result to units of radians per day. Calculation of the transit time, dt. Once we have the transit time, we can compare it to the time difference required by the known times of departure and of arrival. If they match, or very nearly match, then we have found our transfer orbit. If they do not match, we are disgusted at having done so much work for nothing. Case #1: Perihelion at Departure: K=1, J=2. dt = M2/m, if Q1=0 (thus M1=0) Case #3: Perihelion at Arrival: K=2, J=1. dt = -M1/m, if Q2=0 (thus M2=0) Note that hyperbolic orbits do not have aphelia. The proposed elliptical transfer orbit 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 uJ is the eccentric anomaly in the transfer orbit of the non-apsidal endpoint of the intended trajectory, and MJ is the corresponding mean anomaly. When calculating MJ it is necessary that uJ be in radians. Although you generally don’t intend to follow an elliptical transfer orbit around its complete circuit, its orbital period can be found from P = (365.256898326 days) a^1.5 The mean motion (radians per day) of a spaceship in the transfer orbit would be m = 2 pi radians / P Calculation of the transit time, dt. Once we have the transit time, we can compare it to the time difference required by the known times of departure and of arrival. If they match, or very nearly match, then we have found our transfer orbit. If they do not match, we are disgusted at having done so much work for nothing. Apside at Departure: K=1, J=2. Case 1: dt = M2/m, if Q1=0 (thus M1=0) Case 4: dt = (M2-pi)/m, if Q1=pi (thus M1=pi) Apside at Arrival: K=2, J=1. Case 3: dt = (2 pi-M1)/m, if Q2=0 (thus M2=0) Case 2: dt = (pi-M1)/m, if Q2=pi (thus M2=pi) If necessary, correct dt to the interval [0,P) by adding or subtracting the appropriate multiple of P. Plugging in the numbers... Proposed hyperbolic transfer orbit r1 = 2.17821044 e = 6.287121148 a = 0.190141423 AU Q1 = 5.09773397 radians u1 = -1.30600661 radians M1 = -9.44655271 radians m = 0.20747527 radians per day dt = 45.5 days Whoops! This transit time does not equal the required time difference of 225.1 days! So, although the hyperbolic orbit does connect r1 with r2, a spaceship travelling along it would go from r1 to r2 too quickly, and the destination planet would not have arrived yet. So this orbit won't work. Proposed elliptical transfer orbit r2 = 1.00530074 e = 0.651493744 a = 1.318933507 AU Q2 = 4.32697753 radians sin u2 = -0.53583329 cos u2 = +0.36494920 u2 = 5.31030870 radians M2 = 5.84877378 radians P = 553.26446728 days m = 0.0113565676 radians per day dt = 238.4 days Hm, that's what I get for using orbital elements different from those I used in my last episode of working out this problem. Ok, the spaceship pilot sees that his trip will take a little too long, and so he applies a course correction. Jerry Abbott
Jenab Posted May 28, 2004 Author Posted May 28, 2004 We now find the velocity vector, VJ, in the transfer orbit at the non-apsidal endpoint of the intended trajectory. Canonical velocity at non-apsidal endpoint for hyperbolic orbits VxJ’’’ = -(a/rJ) { GMsun / a }^0.5 sinh uJ VyJ’’’ = +(a/rJ) { GMsun / a }^0.5 (e^2 - 1)^0.5 cosh uJ VzJ’’’ = 0 Canonical velocity at non-apsidal endpoint for elliptical orbits VxJ’’’ = -sin QJ { GMsun / [ a (1-e^2) ] }^0.5 VyJ’’’ = (e + cos QJ) { GMsun / [ a (1-e^2) ] }^0.5 VzJ’’’ = 0 Rotation to heliocentric ecliptic coordinates VxJ’’ = VxJ’’’ cos w - VyJ’’’ sin w VyJ’’ = VxJ’’’ sin w + VyJ’’’ cos w VzJ’’ = VzJ’’’ VxJ’ = VxJ’’ VyJ’ = VyJ’’ cos i VzJ’ = VyJ’’ sin i VxJ = VxJ’ cos L - VyJ’ sin L VyJ = VxJ’ sin L + VyJ’ cos L VzJ = VzJ’ By replacing J with K in all the above, you can find the sun-relative velocity in the transfer orbit at the apsidal endpoint, also. Finding the delta-vees on each end of the intended trajectory goes like this: dVx1 = Vx1 - Vx(preburn) dVy1 = Vy1 - Vy(preburn) dVz1 = Vz1 - Vz(preburn) dVx2 = Vx(rendezvous) - Vx2 dVy2 = Vy(rendezvous) - Vy2 dVz2 = Vz(rendezvous) - Vz2 In my example problem, V(preburn) is the orbital velocity of Vesta at the moment of departure, and V(rendezvous) is the orbital velocity of Earth at the moment of arrival. You can work it out explicitly, if you want to. Their magnitudes will be, roughly: dV1 = 9.2 kilometers per second dV2 = 20.4 kilometers per second THE END Jerry Abbott
Recommended Posts
Create an account or sign in to comment
You need to be a member in order to leave a comment
Create an account
Sign up for a new account in our community. It's easy!
Register a new accountSign in
Already have an account? Sign in here.
Sign In Now