Timeline for How accurate is the patched conic approximation when performing a Hohmann interplanetary transfer?
Current License: CC BY-SA 3.0
29 events
| when toggle format | what | by | license | comment | |
|---|---|---|---|---|---|
| Feb 8, 2014 at 22:22 | comment | added | Mark Adler | For that purpose, I recommend the SPICE toolkit, along with a recent ephemeris file such as DE430. You'll also need a recent leap seconds kernel. | |
| Feb 8, 2014 at 21:23 | comment | added | InquisitiveInquirer | Hi Mark. I thought it would be rude not to give you an update on what I've been doing, since you've been kind enough to help me so much. I've since moved to a Lambert solver using a basic bisection method, as the Hohmann transfer was a little too frustrating, and I felt it was far too limiting with regards to transfer opportunities (as I'm sure you're more than aware). I've also started using JPL's Horizons ephemeris for initial conditions which has made the numerical model a lot more accurate. Now I just need to get ephemeris data in a more automated manner, but that's for another topic! | |
| Jan 12, 2014 at 22:26 | comment | added | Mark Adler | The orbit in its plane is described by $a$ and $e$. The orbit plane is rotated from XY to where it actually is by $\Omega$, $i$, and $\omega$ (which are the Euler angles $\alpha$, $\beta$, and $\gamma$ for that rotation respectively). Those three orbital elements aren't "in" a plane. They describe how to transform from one plane to another. | |
| Jan 12, 2014 at 22:12 | comment | added | InquisitiveInquirer | Aha! Sorry if this is a stupid question (I haven't studied Astrodynamics), but does this mean that the eccentricity, and all other orbital elements for that matter, given by Mathematica are in the Perifocal frame, and this "rotation" puts it in the Heliocentric frame? I found a rotation matrix called QxX (took a while to put in all those trig functions!) and multiplied it by {e,0,0} and this is what my output was: {-0.01604561, -0.00466582, 1*10^-8}. Does this look correct? | |
| Jan 12, 2014 at 14:39 | comment | added | Mark Adler |
You don't have to do it that way, but you can make an eccentricity vector by rotating {e,0,0} by the Euler angles $\Omega$, $i$, and $\omega$.
|
|
| Jan 12, 2014 at 11:06 | comment | added | InquisitiveInquirer | Looking at the true anomaly wiki, it requires an eccentricity vector that is dotted with its radius vector, but Mathematica only gives an eccentricity magnitude. Is there a way to get an eccentricity vector without having to use velocity (which seems to be the only way at the moment)? | |
| Jan 11, 2014 at 23:57 | comment | added | Mark Adler |
The first five elements of AstronomicalData["Earth", "OrbitRules"] are the usual non-time-dependent orbital elements. The sixth element, the true anomaly, is a measure of the position of the object in the orbit. You will have to solve for that given the moment in time you want the true anomaly for, and the position of the object at that time, using the "Position" property as you have been doing.
|
|
| Jan 11, 2014 at 22:18 | comment | added | InquisitiveInquirer | Hi there Mark, I hope you're still around, sorry for my disappearance these last few months! I've come back to this problem and tried finding the instantaneous velocity of Earth and Mars using the orbital elements but am having some trouble. I found this link online which gives instructions on how to find velocities using orbital elements: cdeagle.com/omnum/pdf/demosvoe.pdf. The only parts of the puzzle that I can't complete, though, are finding the true anomaly v (nu) and argument of perigee w (omega) and Mathematica's AstronomicalData doesn't seem to have them either. | |
| Nov 28, 2013 at 17:12 | comment | added | Mark Adler | And in case I'm not being clear, don't do it this way. Use the orbital elements instead along with the position and calculate the velocity directly. That will give you a much more accurate result, and will not depend on the vagaries of what $\Delta T$ you pick for a numerical derivative. | |
| Nov 28, 2013 at 17:07 | history | edited | Mark Adler | CC BY-SA 3.0 |
improve function
|
| Nov 28, 2013 at 17:05 | comment | added | Mark Adler |
By the way, you can fix the part specification error by restricting the evaluation of the function to numeric arguments with: f[dt_?NumericQ]:=... I should have done that initially.
|
|
| Nov 28, 2013 at 16:57 | comment | added | Mark Adler |
You need /(2 dt), not /2 dt. The latter is multiplying by dt not dividing by it. Once you fix that, it's still not apples and apples, since now your "dt" is actually half of the real dt. To have a consistent meaning for dt, you need to do 1+dt/2 and 1-dt/2, and the divide the difference by dt instead of (2 dt).
|
|
| Nov 28, 2013 at 15:00 | comment | added | InquisitiveInquirer |
When I use your error formula in the form of a central difference I get a very difference story, although I do get a error saying the following: Part::partd: "Part specification 0[[1]] is longer than depth of object." This is the error checking formula I used for the central difference:f[dt_] := AstronomicalData[ "Earth", {"Position", {2001, 3, 30, 0, 0, 1 + dt}}] - AstronomicalData["Earth", {"Position", {2001, 3, 30, 0, 0, 1 - dt}}] LogLinearPlot[f[dt][[1]]/2 dt, {dt, 0.001, 100}]
|
|
| Nov 25, 2013 at 16:12 | comment | added | Mark Adler |
You need to check each time, since it depends on the data you're using. There is an ND (numerical derivative) function in the NumericalCalculus package you can try. I've never tried it, but it may do the checking for you. Even then I would still check myself to make sure the answers make sense. Doing numerical derivatives is fraught with peril.
|
|
| Nov 25, 2013 at 13:41 | comment | added | InquisitiveInquirer | Oh wow, I always thought that one needed to get ΔT as small as possible. If anomalies like you've shown do occur, how would I know what value of ΔT would be appropriate for use (not just in this case but for any case where I have to choose a ΔT when using a finite difference)? Is there a standardised method, or do I just have to make an ad hoc check using the formula you produced above? | |
| Nov 21, 2013 at 18:27 | comment | added | Mark Adler | I added the plot to the answer. | |
| Nov 21, 2013 at 18:24 | history | edited | Mark Adler | CC BY-SA 3.0 |
add plot of numerical derivative noise
|
| Nov 21, 2013 at 17:57 | comment | added | Mark Adler |
Perfect example. 0.002 seconds is way too small of a $\Delta T$ for the data you are using. Yes, you want a small $\Delta T$, but it is possible to pick too small of a $\Delta T$ that samples noise in the lower significant digits of your data. Run this in Mathematica to see what I'm talking about: f[dt_] := AstronomicalData["Earth", {"Position", {2001, 3, 30, 0, 0, 1 + dt}}] - AstronomicalData["Earth", {"Position", {2001, 3, 30, 0, 0, 1}}] and LogLinearPlot[f[dt][[1]]/dt, {dt, 0.001, 100}]. From that plot, it looks like around 10 seconds is a good $\Delta T$.
|
|
| Nov 21, 2013 at 9:58 | comment | added | InquisitiveInquirer | Large ΔT? I thought we would want ΔT to be as small as possible in order to bring it as close to an infinitesimal as we can, or am I misunderstanding something? Currently I have ΔT = 0.001 and, for example, the central difference I used to calculate Earth's initial velocity is $$v[1] = ((AstronomicalData[ "Earth", {"Position", {2001, 3, 30, 0, 0, 1.001}}] - AstronomicalData[ "Earth", {"Position", {2001, 3, 30, 0, 0, 0.999}}]))/(2 (0.001));$$ | |
| Nov 20, 2013 at 14:03 | comment | added | Mark Adler |
That will work. Make sure that you pick a large enough $\Delta T$ to get a good number of significant digits in the differences. Alternatively you can use the orbital elements from AstronomicalData[], e.g. "SemimajorAxis", "Eccentricity", etc. and propagate the two-body solutions for the Earth, Moon, and Mars in closed form. You really only need to integrate the trajectory of the spacecraft.
|
|
| Nov 20, 2013 at 10:04 | comment | added | InquisitiveInquirer | Well said, going back to the Hohmann transfer then, how would I take into account Earth's velocity at departure in the above equation? Currently I've been using Mathematica's AstronomicalData function to get the initial positions of Earth and Mars in space and then I've found their initial velocities using a central difference on their positions over a relatively small time period. These values are then put in as initial conditions for an N-body simulation. | |
| Nov 19, 2013 at 15:16 | comment | added | Mark Adler | Yes, to find the optimal connection between two real orbits with eccentricities, inclinations, and arguments of periapsis, you would use the Lambert's problem solution. However I recommend that you get the Hohmann transfer between two coplanar elliptical orbits working first, and then graduate to the Lambert's problem solution. | |
| Nov 19, 2013 at 11:40 | comment | added | InquisitiveInquirer | Hmmm, I've been doing a little more reading and it seems like a pure Hohmann transfer is not usually used in interplanetary travel as it's not always practical and instead a solution to Lambert's problem is found. Should I maybe consider moving over to something like that? | |
| Nov 18, 2013 at 0:13 | comment | added | Mark Adler | A second by the way: that equation assumes that the starting orbit is circular. Earth's orbit is not circular (in fact no orbit is really), and so you need to account for the Earth's actual velocity at departure. | |
| Nov 17, 2013 at 17:31 | comment | added | Mark Adler | By the way, you can use $\TeX$ equations in comments, e.g. $\sqrt{\mu_S\over r_1}\left(\sqrt{2r_2\over r_1+r_2}-1\right)$. I think that's the equation you meant, but it looks like you're missing a parenthesis. | |
| Nov 17, 2013 at 17:20 | comment | added | Mark Adler | I can't quite make out your equation, but yes, use Earth's heliocentric radius at departure, and Mars' heliocentric radius at arrival. | |
| Nov 17, 2013 at 12:07 | comment | added | InquisitiveInquirer | When you say I must take into account the eccentricity of Mars and Earth, I'm guessing you mean I must change the equation for v_infinity which is given by v_infinity = sqrt(mu_sun/R_earth(sqrt(2R_mars/(R_mars+R_earth))-1)? Would I set R_earth to Earth's heliocentric radius at spacecraft departure, and R_mars to Mars' heliocentric radius at spacecraft arrival? Currently they are both set to their heliocentric radii at departure. P.S. I apologise for the horrible looking equation, it's tough to make them look nice in a little box like this. | |
| Nov 16, 2013 at 5:15 | history | edited | Mark Adler | CC BY-SA 3.0 |
response to question edit
|
| Nov 15, 2013 at 15:53 | history | answered | Mark Adler | CC BY-SA 3.0 |