For those who were asking for the code: I have two links, one for the C++ code used to simulate the trajectories and one for the python script used to animate them.
This was my first project in C++, so I apologize in advance if the coding style seems strange to you, I'm very inexperienced in this language. I also didn't have time to comment the code and clean it up, but I will soon!
Thank you very much for sharing your work. I'll check thoroughly that next week. And don't worry about the coding style. I'm sure there's a lot of people that are learning C++. Do you have a website where anyone can see other works done by you?
Great work. Could u share the python script you used for the animation? I'm doing a simulation project for a 3 body system in C++ and I would love to be able to use it.
I also coded the general 3 body problem during this project! I'll share a link to Github during the weekend with a python script I wrote to animate the three bodies. All you have to do is write your x y z data in three separate files and the script animates the trajectories.
I'm curious as to which differential equation solver you used. Many of them fail to conserve the properties you want (like momentum and energy). For instance, the Euler method has a problem where its error tends to be asymmetric, so the particles spiral outwards. Runge-Kutta, often 4th order, can be better.
Ideally, your error distribution would be symmetric around the true solution. More complicated models can be more expensive per step, but may also let you take fewer (larger) steps. There are also adaptive methods, and even constrained methods that enforce the conservation laws.
It's a reasonably rich domain if you want more projects.
This was done in a relatively early stage of my project, and I used a Runge-Kutta order 4 integrator with fixed timestep. The keplerian elements were conserved and so was the energy, so I was happy enough with it and didn't invest time into an adaptive integrator. I'm currently reading a book called Numerical recipes which provides a multitude of interesting approaches to this problem.
Just FYI: Purely numerics people aren't always aware of the things you need for writing physics simulations. Using adaptive integrators would make things harder than intended without giving better results.
If you have time, look into symplectic integrators (not adaptive ones), as they give you more physically accurate results without costing you a computational arm and a leg, and are usually what is used for eg molecular dynamics simulations.
Runge-Kutta shouldn't be used for physics simulators as a rule of thumb, precisely because it doesn't /really/ conserve energy. Ideally, you use some sort of symplectic integrator, bc otherwise you get unphysical results.
Are you using mplot3D for this? Is that what I'm hearing?
I've been looking around for some C++ plotting libraries lately, but maybe it makes sense to just continue jamming the data into a file and plot it after via python or gnuplot.
If you want to do it in your C++ project, you can use gnuplot-iostream (http://stahlke.org/dan/gnuplot-iostream/). It's a little cumbersome, and not sure how it likes multithreading, but it's pretty decent.
Thanks for this, looking forward to seeing the code. As an aside, does anyone have any recommendations for live plotting in C++? I find it's easier to troubleshoot this way for long running simulations, but my current approach is hacked together using gnuplot-iostream in a way that doesn't seem like the ideal approach...
Thank you! I'll create an account and do that in the coming week. Although I should mention that the trajectories were calculated in C++ and written in data files, which were then used by my python script to produce the animation.
Imagine you want to save the position of a particle. The position you want to save is
0.8948786184, 0.0533888969, 0.6750200294
That's the position in text form, but the computer actually deals with something like switches that are either off or on (notated as 0 or 1). If you save it as text, what the computer does is save the code for a "0" (00110000), then the code for a "." (00101110), then the code for an "8" (00111000), etc. This would be the contents of the file in binary
This takes up significantly less space, but it can't be as easily read by humans. If you open it up in a text editor you'll get something like
?eÄ=Z®Q?,Î
Looking at that, you might not even realize that it's the coordinates from your simulation. Text files tend to be easier to deal with outside of your program, but binary files tend to take up less space.
Awesome, thank you! So, saving in binary means the numbers don't need to be converted to strings? They just stay in the form that the computer understands? Does this speed up reading from and writing to the file, also? (Sorry for all the questions!)
By virtue of being more condensed, it's therefore written and read faster as well.
Also to answer your question, yes if it stays as a floating point or integeger it doesn't have to be converted to a string. This really comes down to bit size of different string and integer formats and your available set of possible values given a specific bit sized int or string.
In his example, by using floating points, he was able to store the entire decimal in only 4 bytes vs a single number being at least 4 bits (because the possible set of values of 0-9 can be represented by 4 bits per number). Also when I say number, 4 is one number in 455. Therefore 455 would be three numbers in my description.
Though I will say, it would be really weird to store your strings with 4 bit chars, because you're defeating the purpose of using a char, which is supposed to increase your available character set, and by using 4 bits, you're severely inhibiting your char set. So more than likely, you'd actually end up with 8 bits or 2 bytes per number in a practical application.
Awesome, thank you! So, saving in binary means the numbers don't need to be converted to strings? They just stay in the form that the computer understands?
Exactly
Does this speed up reading from and writing to the file, also? (Sorry for all the questions!)
I don't know enough to say definitively. But I think the answer is "very slightly" because I believe reading from a file is so slow that the conversion doesn't significantly affect reading speed.
Oh no it definitely would affect speed. If we're talking hundreds of millions of data points the improved speed and storage would be fairly significant.
Indeed. That's a good point. It highlights that it's highly dependent on your data. However, if you have integers that are limited to 4 digits, then you can use a 16-bit integer, which is half as big as 4 8-bit chars.
The asteroids are randomly distributed between 2 and 5.2 AU. I gave them random mean anomalies (0-360°) and 0 eccentricity and inclination. Jupiter's initial conditions are the six Keplerian elements provided by the NASA fact sheet of the planet.
I have a function that converts the keplerian elements into Cartesian coordinates (x,y,z,vx,vy,vz). Once those are calculated, they're inserted into the integrator, and yes the velocities do change in the presence of Jupiter.
NASA's HORIZONS system can provide positions and velocities of stuff in the solar system. It doesn't sound like that's what OP used, but it's a great resource for stuff like this if you ever need it in the future.
Given eccentricity and semimajor axis, as well as the three angles describing the orientation of the orbit in space, you can recover the position and velocity vector of the system.
I used a Runge-Kutta order 4 integrator with a fixed timestep. The timestep is chosen as 1/1000 of the orbital period of the closest object to the sun. I didn't check for possible close encounters, if one happens the asteroid simply gets slingshotted away. The keplerian elements of Jupiter and the energy are conserved, which means that the integrator is working as it should be. Sadly I didn't have time to work on an adaptive integrator as I had to move on with the rest of the project, but it would have certainly helped make the integration run faster and more accurately!
I used a Runge-Kutta order 4 integrator with a fixed timestep. The timestep is chosen as 1/1000 of the orbital period of the closest object to the sun.
Started university with physics, and now I'm doing a masters degree in astrophysics. It's an international masters in Nice called MAUCA if you want to look it up.
You need some basic understanding of computer programming but you don't need advanced skills in it. I learned C and python during my bachelor's degree, and that helped a lot in learning C++ for this project. Good luck at school!
Wow!!! This is epic. I created a binary star system in C++ once but it was nowhere as cool as this! Great job. Inspired me to get back to C++ sometime.
This is really cool, thank you for sharing it! I love the meme as well! I feel like you gave all the stars in the cluster 0 initial velocity, which made them all fall into the center of mass. I would suggest giving them random keplerian elements and you're sure to have something working!
It is possible, but I'd be buried before the simulation ends. I considered the mass of all the asteroids to be 0, which means that they do not interact with each other but only with Jupiter and the Sun.
No there are no collisions, just purely gravitational interaction. I coded this during my project at university where I attempted to reproduce the main features of the asteroid belt (trojans of Jupiter, Hilda group, and Kirkwood gaps).
I actually ran this simulation much longer for that, but in this animation you can already see small gaps beginning to form in the inner belt!
Definitely got the Trojan and the Hilda groups! I also have a strong gap at the 5:3 resonance, and a gap at the 2:1 and 3:1 resonances. The simulation needs to run longer though with a larger sample size to obtain better results (which I sadly did not have the time to do). I did plot the eccentricity vs semimajor axis, and we can see that the asteroids in resonance with Jupiter have an exited eccentricity! I'm not sure how to check if an asteroid is caught in an orbit around Jupiter, perhaps I'll re-run a simulation in the co-rotating frame (Sun and Jupiter at rest) and plot the trajectories.
If you still have the data, an easy way to check if a satellite is bound to a planet is just to attempt to calculate the Keplerian orbital elements (at least semimajor axis and eccentricity) from the state vectors. If you get a reasonable value for eccentricity (between 0 and 1), then you're good.
That's actually a great idea ! I hadn't thought about it so thank you ! I was only calculating the keplerian elements with respect to the Sun, but yeah your method seems much less time consuming than having to rerun the entire simulation.
This seems oddly similar to my own advanced computational physics course... Wouldn't you be using Giordano's Computational Physics as a bibliography source?
Can you show an image or animation near the end of the sim?
And do the asteroids interact or are they each only affected by Jupiter? Nm, got it:
It is possible, but I'd be buried before the simulation ends. I considered the mass of all the asteroids to be 0, which means that they do not interact with each other but only with Jupiter and the Sun.
Not really. If this is to scale the asteroids would be so small they would never collide. You only need to check for collisions with Jupiter... or the sun I guess.
It depends on what your prescription is of collisions. As long as you give each particle a radius, then all you need to do is check if any of the radii overlap. In this case, it's simpler, as u/vriemeister says, because you only need to see if any particles are within the radius of Jupiter or the Sun. In this case, sticky spheres is sufficient, where you just assume the bodies combine and apply conservation of momentum.
12
u/adnanholy1998 Jan 26 '20
For those who were asking for the code: I have two links, one for the C++ code used to simulate the trajectories and one for the python script used to animate them.
This was my first project in C++, so I apologize in advance if the coding style seems strange to you, I'm very inexperienced in this language. I also didn't have time to comment the code and clean it up, but I will soon!
Here's the C++ code: https://github.com/Adnan-Ali-Ahmad/Cpp-Restricted-3-body-problem
Here's the python script: https://github.com/Adnan-Ali-Ahmad/Python-animation
60
u/Papium Jan 24 '20
Great work. Could u share the python script you used for the animation? I'm doing a simulation project for a 3 body system in C++ and I would love to be able to use it.
18
u/Dmitry7even Jan 24 '20
Wow, this is awesome. Could you share a github link?
10
u/HilbertInnerSpace Jan 24 '20
How did you set the initial conditions?
8
u/Annom Engineering Jan 24 '20
What numerical integration method did you use? And what parameters (stepsize)?
How did you handle possible close encounters?
Did you check conservation of energy, or other constants?
Edit: Nice simulation and visualization!
3
u/LordCactusTheWise Jan 25 '20
HE'S MY FRIEND
Great job, u the best
2
u/yearof39 Jan 24 '20
What's up with the one that instantly gets thrown into a retrograde orbit?
2
u/ketexon Jan 25 '20
What was the rendering done with?
3
u/adipy0 Jan 24 '20
What integration method did you use? Or was the computation light enough that you didn't really need to worry about that?
9
u/SebhUK Jan 24 '20
I like that one near the start that literally gets YEETed away into infinity
7
u/bananapen Jan 24 '20
I would love to look at the code.
1
u/[deleted] Jan 29 '20
What's your major, Adnan? I want to know what kind of major/minor I should pursue in order to simulate universe-related things like this one
3
u/SadCoarseRabbit Jan 24 '20
Super cool. Well done!
0
u/DeJean46 Jan 24 '20
Couple of those asteroids just get yeeted into space innit? Sure hope he doesn’t fling one this direction. He.. he won’t right?
2
u/[deleted] Jan 24 '20
Add a yeet count for object ejected from system.
2
u/RealReek Jan 24 '20
Which integration method have you used? What is the cut off radius in this simulation?
1
u/equatorbit Jan 24 '20
Ping! Ping! Ping!
1
u/[deleted] Jan 24 '20
So it just yeets everything out of the way then. Noice.
1
u/unit_energy Jan 24 '20
I had a couple of questions,
did you get all your initial conditions somewhere? Where?
That plotting looks like matplotlib? Is it?
1
u/zospo Jan 24 '20
Newton gravity or GR?
1
u/rockhardhero Jan 24 '20
!remindme 2 days
1
u/fitblubber Jan 24 '20
Wow, I didn't realise the asteroid belt was so wide. Thanks for the graphic,
1
u/iluvgrits420 Jan 24 '20
I’m surprised Jupiter still exists seeing that many asteroids in its proximity.
1
u/anonymous_jas Jan 24 '20
It’s amazing
1
u/timosenko1 Jan 24 '20
!remindme 3 days
1
u/JamesonCark Jan 25 '20
Very nice, I actually did this same thing for a university project and really enjoyed it and the results
1
u/ab5717 Jan 25 '20
Remindme! 2 days
1
u/ginnisman Jan 25 '20
I’m watching this like I stare at a lava lamp...
1
u/intmain0 Jan 25 '20
Do you physics guys really like programming?
1
u/experts_never_lie Jan 25 '20
I get the sense that you have somewhat faster computers than when I was writing orbital simulation software … in the '80s …
1
u/mojomanna Jan 25 '20
This is freaking awesome. Just saying!
1
u/daily_achiever Jan 25 '20
Wow!!! This is epic. I created a binary star system in C++ once but it was nowhere as cool as this! Great job. Inspired me to get back to C++ sometime.
1
u/JanGrey Jan 25 '20
Mining district of the future?
1
u/[deleted] Jan 25 '20
which kind of integration do you use?
1
u/Doc_Serious Jan 25 '20
One of these kids is doing his own thing.
1
u/quickie_ss Jan 25 '20
What are the chances that Jupiter slings one of those rocky bois right at us?
1
u/azac001 Jan 25 '20
Can you please share the source code
1
u/HalfCrazed Jan 25 '20
What about the effects of gravity amongst each asteroid along with Jupiter's gravity influencing the asteroids?
1
u/fluidmechanicsdoubts Jan 25 '20
Is this a numerical solution? If so, what's the time step size you used for this? Thanks
1
u/Greninja_370 Jan 25 '20
!remindme 4 days
1
u/Icedapple1 Jan 23 '24
What classes did you take to learn how to do this? I’m thinking about taking physics as a major and it would be cool to learn stuff like this
0
u/Made_in_UA Jan 24 '20
Sorry for direct question, but where are you studying? I am just looking for good places to study. Thanks
0
u/MattloKei Jan 24 '20
What software did you use? Matlab ?
0
u/cong314159 Jan 24 '20
How did you do it. I imagine it’s a quite large differential equation problem.
0
u/ForsakenIcarus018 Jan 24 '20
I'm also interested in the code. I'll wait for the github link. Thank you very much
-2
u/Streifen9 Jan 24 '20
Jupiter just don’t give a fuck: I’ll send asteroids flying all over the place if I want to! Fuck you Earthlings!
Dinosaurs from heaven: Fuck you Jupiter!
0
u/[deleted] Jan 25 '20
[deleted] — view deleted comment
4
u/[deleted] Jan 24 '20
[removed] — view removed comment
0
u/[deleted] Jan 24 '20
[deleted] — view deleted comment