Showing posts with label longitude of perihelion. Show all posts
Showing posts with label longitude of perihelion. Show all posts

Saturday, 13 November 2010

Week 6 - Initial Coding - Orbital Calculations

Coding has begun! This week I have set up the initial class framework in visual studio (I have switched to visual studio 2008 as this is what is used in the University labs, in case anything happens to my laptop), and begun programming the orbital calculations class.

Firstly I'd like to mention some initial problems I have come across and what I have done to fix them, and changes I have had to make to the UML plan.

In the CreatePlanet class, i have added longitudeAscendingNode and longitudePerihelion as private instance variables (doubles). This was an oversight, as these values are required in the calculations for x1, y1 and z1. I have also removed periodOfOrbit as a variable in CreatePlanet, and moved it to the OrbitalCalculations class. I found it was easier to calculate the period of orbit using the semi-major axis, instead of hard code in a value for each planet.

The bulk of changes has occured in the OrbitalCalculations class. Firstly, meanAnomaly has been added as a private instance variable. I'm not sure why I left this out, but the variable is required to calculate the mean anomaly in the drawOrbit method. Secondly, I found I had neglected to consider the need of parameters in the majority of the methods. This was a mistake as of course most of the calculations utilise data which will, in the end, be taken from the created planets. These changes are as follows:



As you can see from the image, I have also renamed the vector class to Vector3D. I was a bit worried about vector being a reserved C++ word, so changed it to be on the safe side. Another change is to make the drawOrbit method Vector3D, instead of void. This simply means that the method returns the coordinates in vector form, instead of having another method to do this.

One of the major issues I knew I would have is calculating the eccentric anomaly. I tried out a few iteration methods that didn't work, but managed to find C# iteration code on this website: http://www.jgiesen.de/kepler/kepler.html. I adapated this code into C++, and tested it to give the same results as the calculator on that website. A screenshot of this method is below:




The final drawOrbit method utilises all the previous methods in the class to calculate the final x,y,z coordinates, here is a screenshot of this method:


I have also included test printing in this method, which will print the various stages to the console when the method is called in DrawSolarSystem. I had an issue with printf, the printing mechanism I have previously used in C++, and so decided to use cout as my test printing method throughout the program.

I had one final major oversight at this stage, which became apparent when I went to create the planet and test the drawOrbit method in DrawSolarSystem. The mean anomaly calculation requires a current time and a time since the planet was last at the perihelion point to be defined - what format was that in? I hadn't previously considered this fully. With some research, I discovered that astronomers use the Julian dating system (defintion), and thus my program would have to have some system of converting standard (Gregorian) dates into Julian dates. What I did was to add another class, called DateTime. In this class a date with a year, month and day can be defined, then a method called to convert it into a Julian date. It also has a method to convert a Julian date back into a Gregorian date, should it be necessary. The calculations for these methods were adapted from the following website: http://www.usno.navy.mil/USNO/astronomical-applications/astronomical-information-center/julian-date-form. I also used this website which calculates Julian dates to confirm that the method works correcetly. A screenshot of the class is below:


Now i was able to create a planet in the InitGL class of DrawSolarSytem, and call the drawOrbit method of OrbitalCalculations in order to check the test prints. Here are some screenshots of the method, and the resultant console print out:





As you can see, data is given out for each calculation, so it's definately doing something!

Now onto the unresolved problems. As of yet, I haven't been able to find any data to verify and compare the majority of my answers with, aside from the Julian date and obvious ones like period of orbit.

In addition, you may have noticed commented out code to test Mercury as well as the Earth. I've searched for quite some time today trying to find dates for when Mercury was at it's perihelion point, but have so far been unable to find any data. This was easily sourcable data for the Earth, but appears much harder to find for other planets, and is something I will be consulting my supervisor about.

That's pretty much it for this week, some useful links I've found and used this week to follow.

Useful Links:

Perihelion dates for the Earth - dates for when the Earth was/will be at it's perihelion point. Would be lovely to have a similar table for other planets!

http://www.met.reading.ac.uk/~ross/Documents/OrbitNotes.pdf - contains slightly different orbital elements data than the JPL source I have been using. Something to discuss with Rob?

http://www.devmaster.net/forums/showthread.php?t=5059 - How to load 3Ds max objects into OpenGL. Could be useful if my current method fails.

Saturday, 23 October 2010

Blog Week 3 - Orbital Elements and Project Spec

This weeks blog will be relatively short, as I will also be writing a Project Specification, which is due in a week today. Here is an overview of what it will contain:

  •   The project specification is a more detailed description of the work to be done. Two sides of A4 should suffice. 
  • This should document the initial avenue of research / development which the student follows; indicate any prototype developments which are envisaged; suggest comparative studies of different design approaches etc. The initial project plan should be drawn up after a preliminary investigation of the literature. 
  • This deliverable is not formally marked, but it will inform the supervisor’s assessment of conduct. Failure to submit on time (without good reason) will have a negative effect on the conduct grade.

Following my supervisor meeting this week, it became clear that I was looking in the wrong area for the equations I needed to calculate an orbit. What I needed were Orbital Elements. These are 6 elements, some detailed in the diagram below, from this website




The actual 6 orbital elements are:

Eccentricity (e)
Semimajor Axis (a)
Inclination (i)
Longitude of the ascending node (Ω)
Argument of periapsis (\omega\,\!)
Mean anomaly (M)

First, a couple of equations are needed to calculate the planets polar coordinates, (r, v). Kepler's equation allows us to solve the Eccentric anomaly, once we have calculated the Mean anomaly. The two equations are shown below:





Once E has been calculated, r and v can be:









I can make use of this data to calculate the Heliocentric Ecliptic coordinates (x,y,z) for a planet, which are an X,Y,Z coordinate system centered on the Sun in which the ecliptic lies in the X-Y plane. The equations to calculate these points are:



Now, what I really need are the Heliocentric Equatorial coordinates (X,Y,Z) for a planet, in which the coordinate system centred on the Sun with the Sun's equator lying on the X-Y plane - the type I will base my system around. The equations to calculate these points are:





Onto the most important part of the blog - this data is confirmed in a PDF document on the JPL Solar System Dynamics website (a NASA source). The PDF also lists the Keplerian Elements for the time period 1800 AD - 2050 AD, an invaluable resource for me to be able to calculate the planetary orbits for this project. It should be noted that the website states that the data for these positions is only approximate, but it should be efficient enough for this project.