Orbit Visualization and Ephemeris Generation:



Orbit Visualization and Ephemeris Generation:

One of the first steps in orbit determination is ephemeris generation. An ephemeris is a table of values that gives the position of an object in the sky as a function of time. With a known orbit we want to tabulate the future positions and velocities of the orbiting object.

The big picture of this assignment is that given the classical orbital elements (a,e,v,i,Ω,ω) we wish to display the orbit of this asteroid using Newtonian gravitation and a central force model and calculate the position and velocity of that object at that time.

Recall that previously you placed the earth, sun and your asteroid in a solar system. You drew the orbit based on an entered eccentricity, semi-major axis and inclination. You drew arrows that connected the sun to the earth (Something we obtain from JPL from the Horizon system), the earth to the asteroid (What we measure from the CCD) and the sun asteroid vector (Determined as the sum of the earth sun and earth asteroid vector.) The output of your program should have looked something this:

Now we would like to implement the other orbital elements into our model. If we can model these elements completely, we can run our model forward and backward in time and find the position of the asteroid at any time. This is the process of ephemeris generation.

1. Make any necessary modifications to your program from Monday so that you can display the orbit of your asteroid given any (e,a and i)

2. To implement the inclination rotation, you are rotating about the x axis. You probably implemented the inclination by assuming that the initial momentum was all in one direction, say the z direction and then multiplying the y and z components of the momentum by -sine(i) and cosine(i) respectively. A more elegant way to do this is to apply a rotation matrix.

3. Vpython provides a simple way to apply these rotation matrices using the “rotate” command.

newvector = rotate(oldvector, angle=pi/4., axis=axis, origin=pos)

The rotate function applies a transformation to the specified object oldvector.). The transformation is a rotation of angle radians, counterclockwise around the line defined by origin and origin+axis. By default, rotations are around the object's own position and axis. Modify you program so that it uses the “rotate” command to implement the inclination.

HINT: You might want to display a set of X, Y and Z axis centered on the Sun to help you visualize the rotations.

Make sure you understand how the rotate command works before proceeding!

4. Now you need to implement the longitude of the ascending node(Ω). Based on the diagram at the very top, you should see that this is a second rotation of the orbit, but this time about the y axis. Use the vector rotate command to implement the longitude of the ascending node. (Hint: you will need to rotate both the momentum vector and the position vector)

5. Draw the line of nodes. This is a straight line from the ascending node to the descending node. Also draw a line to indicate the origin of the system. Call Prof Mason or a TA over when you are done so that they can see you work.

[pic]

6. Now implement the argument of periapsis (ω). This is yet another rotation. However this one is a bit trickier. By inspection of the diagram at the top, it is clear that we need to rotate the ellipse about an axis that is perpendicular to the plane defined by the orbital plane of the ellipse. Use another rotation of the position vector and momentum vector about an axis perpendicular to the plane of the ellipse to implement the argument of periapsis. Hint: you might want to think about the definition of the cross product. If you get stuck, talk to your neighbor or ask for help. Keep checking your progress against the Wikipedia picture at the vey top and the picture above.

If you use the following values you should produce a picture that pretty closely matches the Wikipedia picture and the picture below:

e = 0.4,a = 3e11, i = .80,Ω = 4.72,ω = 1.57

Call Prof Mason or a TA over when you have finished your work.

[pic]

7. The last element is the “starting point” of the asteroid along its orbit. We can’t simply do another rotation, because the asteroid needs to propagate forward along the ellipse, not on a circle! This will require some thinking!!!!

The simplest thing to start with is the true anomaly. The true anomaly is the angle between the direction of perihelion and the current position of the asteroid as seen from the sun.

If we know the true anomaly, we can write the rectangular coordinates of the ellipse as

[pic]and taking the derivative [pic]

Where:

[pic]

Note that in the case that ν = 0, the expression reduces to[pic] or

[pic]

Which simplifies to:

[pic]

This is equivalent to the perigee velocity we derived in class.

Implement these conditions as the initial conditions for the position and velocity of your asteroid. Your asteroid should now start off rotated on the ellipse.

Mean Anomaly:

Unfortunately, True Anomaly is not typically reported in the description of orbits! What is typically reported is Published Mean Anomaly. I didn’t make this stuff up so please bear with me for a bit! The first thing we have to do is get the Mean anomaly(M) from the published Mean anomaly(Mp). In order to do this we also need to know the ephemeris time. This is simply a time when astronomers agree to start their results. The MPC recomputes this every year and most recently, ephemeris time (TE) is June 18,2009. Since time is linear, we just calculate the Mean Anomaly (M) by how much time has elapsed since the ephemeris time.

[pic]

Where t is the current Julian Date and Mp is the published Mean Anomaly

NOTE: Mp is usually published in DEGREES. However the units in this equation assume that you are working in radians.

Once you have the Mean Anomaly, you can calculate the eccentric anomaly by inverting the Kepler Equation using Newton’s Method: M = E – e sin E. Recall that we discussed this in class on Monday.

Exercise: Write a function to take a value for M and a tolerance and return the value for E. Here is some sample data to help:

e = 0.5 M = 45o E = 72.29o

Note: You can also check you function by plugging E back into M = E – e sin E

Once we have the Eccentric Anomaly, we are in the homestretch. Dr Ran derived the rectangular coordinates of the initial position and velocity of the asteroid as a function of a, e and E as:

[pic] and [pic]

Where this position is in ecliptic coordinates with the x axis coinciding with the perihelion.

Now change your program to accept Published Mean Anomaly and Ephemeris Time instead of True anomaly.

[pic]

Be sure to get Prof Mason or a TA to check you off at this point.

Congratulations! The hard part is done. We can generate the position and velocity vectors in ecliptic coordinates of an asteroid at a given (a,e,Mp,i,Ω,ω) and Te.

These are just the position and velocity vectors your program generates after applying the Mean Anomaly and the three rotations.

8. First we need to convert the coordinates into ecliptic(The earth velocity vector given below has a non zero z component!) and then equatorial coordinates.

Convert r to equatorial by rotating the coordinates by the tilt of earth relative to the ecliptic ε:

[pic]

This is a rotation around the x axis. Is that still appropriate in our coordinate system? Where ε = 23.439281o

9. Now determine the new earth asteroid vector (also known as the range vector or ρ) by adding the earth-sun vector to the new r vector. To do this accurately, we need to put the earth in the right place relative to the sun! Modify your program so that it places the earth in the right place in the sky relative to date and time by incorporating your program from the homework.

10. Now solve for your RA and DEC using [pic] Be sure to worry about angle ambiguity!

You ephemeris project is to implement all of these ideas in a program that takes a set of orbital elements plus a JD and generates the RA and DEC of the object at that time.

-----------------------

Wikipedia 2008

[pic]

Origin Line

Line of Nodes

Wikipedia 2008

................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download