Making an N body gravity simulator in C#
Why?
I wanted to model how to fire something at the earth-moon L1 Lagrange point and have it enter a stable earth orbit. Also for fun.
Goals:
- Use high school maths plugged into C# with some kind of UI that lets me specify a speed, angle, and height to fire a projectile to be caught by the moons gravity and dragged into a stable earth orbit.
- Add frictional simulation to model the effects of the earths atmosphere.
Planned steps:
- Create a Newtonian calculator that returns gravitic forces given object of specified mass and distance. Use temporary unit tests to verify outputs.
- Create a vector calculator that can calculate a final vector given a number of input vectors. Use temporary unit tests to verify outputs.
- Create models to represent bodies including vectors and influences.
- Create a UI and controller classes that can specify bodies and vectors
- Create a calculator class that loops through the specified bodies and calculates new positions given the starting vectors and influences.
GitHub:
Diary:
Step one completed.
Unit tests for Newtonian calculator returning 9.8m/s force +- 0.000...1.
It was very satisfying to see an object with mass 1-100000 feeling a force of 9.8ms at the radius of an object with earth mass. And also seeing the earth mass feel the same force.
Lessons learnt: At first the force was way off, I had missed out the fudge factor "G", which wasn't present in the first web page I looked up. This got me thinking about what the G constant actually means and why its needed. Guesses: 1. Newtons calculations are a simplification of relativity. Would G be different for objects moving nearer the speed of light? 2. The output (force) is measured in newtons, which is directly connected to the unit of mass (in this case kgs), but completely unconnected to the distance unit (in this case meters). Some fudge factor will be needed to connect the two.
Step two/three completed
I totally cheated for the vector calculator as I had already written one for a previous project and just robbed that. I've built the vector calculator into the same class as the Newtonian one for now. This is an early optimisation as many of the calculations the vector calculator needs are already being done by the Newtonian one so it made sense to keep them together. Dumped the previous unit tests as they are no longer of use. The vector logic had already been tested in my previous project so didn't bother with any unit tests for that.
My body model has current and future Coordinates and Vector, in addition it holds the mass, radius, an Id, and a collection of influencing vectors.
Step four/five completed.
Again I borrowed bits from a previous project, a simple winform with text boxes to specify parameters and display outputs.
For some reason I've taken to calling my 'fired' object as the 'payload'. It could be described as a 'satellite' but for some reason my brain prefers payload.
My first order of business was to test it with a payload just above the surface and see what happened. I added the ability to step through iterations with a step button so I could watch the calculation results for each iteration. The object accelerated and fell but not quite right.
I probably should've written integration tests but I was too excited for all that. Manual tests ahoy!
Some debugging later and I realised my influencing forces were not being reset after each iteration so the forces were compounding.
In addition my decision to use longs for the x,y of the Coords is now starting to seem like a really bad idea. In fact the vector calculator is using doubles anyway so any benefit of using longs is lost in the calculation. Quick refactor...
Once the calculation looked ok for a simple falling body I plugged in the numbers for the ISS. 408km altitude, 7.66kms.
At first it looked like it was wrong, the payload seemed to fall when I expected it to stay constant. But the UI display was giving all the figures as pure integers. I realised this made interpreting the data hard so I made all the data fields go though a formatter that added commas, and changed the units when past 1000. Much easier to see whats going on. In fact, the payload is dropping a little, but then going back up. The orbit is slightly elliptical. Not sure why. Maybe its a rounding error somewhere. But in any case the orbit is only elliptical by 20ish km out of 400 so I'm happy enough.
I decided to add a T-minus counter to the output to see how much time has been simulated. This also allowed me a quick guestimate of the orbital period. This was by manually pausing the simulation after once orbit so wasn't an accurate measurement but came out at 92.5 minutes. Quickly checking google for the actual orbital period... 92.65 minutes. Wooo. happy with that :)
 |
The snapshot shows the simulator after 14 days 9 hours, 38 minutes and 41 seconds. The payload is still orbiting +-20km from the initial 408k. Also in this screenshot the moon has been added :)
Adding the moon.
With the excitement of the ISS behaving as expected I was eager to add the moon. Unfortunately the earth did not behave as expected. It accelerated towards the moon much faster than anticipated and went much further before changing direction. In short the 'wobble' the earth feels from the moon was far too high.
It turns out I had made the school boy error of putting a degrees value into a radians field. The moon was moving in a highly elliptical orbit which was yanking the earth around. Once fixed however there was still too much wobble. Not nearly as much but over 7000 km whereas it should be around 4000. I began to wonder, "Is this what happens if you teleport in a (correctly moving) moon next to a stationary earth?" So I checked what speed the earth moves around the barycentre, 12.5ms, then I thought about which direction it should move, the opposite to the moon. Plugging these values into the startup has mostly fixed the wobble.
However now my payload's orbit is too elliptical. Of course, if I'm starting the simulation with a moving earth, I need to add the same vector to any vector I start my payload with.
I decided the easiest way to do this is add the earths current vector to the payloads influences at simulation start. A little refactoring and all is good again (I had to move the influence clear-down to after each iteration vs the beginning,).
There is still too much wobble in the earth though. The speed is staying around 12.5 whereas before it was varying a lot more, however as can be seen in the above screenshot it still gets too high. The screenshot shows the moon just after 14 days. Its moved from directly above (the simulator is in 2D as viewed from above) the earth to almost directly below (just past, moving clockwise). The earth is a little bit to the left but mostly up, which is good as the earth and moon should stay opposite each other in respect to the barycentre. It's just it's about 5000km too 'up'. I've also noticed the moons height fluctuating regularly over a period of minutes/hours rather than days. This is odd, and I've no idea why. It could be precision loss from the double data type used by all the internal calculation.
Update: I figured it out whilst trying to get to sleep. I'm starting the earth at point 0,0, which would make its barycentre x=0, y=4671. This means after half a moon orbit (14 days) the earth should be at x=0, y=9342. Its slightly above that but its close enough for me to rest easier. I could make the barycentre closer to 0,0 by starting the earth at 0, -4671. I would then also have to translate the payload and moons' starting coordinates. Its probably easiest to leave it for now.
Todo:
Much.
- Graphics! I have a big blank canvas below the text output. I need to write a transform to that I can draw the positions and trace out the movements visually.
- Lagrange tests! I want to try starting the payload in each of the L points to see if it can stay stable. L3 and L5 should be stable for long periods of time. The others should be stable-ish.
- Atmosphere! Try to find resources on how the density of the atmosphere changes with height, add calculations for drag and coefficient etc.
Also:
The simulation is quite slow, there are lots of optimisations to be done.
- Firstly I should decouple the UI updates from the intervals. The CPU is churning out a lot of intervals per second, way more than the 60fps the screen can show. This is very wasteful as UI updates are costly.
- In addition with some refactoring all 3 bodies could have their calculations done in parallel. I suspect this wouldn't speed things up as much as fixing the UI however.
- Profile! Many would argue you should do this first, but experience can shortcut you to the biggest gains. Like knowing that calling into a DB or making an HTTP once with a big-ish data packet will be faster than making thousands of smaller calls. The other thing is, how fast do I need it to go? For testing Lagrange point stability many moon orbits would be nice, but just for the primary objective only a few days of simulation is needed. Therefore I'm going to avoid heavily optimising the code until I need to.
Comments
Post a Comment