Orbit Determination Using Matlab

[stupid6625]

I am brand new to the physics forum community, so let me start off by saying hello to everyone. I am interested I am coming up with a MATLAB code for orbit determination. I saw a thread recently discussing this type of topic but i think it was orbit determination using RA and Dec measurements only. I was thinking more along the lines of using a given .TLE file (Two Line Element) and going from there. This said, i am not too familiar with TLE's or how they are used. Any help would be appreciated as well as any example MATLAB codes as i am not greatly proficient with MATLAB or any other language for that matter. Thanks guys!

 

[Jenab]

I am brand new to the physics forum community, so let me start off by saying hello to everyone. I am interested I am coming up with a MATLAB code for orbit determination. I saw a thread recently discussing this type of topic but i think it was orbit determination using RA and Dec measurements only. I was thinking more along the lines of using a given .TLE file (Two Line Element) and going from there. This said, i am not too familiar with TLE's or how they are used. Any help would be appreciated as well as any example MATLAB codes as i am not greatly proficient with MATLAB or any other language for that matter. Thanks guys!

I wrote a post on how to determine the orbital elements of an asteroid from three angle-only (i.e., RA and DEC) positions. All that rigmarole by which the distances are successively approximated is avoided when you know the distances in advance. Given the distances as well as the angular positions, orbit determination is easy.

I'm not sure what "two line element" means, since I'm not a professional astronomer; I'm just a goatherd who had some astronomy classes in school.

If you're given the elements at the start and are calculating positions in space, then you're not doing orbit determination. You're using an already determined orbit to construct an ephemeris.

Jerry Abbott

 

[stupid6625]

Ok, so what are the measurements that i would need in order to predict an orbit? For example, i read through the posts about determining an orbit using just RA and DEC. The end result for that would be a complete analysis of the orbit correct? You would then know the inclination, the eccentricity, the RAAN, etc.?? If so, then maybe that is the route i should take. I am thinking of this more along the lines of a grad school type project and am really just trying to lay the groundwork for what would be involved. I would still love to see any MATLAB codes that anyone has just so i can get an idea of the logic behind it all. Thanks again!

 

[Jenab]

Ok, so what are the measurements that i would need in order to predict an orbit? For example, i read through the posts about determining an orbit using just RA and DEC. The end result for that would be a complete analysis of the orbit correct? You would then know the inclination, the eccentricity, the RAAN, etc.?? If so, then maybe that is the route i should take. I am thinking of this more along the lines of a grad school type project and am really just trying to lay the groundwork for what would be involved. I would still love to see any MATLAB codes that anyone has just so i can get an idea of the logic behind it all. Thanks again!

Well, if my thread at
https://www.physicsforums.com/showthread.php?t=36657
shows how to go from three angle-only (RA & DEC) sky positions at three known times, provided that the geocentric positions of the sun at those times are also known.

Input.
{RA(i), DEC(i), t(i), Xs(i), Ys(i), Zs(i)}
Where RA(i) is the local right ascension of the asteroid at t(i).
Where DEC(i) is the local declination of the asteroid at t(i).
Where t(i) is the time of observation.
Where [Xs(i), Ys(i), Zs(i)] is a vector from observer to the sun at t(i).

Intermediate Quantities.
tc(i), r(i), p(i)
Where tc(i) are the times of observation corrected for speed of light lag.
Where r(i) are the distances from sun to asteroid at the corresponding tc(i).
Where p(i) are the distances from Earth to asteroid at the cooresponding tc(i).

Output.
{ a, e, i, L, w, T}
The elements of the asteroid's orbit,
Where a is the semimajor axis.
Where e is the eccentricity.
Where i is the inclination to the ecliptic.
Where L is the ecliptic longitude of the ascending node.
Where w is the argument of the perihelion.
Where T is the time of perihelion passage.

However, if you want to assume that you already know the distances r(i) and p(i), then you can skip all the first stuff and proceed directly to the procedure given in the 3rd and 7th posts in the thread.

Jerry Abbott

 

[enigma]

"Fundamentals of Astrodynamics and Applications" by Vallado has pseudocode for practically everything you'd want. The downside is it's poorly edited (IMHO) and you'll have to hunt and peck through the entire book. Take a look through Jenab's hillbilly tutorials. He wrote up several examples for doing various problems in celestial mechanics. You should be able to withdraw some pseudocode from that.

If you have a specific problem you're trying to work through in Matlab, I can go through my old orbital dynamics code and see if I can find what you're looking for. The files much too long (and open up potential for lurker cheating) to simply post them all here.

Welcome to the forums, BTW!

 

[stupid6625]

Ok, that actually helps out quite a bit. The Thread about finding the orbit based on RA and Dec is for an object in a heliocentric orbit correct? What if i was just interested in orbits just around Earth? This would change things quite a bit right?
Also, in the example calculations, what are the "tau" "nc" "nu" "K0" "L0" representative of? I think i know what tau and nu are representative of, but i am unsure about a lot of the variables in that code. And in that thread it looks like there are several additions or corrections to the method...What would the final method look like?

 

[BobG]

I know absolutely nothing about MATLAB, but I do know about orbital elements.

This gives the format of two-line elements – what each element means.

http://celestrak.com/NORAD/documentation/tle-fmt.shtml

The fourth element is the epoch time. In other words, the satellite was in the position given by the elements at the time specified by the epoch time.

The first derivative, second derivative, and B*star are for orbit propagation. They define atmospheric drag characteristics of the satellite. To actually predict where the satellite will be in the future, you need these, but that’s something you add on to your program (or spreadsheet, in my case) after you’ve got the basics down.

The second line contains info on the satellite’s position and velocity. Three of the elements fix the satellite’s position and velocity within the orbital plane (mean motion, eccentricity, and mean anomaly). Three rotate the satellite’s position and velocity from a perifocal coordinate system (position/velocity within the orbital plane) to the geocentric equatorial coordinate system (the satellite’s position and velocity relative to the center of the Earth (right ascension, inclination, argument of perigee).

You have two of the second line elements that have to be changed right off the bat. If you know Kepler’s third law, you can calculate the semi-major axis from the mean motion. Mean motion in your two-lines are expressed in revolutions/day instead of radians/second, so you have to convert the units to radians/second. The true anomaly is needed instead of the mean anomaly. This is a two step method, normally using the Newton-Raphson method to go from mean anomaly to eccentric anomaly (there is no analytic solution – you have to use a trial and error method to solve the equation), then using a more standard analytic equation to go from eccentric anomaly to true anomaly.

If you don’t care about understanding how you get from two-lines to geocentric coordinates, you can do it in about 6 equations. If you want to understand what’s happening, you need a book’s worth of equations so you can step things through one step at a time. As enigma kind of hinted at, knowledge of orbital mechanics and an ability to express oneself coherently seem to be mutually exclusive talents.

Angles are expressed in greek letters, linear numbers in western letters.

a= semi-major axis
p= semi-latus rectum, or p=a(1-e^2)
r= satellite radius (the velocity is wholly dependent on the satellite's position, so you never see the velocity in the equations)
$$\Theta$$ or Theta (with a capital letter) is the sum of argument of perigee and true anomaly (i.e. - the satellite's location along the orbital ellipse relative to the equator)
$$\nu$$ or nu is the true anomaly
$$\iota$$ or i is the inclination
$$\Omega$$ or Omega (with a capital letter) is right ascension of ascending node
$$\omega$$ or omega (with a small letter) is argument of perigee
$$\tau$$ or tau is usually orbital period

All the definitions depend on the book you happen to be using, but these are the most common.



$$x=r(cos\Theta cos\Omega - sin\Theta sin\Omega sin\iota)$$
$$y=r(cos\Theta sin\Omega + sin\Theta cos\Omega cos\iota)$$
$$z=r(sin\Theta sin\iota)$$

$$\Theta = \omega + \nu$$

$$\dot x = - \sqrt{\frac{\mu}{p}} [cos\Omega \left(sin \Theta + e sin \omega) + sin \Omega cos \iota (cos \Theta + e cos \omega)]$$
$$\dot y = - \sqrt{\frac{\mu}{p}} [sin\Omega \left(sin \Theta + e sin \omega) - cos \Omega cos \iota (cos \Theta + e cos \omega)]$$
$$\dot z = \sqrt{\frac{\mu}{p}}[sin \iota ( cos \Theta + e cos \omega)]$$

$$\mu$$ is the geocentric gravitational constant (3.986 x 10^5 km^3/sec^2)