Orbits - a demonstration of planetary orbits for the hp 39g+

Introduction

Astronomers know that the orbits of planets around the Sun (or of satellites around the Earth) are elliptical, but it is not always easy to visualise their shapes or to appreciate the effect of the 'eccentricity' of orbits. This aplet draws the shapes of typical orbits and also demonstrates Kepler's 2nd law of planetary motion.


Orbits on the hp 39g+ Graphing Calculator

Orbits uses the general equation for a 'conic section' to draw orbits with various eccentricities. As well as seeing the shape of orbits, the Numeric view can be used to provide a more quantitative feel for the difference that eccentricity makes.


Loading Orbits

Orbits is an aplet. Copy the files HP39DIR.000, HP39DIR.CUR, ORBIT4T0.000, ORBITECC.000, ORBITKP2.000, ORBITKXY.000, ORBITR00.000, ORBITS00.000, ORBITS01.000 and ORBITSV0.000 to an empty directory on a PC, start the HP39G Connectivity program and point it to the directory holding the files. Transfer Orbits from the computer into the APLET catalog part of the hp 39g+'s memory using the RECV option on the calculator and selecting Orbits as the file to download. The other files are linked as part of the aplet and should automatically transfer into the PROGRAM catalog.


Deleting Orbits

Delete Orbits from the APLET catalog. The constituent programs then need to be individually deleted from the PROGRAM catalog. The files are: .Orbit.4T, .Orbit.ECC, .Orbit.KP2, .Orbit.S, .Orbit.R, .Orbit.SV and .Orbit.KXY. (I have read that the original design for the hp39 called for programs linked to an aplet to be automatically deleted with it, but for some reason this was not done.)


Theory of Orbital Motion

Textbooks on mechanics and kinematics (e.g. chapter 9 of 'An Introduction to Mechanics' by Kleppner and Kolenkow, McGraw-Hill 1981) show that the gravitational force between any two objects is given by:
F = GMm/R²
where F is the force, G is the constant of gravitation, M and m are the masses of the two objects, and R is the distance between them. It is usually assumed that the mass of M (e.g. the Sun) is much larger than the mass of m (e.g. a planet) so that M can be thought of as effectively stationary with m moving around it.

It can further be shown that the position of one object relative to the other can be represented by the very simple equation:
R = R0/(1 - e*cosq)
where R is the distance between them, R0 determines the overall size of the orbit (the value of A below), e is the eccentricity which is a measure of how much the orbit differs from a circle, and q is the angle measured around the orbit from the direction of greatest separation, as in Figure 1:
Definition of orbital terms
Figure 1 - Definition of orbital terms

The shape of the orbit is fixed purely by the value of the eccentricity, e.
A, called the semi-major axis, only enlarges or reduces the size.
A is related to R0 and the eccentricity by A = R0/(1 - e²).

Notice that the larger mass is not at the centre of the ellipse but offset from it, at the focus. The distance by which the focus is offset from the centre is equal to A*e so that the larger the eccentricity the further away the focus is from the centre.


Conic Sections

Picture a circular-base cone, like an inverted funnel. Now imagine taking slices through it at various angles and looking at the shape of the exposed surface. It turns out that the shape of these so-called 'conic sections' can also be represented by the equation R = R0/(1 - e*cosq), where R0 is now a measure of how far down from the tip the slice is made. There are four distinctly different situations, as drawn in Figure 2:
Conic sections
Figure 2 - the four conic sections. The direction of the slice is shown by the yellow rectangle and the exposed surface is outlined in red.


The four possibilities are:

  1. If the slice is parallel to the base of the cone then the surface is a circle, with a constant radius, say equal to A. This corresponds to eccentricity e=0.
    Circular orbitA circular orbit, with constant radius R=A

  2. If the cut is made at a slight angle to the base, but less than the angle of the sides of the cone, the resulting shape is an ellipse whose diameter along the long axis is usually represented by 2A. The eccentricity in this case is more than zero but less than one, being larger the steeper the angle of the cut, i.e. 0<e<1.
    The radius of the ellipse along the shorter axis is normally written as B and can be calculated from B=AÖ(1-e²).
    Clearly the larger e is the smaller B is relative to A and the more 'squashed' the ellipse becomes.
    Elliptical orbitAn elliptical orbit, eccentricity approx 0.5.

  3. When the cut is made exactly parallel to the side of the cone, one side of the cut never leaves the cone (which in theory extends to infinity) so the edges of the surface never curve back onto themselves. The shape is a parabola and the corresponding eccentricity is e=1. When considered as an orbit, because it is not 'closed' the object would make only a single pass by the Sun and then recede to infinity, never to return.
    Parabolic orbitA parabolic orbit.

  4. Finally, if the cut is made at a steeper angle than the side of the cone the resulting shape is a hyperbola. This looks quite similar to a parabola near the curved end but at the open end the two arms diverge more rapidly than for a parabola. The eccentricity for a hyperbola is more than 1, that is e>1.
    Hyberbolic orbitA hyperbolic orbit.

When considering the movement of real objects under gravity, it is often useful to think in terms of energy and escape velocity. Escape velocity is defined (Oxford Dictionary of Astronomy) as the minimum velocity required to escape from a gravitational field, and can be calculated as
V = Ö(2GM/R) where as before G is the gravitational constant, M is the mass to be escaped from and R is the distance of the escaping object from the centre of mass.

What are generally described as 'planets' orbit the Sun continuously since they have too little energy to escape from the Sun's gravity and always travel at less than escape velocity. Their orbits are 'closed' and are ellipses (where a circle can be thought of as just a special case of an ellipse).

'Comets' on the other hand are travelling faster at a given distance from the Sun than planets are and thus have more energy. If their speed is just less than escape velocity they will have a highly elliptical orbit with eccentricity almost 1. (Halley's comet has an eccentricity of 0.967.)
If their speed is exactly equal to the Sun's escape velocity then their orbit is parabolic and after passing by the Sun they will continue to move away from it but at a decreasing rate so that eventually they will in theory become stationary relative to the Sun.
For a speed greater than escape velocity the orbit is hyperbolic and the comet will continue to move away from the Sun without ever coming to rest.
The dividing line between highly elliptical/parabolic/hyperbolic is quite fine and it can be difficult to distinguish between a comet with a very long period elliptical orbit and one which will never return.


More on Elliptical Orbits

Since all planets, asteroids, moons and artificial satellites follow elliptical orbits this type will be described in more detail. Refer to Figure 3.
Parts of an ellipse
Figure 3 - Key parts of an ellipse


The distance from the centre of the ellipse to the circumference along the longer axis is the semi-major axis, = A.
Half the width across the shorter axis, called the semi-minor axis, = B.
The eccentricity e can be calculated as e² = 1 - B²/A²
This can be rearranged as B² = A²(1 - e²) from which it can be seen that if the eccentricity is zero, B=A and the ellipse becomes a circle.

For a planetary orbit the Sun is at one focus of the ellipse rather than at the centre; there is nothing at the other focus.
The distance between the centre and either focus, labelled as C in Figure 3, is given by C = Ae.
The distance marked as D, between the focus containing the Sun and the closest part of the orbit, is called the perihelion ('near sun') distance and is calculated as D = A(1 - e).
Similarly the maximum distance from the Sun to the planet, when it is at the opposite end of the major axis, is called the aphelion ('away from sun') and is equal to A(1 + e).

As stated earlier, the general equation for an ellipse using polar coordinates is
R = R0/(1 - e*cosq) where the centre of the coordinate system is at one focus of the ellipse. If the length of the semi-major axis is known then the general polar equation for an ellipse with semi-major axis = A becomes
R = A(1 - e²)/(1 - e*cosq)     since A=R0/(1-e²).

Alternatively, if cartesian coordinates are used, the equation for an ellipse can be specified as
X = A*cosq, Y = B*sinq, where as before A and B are the semi-major and semi-minor axes. Note that in this case however the zero of the coordinate system is at the centre of the ellipse.

The 'height to width ratio' of the ellipse, that is the ratio of semi-minor to semi-major axis, is calculated as simply
B²/A² = (1 - e²), thus B/A = Ö(1 - e²).
The larger the value of the eccentricity e (i.e the closer to 1) the smaller is the square root of (1-e²) and so the flatter is the ellipse. The shape of various ellipses, all with the same length of semi-major axis A but different eccentricities, is illustrated in Figure 4.
Effect of eccentricity on shape of ellipses
Figure 4 - Ellipses with eccentricities of 0, 0.2, 0.4, 0.6, 0.8 and 0.95

From Figure 4 it can be seen that the effect of eccentricity on the shape of the orbit is not very marked for eccentricities below about 0.5, but the shift of the perihelion point towards the focus is clear. Some typical figures for perihelion distance and ratio of B/A versus eccentricity are:
EccentricityRatio B/APerihelion distance if
A=1
0.01.0001.00
0.10.9940.90
0.20.9800.80
0.50.8660.50
0.90.4360.10
0.950.3120.05
0.990.1410.01


Application to Planets in the Solar System

Most of the planets in our solar system have small eccentricities:
Planet             Eccentricity
Mercury0.206
Venus0.007
Earth0.017
Mars0.093
Jupiter0.048
Saturn0.055
Uranus0.051
Neptune0.007
Pluto0.252
Halley's comet0.967

The early Ptolemaic theory which proposed that the planets moved in perfect circles was reasonably successful because all the planets except Mercury and Pluto do have orbits which are nearly circular. Mercury is always close to the Sun and not often visible, and Pluto was undiscovered at the time.
It is probably fortunate for us that the Earth's orbit has a low eccentricity since a high eccentricity, like Mercury's, would greatly alter the intensity of sunlight at different times of the year and either strengthen or weaken seasonal temperature differences. Even now the Earth receives a few percent more solar heat at perihelion (in January) than at aphelion (in July). It is thought that long-term variations in the eccentricity of Earth's orbit may be a factor in the cycle of ice ages and warm spells over the last few million years.


Kepler's Laws

The German mathematician and astronomer Johann Kepler used the accurately measured positions of the planets determined by Tycho Brahe to calculate their orbits. He found that circular orbits could not be fitted to the observations and deduced that instead the planets move in ellipses. Between 1609 and 1621 he derived his three laws of planetary motion:

  1. The planets move in an ellipse with the Sun at one focus.
  2. The imaginary line joining a planet to the Sun (the 'radius vector') sweeps out equal areas in equal times. (Figure 5)
  3. The square of the orbital period of a planet (in years) is proportional to the cube of the semi-major axis of its orbit,
    T² µ A³   or A³/T² = constant.

Kepler's second law
Figure 5 - Kepler's 2nd Law. The blue arrows represent the distance around its orbit that a planet moves in the same particular time interval T. It moves more rapidly and therefore along a greater length of orbit when it is close to the focus than when it is further away. The green and red shaded areas are the parts of the orbit swept out by a line from the planet to the focus in time T, around perihelion and around aphelion. Although the shapes of the two ellipse segments are very different they have the same area.

The third law was particularly important since it would allow the orbital distances of any newly-discovered planets or asteroids to be calculated from their periods.
In the 1660s Isaac Newton developed the theoretical basis for the Kepler's first and second laws, and showed that the third law should be amended to
A³/T² = G(M+m)/(4p²) where M is the mass of the Sun and m is the mass of the planet.
In other words the orbital period of a planet depends on the total mass of the Sun plus the planet, which is not quite a constant for all planets. In practice the masses of the planets are much smaller than that of the Sun and Kepler's law is often a good enough approximation. An especially useful result of Newton's version is that it is possible to calculate the mass of a planet, such as Jupiter, if the orbital period and orbital size of one of its moons is measured. The mass of the planet is very much larger than that of the moon so the combined mass is almost entirely that of the planet.


Running Orbits

Start Orbits from the calculator's aplet catalog. This brings up the VIEWS menu below.

Orbits VIEWS menu
The main menu

Types of Orbit

The first menu option plots an illustration of the four different types of orbit: circular, elliptical, parabolic and hyperbolic. The centre of the axes represents the focus, where the Sun is, and it can be seen that the perihelion distance decreases for the four types, as below.

Four orbit types
The plot of four orbit types

Pressing SYMB displays the equations used. The eccentricities are 0.0, 0.5, 1.0 and 1.5 for R1 to R4.

Pressing NUM for the numeric view allows the way the radius of the orbit varies with angle, measured from the +X axis, to be investigated for the four types of orbit.

  • For the circular orbit R1 is always 1.
  • For the elliptical orbit R2 varies from 2 at q=0 to 0.667 at q=180 degrees.
  • For a parabolic orbit R3 is undefined at q=0 degrees because the parabola extends to infinity, but the radius then reduces rapidly down to 0.5 at q=180 degrees.
  • In the case of the hyperbolic orbit, R4 is undefined for all angles less than 60 or more than 300 degrees because the open shape of the hyperbola means the object moving in a hyperbolic orbit never passes through these regions. The actual angle limits depend on the eccentricity and can be calculated as q=arctanÖ(e²-1), or 360 minus this value. At 180 degrees the radius reduces to 0.4, demonstrating that a hyperbolic orbit means passing close by the Sun.
Press VIEWS again to return to the menu.

Orbital Eccentricity

The second menu option shows how the value of the eccentricity affects the shape of elliptical orbits. Four orbits are plotted, all with the same length of semi-major axis (1 unit) but with eccentricities of 0.1, 0.3, 0.7 and 0.9. See below.

Effect of eccentricity
The effect of eccentricity on orbit shape

All these orbits would have the same orbital period (the same 'year' for the planet) but obviously differ considerably in perihelion and aphelion distances. The e=0.1 orbit looks almost circular whereas the e=0.9 one is very flattened.

Pressing SYMB shows the equations but if editing these note that the eccentricity (e.g. 0.1 for R1) appears twice, in order to keep the major axis length the same.

Press NUM to get the numeric view and see how the aphelion distance (at q=0 degrees) varies from 1.1 to 1.9 and similarly the perihelion decreases from 0.9 to 0.1. The VIEWS key returns to the menu.

Kepler's 2nd Law

The third VIEWS menu option gives an animated portrayal of Kepler's 2nd law. An elliptical orbit with eccentricity=0.5 is drawn then pairs of lines are drawn from the position of the Sun to points on the orbit representing the location of the planet at the start and end of a given time period. The time period between adjacent lines is equivalent to 1/12 of a full orbit and the second pair of lines represents the same time interval half an orbit later. See below.

Animation of Kepler's second law
Illustration of equal areas in equal times.

The pairs of lines, which are radius vectors, are repeatedly redrawn in steps equal to 1/36 of an orbital period so that the way they move around the orbit can be seen. The lines are widest apart close to perihelion, when the planet is moving rapidly, and closest together at aphelion when the planet moves at its slowest speed. However the area enclosed by the two lines and the arc of the ellipse remains constant all the way round the orbit, equal to 1/12 of the total area of the ellipse. Thus the radius vectors sweep out 'equal areas in equal times'.
After a full orbital period the demonstration returns to the main menu.

The final two menu options simply duplicate the effect of the RESET and START aplet buttons.


Hardware Requirements

Orbits runs on a Hewlett-Packard 39g+ calculator. It should also work on the 39g, 39gs and 40gs, and may work on the 48g, 49g and 50g but has not been tried. If the four statements DISPXY -1;1.8;2;"Kepler's 2nd Law":DISPXY 2.5;1.2;1;"Equal areas":DISPXY 3.2;0.8;1;"in":DISPXY 2.5;0.4;1;"equal times": are deleted from program .Orbit.KP2 then it ought to run on the 38g as well, without any loss of function but more slowly.

The aplet and associated programs take up 7.5 kilobytes of RAM.


Variables Used

Orbits uses most of the HOME single-value variables including q, and will therefore overwrite any existing information stored in them.


To-Do

A demonstration of Kepler's 3rd law would be nice, and possibly a discussion of transfer orbits. An explanatory note could be attached to the aplet but it would be difficult to include a useful amount of information without taking up too much RAM. It is probably better to leave the explanation in this external file.


Known Bugs

If the program is used 'as-is' there should not be any faults, but if the equations are modified the axes scales will need to be adjusted too. Auto scaling is not really suitable because it does not keep the scales the same on both axes and circles will no longer be circles.
The calculation of the true anomaly when demonstrating Kepler's 2nd law is a little slow, but the alternative method using the equation of the centre is not accurate enough for eccentricities greater than about 0.3, and smaller eccentricities do not show the variation in length of radius vector clearly enough.


Program Listings

Note: the character Þ represents the 'STO' arrow.
True line breaks in the program are shown by ¿, other line breaks are just to fit the listing into the table.


Program .Orbit.SVExplanation
Defines the VIEWS menu. Only needs to be run once, when the aplet is first written, or if the views menu is changed.
SETVIEWS
"4 Orbit types";".Orbit.4T";1;
"Eccentricity";".Orbit.ECC";1;
"Kepler's 2nd";".Orbit.KP2";7;
"Start";".Orbit.S";7;
"Reset";".Orbit.R";10;
" ";".Orbit.SV";10;
" ";".Orbit.KXY";10:¿		Define the menu items and specify which program is to be run when they are selected, and what view to return to afterwards. Also lists programs .Orbit.SV and .Orbit.KXY with a menu name of a single space, so that these programs are linked to the aplet but no menu item appears.

Program .Orbit.RExplanation
Called when RESET is selected and to initialise the display and polar equation definitions at other times.
ERASE:¿Clear the display
1ÞAngle:0ÞSimult:0ÞLabels:
1ÞAxes:1ÞConnect:
0Þqmin:360Þqmax:1Þqstep:
0ÞNumStart:10ÞNumStep:¿		Set angle measure to degrees, sequential graph drawing, labels off, axes on, polar angle range 0-360 with a step of 1, and numeric view to start at 0 with steps of 10.
0ÞR1(q):0ÞR2(q):0ÞR3(q):
0ÞR4(q):0ÞR5(q):0ÞR6(q):¿
UNCHECK 1:UNCHECK 2:UNCHECK 3:
UNCHECK 4:UNCHECK 5:UNCHECK 6:¿		Set any existing equations in polar expressions R1 to R6 to zero (there does not seem to be any way to delete them entirely from a program) and deselect them for graphing. This avoids unwanted plots appearing.

Program .Orbit.SExplanation
Called when START is selected. Initialises the display and, because of the SETVIEWS command used, brings up the VIEWS menu when it returns.
RUN ".Orbit.R":¿Just calls the reset routine and returns.

Program .Orbit.4TExplanation
Draws plots corresponding to circular, elliptical, parabolic and hyperbolic orbits. Once the program has run the user is in the PLOT view.
RUN ".Orbit.R":¿Make sure the display is correctly set up.
'A/1'ÞR1(q):CHECK 1:¿The equation for a circle. Just storing 'A' into R1 does not work as the current value of A is inserted instead. Hence use of 'A/1'. The equation is selected.
'A/(1-0.5COS(q))'ÞR2(q):CHECK 2:¿The equation for an ellipse of eccentricity 0.5.
'A/(1-COS(q))'ÞR3(q):CHECK 3:¿A parabola
'A/(1-1.5COS(q))/
((q³ATAN(Ö(1.5²-1)))*
(q£360-ATAN(Ö(1.5²-1))))'
ÞR4(q):CHECK 4:¿		The equation for a hyperbola. The tests are used to prevent the arms being drawn within the range where the radius of the hyperbola is in reality infinite. Without the tests a second mirror-image hyperbola is drawn.
1ÞA:-6AÞXmin:6AÞXmax:
-3AÞYmin:3AÞYmax:
AÞXtick:AÞYtick:2Þqstep:¿		The 'size' of all the curves is set to 1 unit and suitable axes ranges are defined so that the plots fit on and are 'square'. A step size of 2 rather than 1 speeds up plotting with no loss of resolution.

Program .Orbit.ECCExplanation
Draws four ellipses with different eccentricities. The plots are scaled using A*(1-eccentricity²) so that they all have the same length of semi-major axis. Goes to the PLOT view after running.
RUN ".Orbit.R":¿Make sure the display is correctly set up.
'A*(1-0.1²)/(1-0.1COS(q))'
ÞR1(q):CHECK 1:¿		An ellipse with eccentricity 0.1 and semi-major axis=A. Select the equation so it is plotted.
'A*(1-0.3²)/(1-0.3COS(q))'
ÞR2(q):CHECK 2:¿		An ellipse of eccentricity 0.3
'A*(1-0.7²)/(1-0.7COS(q))'
ÞR3(q):CHECK 3:¿		Eccentricity 0.7
'A*(1-0.9²)/(1-0.9COS(q))'
ÞR4(q):CHECK 4:¿		Eccentricity 0.9
1ÞA:-2AÞXmin:3AÞXmax:
-1.3AÞYmin:1.2AÞYmax:
AÞXtick:AÞYtick:2Þqstep:¿		The semi-major axis of all the ellipses is set to 1 unit and suitable axes ranges are defined so that the plots fit on and are 'square'. A step size of 2 rather than 1 speeds up plotting with no loss of resolution.

Program .Orbit.KP2Explanation
Produces an animated illustration of Kepler's second law, that radius vectors sweep out equal areas in equal times, by drawing two pairs of radius vectors. The two members of each pair are the same distance apart in terms of orbit time (1/12 of an orbit) and the two pairs are half an orbit time apart. They are redrawn every 1/36 of an orbit to give an animated effect showing that the area enclosed between each pair of radii remains the same as they rotate..
ERASE:¿Clear the display
-1.3ÞXmin:4.7ÞXmax:
-1.2ÞYmin:1.8ÞYmax:0.5ÞE:¿		Set up a suitable 'square' scale for the display. Set the eccentricity to 0.5.
DISPXY -1;1.8;2;"Kepler's 2nd Law":
DISPXY 2.5;1.2;1;"Equal areas":
DISPXY 3.2;0.8;1;"in":
DISPXY 2.5;0.4;1;"equal times":¿		A brief description of what the animation will show.
0Þq:RUN ".Orbit.KXY":¿Get initial points for the orbit plot.
FOR Z=0.1 TO 360.1 STEP 9;¿Go round the orbit in fairly large steps to speed up plotting. Display is still accurate enough. Starting at 0.1 instead of 0 is to avoid trying to calculate arctan of 90 degrees, which would give an error.
XÞA:YÞB:ZÞq:RUN ".Orbit.KXY":¿Store the X,Y coordinates of the previous point and calculate the coordinates of the next point around the ellipse.
LINE A;B;X;Y:¿Draw a line from the last point to the current point.
END:¿Orbit is plotted.
-9.9Þq:RUN ".Orbit.KXY":XÞN:YÞO:
q+30Þq:RUN ".Orbit.KXY":XÞP:YÞQ:¿		Calculate and store initial values for a radius of the orbit just before the start and a further 30 degrees round.
170.1Þq:RUN ".Orbit.KXY":XÞR:YÞS:
q+30Þq:RUN ".Orbit.KXY":XÞT:YÞU:¿		Calculate initial values for a radius of the orbit half way round and a further 30 degrees round.
FOR Z=0.1 TO 360.1 STEP 10;
NÞF:OÞG:PÞH:QÞI:
RÞJ:SÞK:TÞL:UÞM:¿		Move round the orbit in steps of 10 degrees and as before avoid hitting exactly 90 degrees. Store the coordinates of the points where the previous set of radii met the orbit.
TLINE 0;0;F;G:TLINE 0;0;H;I:
TLINE 0;0;J;K:TLINE 0;0;L;M:¿		Draw in the two pairs of radii.
ZÞq:RUN ".Orbit.KXY":XÞN:YÞO:
Z+30Þq:RUN ".Orbit.KXY":XÞP:YÞQ:¿		Calculate the coordinates for the next set of radii, 30 degrees apart.
Z+180Þq:RUN ".Orbit.KXY":XÞR:YÞS:
Z+210Þq:RUN ".Orbit.KXY":XÞT:YÞU:¿		Calculate the ends of a second pair of radii 180 degrees further round.
TLINE 0;0;F;G:TLINE 0;0;H;I:
TLINE 0;0;J;K:TLINE 0;0;L;M:¿		Draw over the previously drawn radii again using the 'toggle bit' option, and so delete them.
The parameters for the next four lines to be drawn have already been calculated and stored so that on the next pass round the loop they can be drawn immediately and there is a minimum of time when no lines are visible.
END:WAIT 2¿Keep going for a full circle, then pause briefly before returning (normally to the VIEWS menu).

Program .Orbit.KXYExplanation
A subroutine which converts a given angle q (called the mean anomaly around an equivalent circular orbit to the angle measured from the focus of a planet in an elliptical orbit after the same fraction of an orbit.
This is done by solving Kepler's equation for eccentric anamoly E, found from E - e*sinE = q, where e is the eccentricity and all angles are in radians. The eccentric anomaly is the angle to the true position of the planet as seen from the centre of the ellipse. E can only be found via numeric approximation.
The true anomaly V, which is the angle to the planet from the focus, is then calculated from
tan(V/2)=Ö((1+e)/(1-e))*tan(E/2).
Next the distance of the planet from the focus for this angle of true anomaly is found and the resulting polar coordinates converted to X,Y cartesian coordinates.
FNROOT(DEG®RAD(X)-E*SIN(X)-
DEG®RAD(q),X,q)ÞX:
2*ATAN(Ö((1+E)/(1-E))*TAN(X/2))+180ÞV:¿		Uses the Root Finder command to numerically find a solution to Kepler's equation for the eccentric anomaly for a given mean anomaly. Angles have to be in radians for this to work.
This is converted to the true anomaly and stored in V.
180 degrees is added because true anomaly is usually measured from perihelion direction whereas the rest of the program uses measurements from aphelion.
(1-E*COS(V))-1ÞW:¿The length of the radius vector for this true anomaly is calculated, assuming the semi-major axis=1.
W*COS(V)ÞX:W*SIN(V)ÞY:¿The length of the radius and angle are converted to X,Y coordinates.