So I decided to make a solar system...

[Eriksonn]

I have been playing around recently with some n-body gravity and such, and i made this nifty thing. It starts with 800 "asteroids" orbiting a star, and after some time they collide and form planets that then merge to make a small solar system with stable orbits.

Now that i have got it into a mostly finished state, i thought i would share it here to see what people think about it. 

It is in 3d so they have inclination also but i am showing it top-down so it is easy to see.

So, lets get into it:

It starts with 800 asteroids around a star.

Then the first planets form:

When more planets form, they start pulling on eachother and some get unstable. If you look closely some have much higher eccentricity than the rest, like that purple one in the center (he collided soon after).

After alot more timewarping(it goes faster with fewer objects), they settle down into 3 planets, with a few "comets" on the outside. (You might need to magnify the image to see them)

It usually end up with 3-4 planets in the end, but this one almost had 5 until the 5th one got wonky and messed everything up.

 

Edit: Here is the code for anyone who wants to play around with it. You are free to ask questions about any problem you might have.

https://processing.org/download/

I couldnt figure out how to put the actual file here so i did this instead. If anyone knows how to put files here, let me know. Just make a new processing program and copy this in there and it should work.

ArrayList<Planet> Planets=new ArrayList<Planet>();
int PathLength=60;
int pathStep=20;
Sun sun;
float LineLimit=0.3;//only draw lines when mass is above
float G = 0.1;
boolean SpedUp=false;
boolean HasPressed=false;
int Warp=1;
PGraphics G3D;
PGraphics G2D;
float rotY=0;
float rotX=0.5;
float cameraSize=4;
int OldSize=0;
float Year=0;
float YearSpeed;
void setup()
{
  size(800,800,P3D);
  G2D = createGraphics(90,70);
  G3D = createGraphics(width,height,P3D);
  OldSize=width+height;
  colorMode(HSB);
  sun = new Sun();
  
  //you can manually add planets like this if you want, and remebmer to comment out the cloud function
  //Planets.add(new Planet(new PVector(sun.Pos.x-550,sun.Pos.y),new PVector(),0.5));
  
  GenerateCloud(sun);
}
void GenerateCloud(Sun sun)
{
  float R=1000;
  float a = (R+R*0.1)/2;
  YearSpeed=1/(PI*sqrt(a*a*a/(G*sun.Mass)));
  println(YearSpeed);
  int Amount=800;
  float Ratio=600;//sun-disc mass ratio
  float SolarMass=sun.Mass;
  float AverageMass=SolarMass/(Ratio*Amount);
  println(AverageMass);
  for(int i =0;i<Amount;i++)
  {
    float r=R*sqrt(random(0.1,1));
    float A=random(2*PI);
    float Mass = AverageMass*random(0.7,1.3);
    Planets.add(new Planet(new PVector(sun.Pos.x+r*cos(A),sun.Pos.y+r*sin(A),random(-r,r)*0.05),new PVector(),Mass,color(random(255),255,255)));
  }
}
float getGrav(PVector D)
{
  return 1/(D.magSq());
}
float getSize(float Mass)
{ 
  return sqrt(Mass)*2;
}
void draw()
{
  
  
  if(OldSize!=width+height)
  {
    OldSize=width+height;
    G3D = createGraphics(width,height,P3D);
    
  }
  for(int i =0;i<Warp;i++)
    {
      Year+=YearSpeed;
      update();
      UpdateSunGrav();
    }
  G3D.beginDraw();
  G3D.camera(200*cameraSize,0 ,0, 
         0.0, 0.0, 0.0, 
         0.0, 0.0, -1.0);
  G3D.background(0);
  
  G3D.rotateY(rotX);
  G3D.rotateZ(-rotY);
  G3D.strokeWeight(2);
  show();
  G3D.endDraw();
  image(G3D,0,0);
  G2D.beginDraw();
  G2D.background(0);
  G2D.text("Amount:"+Planets.size(),20,20);
  G2D.text("Warp:"+Warp,20,40);
  G2D.text("Year:"+round(10*Year)/10.0,20,60);
  G2D.endDraw();
  image(G2D,0,0);
}
void update()
{
 
 Planet[] Crashed=new Planet[2];
 boolean crash=false;
 for(int i =0;i<Planets.size()-1;i++)
 {
   for(int j =i+1;j<Planets.size();j++)
   {
     Planet P1 = Planets.get(i);
     Planet P2 = Planets.get(j);
     PVector D = PVector.sub(P1.Pos,P2.Pos);
     boolean P1Line=P1.Mass>LineLimit;
     boolean P2Line=P2.Mass>LineLimit;
     if(D.magSq()<sq(P1.R+P2.R))
     {
       Crashed[0]=P1;
       Crashed[1]=P2;
       crash=true;
     }else if(P1Line||P2Line)
     {
       float F = G*P1.Mass*P2.Mass*getGrav(D);
       D.normalize();
       D.mult(F);
       if(P2Line)
       P1.Speed.sub(PVector.div(D,P1.Mass));
       if(P1Line)
       P2.Speed.add(PVector.div(D,P2.Mass));
     }
   }
 }
 if(crash)
 {
   Planet P1 = Crashed[0];
   Planet P2 = Crashed[1];
   int id=0;
   if(P1.Mass<P2.Mass)
   id=1;
   float mass=P1.Mass+P2.Mass;
   PVector pos = PVector.add(P1.Pos,P2.Pos).div(2);
   PVector spd=PVector.add(PVector.mult(P1.Speed,P1.Mass),PVector.mult(P2.Speed,P2.Mass)).div(mass);
   Crashed[id].Pos=pos;
   Crashed[id].Mass=mass;
   Crashed[id].R=getSize(mass);
   //println(Crashed[id].R);
   Crashed[id].Speed=spd;
   Planets.remove(Crashed[1-id]);
 }
 crash=false;
 
 for(Planet P:Planets)
 {
   float KillDist=(width+height)*3;
   PVector D = PVector.sub(P.Pos,sun.Pos);
   float DistSq=D.x*D.x+D.y*D.y;
   if(DistSq<sq(P.R+sun.R)||DistSq>sq(KillDist))
   {
     crash=true;
     Crashed[0]=P;
   }else
   {
     float F = G*P.Mass*sun.Mass*getGrav(D);
     D.normalize();
     D.mult(F);
     P.Speed.sub(PVector.div(D,P.Mass));
     P.Pos.add(P.Speed);
     if(P.Mass>=LineLimit)
     {
       P.PathTimer++;
       if(P.PathTimer>pathStep)
       {
         P.Path[0]=new PVector(P.Pos.x,P.Pos.y,P.Pos.z);
         P.PathTimer=0;
         for(int i = 0;i<PathLength-1;i++)
         {
           int i2=PathLength-i-1;
           P.Path[i2]=P.Path[i2-1];
         }
       }
     }
     
   }
 }
 if(crash)
 {
   Planets.remove(Crashed[0]);
 }
}
void UpdateSingleBody(Planet Main)
{
  for(int i =0;i<Planets.size();i++)
  {
    Planet P1 = Planets.get(i);
    if(P1!=Main)
    {
     PVector D = PVector.sub(P1.Pos,Main.Pos);
       float F = G*P1.Mass*Main.Mass*getGrav(D);
       if(F<0.01)
       {
         D.normalize();
         D.mult(F);
         P1.Speed.sub(PVector.div(D,P1.Mass));
         Main.Speed.add(PVector.div(D,Main.Mass));
       }
    }
    
  }
}
void UpdateSunGrav()
{
  for(int i =0;i<Planets.size();i++)
  {
    Planet P1 = Planets.get(i);
    PVector D = PVector.sub(P1.Pos,sun.Pos);
    float F = G*sun.Mass*getGrav(D);
     //D.normalize();
     D.setMag(F);
     P1.Speed.sub(D);
     P1.Pos.add(P1.Speed);
  }
}
void show()
{
  G3D.strokeWeight(1);
  for(Planet P:Planets)
 {
   G3D.stroke(P.C);
     for(int i =0;i<PathLength;i++)
     {
       if(P.Path[i]!=null)
       {
         G3D.stroke(P.C,lerp(255,20,i/(float)PathLength));
         if(i==0)
         G3D.line(P.Path[i].x,P.Path[i].y,P.Path[i].z,P.Pos.x,P.Pos.y,P.Pos.z);
         else
         G3D.line(P.Path[i].x,P.Path[i].y,P.Path[i].z,P.Path[i-1].x,P.Path[i-1].y,P.Path[i-1].z);
       }
     }
     G3D.noStroke();
     //stroke(0);
     
     if(P.Mass<LineLimit)
     {
       //fill(P.C);
       G3D.stroke(P.C);
       G3D.point(P.Pos.x,P.Pos.y,P.Pos.z);
       G3D.noStroke();
     }else
     {
       G3D.fill(P.C);
       G3D.translate(P.Pos.x,P.Pos.y,P.Pos.z);
       G3D.sphere(P.R);
       G3D.translate(-P.Pos.x,-P.Pos.y,-P.Pos.z);
     }
     
     
 }
 G3D.fill(sun.C);
 G3D.sphere(sun.R);
}
void mouseWheel(MouseEvent event) {
  float e = event.getCount();
  if(e>0)
  {
    cameraSize*=1.2;
  }else if(e<0)
  {
    cameraSize/=1.2;
  }
}
void mouseDragged()
{
    rotY+=(mouseX-pmouseX)/200.0;
    rotX+=(mouseY-pmouseY)/200.0;
    rotX = constrain(rotX,-PI/2,PI/2);
}
void keyPressed()
{
  if(key=='+'&&keyPressed)
  {
    Warp=(int)(Warp*1.3)+1;
  }else if(key=='-'&&keyPressed)
  {
    if(Warp>1)
      Warp=(int)(Warp/1.3)-1;
  }
}
class Planet
{
  PVector Pos;
  PVector Speed;
  PVector[] Path;
  int PathTimer;
  float Mass;
  float R;
  color C;
  Planet(PVector Pos,PVector Speed,float Mass)
  {
    Path=new PVector[PathLength];
    C=color(random(255),255,255);
    R=getSize(Mass);
    R=max(R,2);
    this.Pos=Pos;
    if(Speed.magSq()==0)
    {
      PVector D=PVector.sub(Pos,sun.Pos);
      float R = D.mag();
      D=new PVector(D.y,-D.x);
      D.setMag(sqrt(G*(Mass+sun.Mass)/R));
      this.Speed=D;
    }else
      this.Speed=Speed;
    this.Mass=Mass;
  }
  Planet(PVector Pos,PVector Speed,float Mass,color C)
  {
    this.C=C;
    Path=new PVector[PathLength];
    R=getSize(Mass);
    R=max(R,2);
    this.Pos=Pos;
    if(Speed.magSq()==0)
    {
      PVector D=PVector.sub(Pos,sun.Pos);
      float R = D.mag();
      D=new PVector(D.y,-D.x);
      D.setMag(sqrt(G*(Mass+sun.Mass)/R));
      this.Speed=D;
    }else
      this.Speed=Speed;
    this.Mass=Mass;
  }
}
class Sun
{
  PVector Pos;
  float Mass=36000;
  float R = 25;
  color C =color(50,255,255);
  Sun()
  {
    Pos=new PVector(0,0);
  }
}

 

Edited August 17, 2018 by Eriksonn
Added the actual code here aswell

[Jaelommiss]

Looks interesting. How much can be tweaked between simulations? Can a user start with more or less asteroids, change the range of masses they start with, their initial velocity, etc? Does it give you the option to automatically save a screenshot periodically (say, every year) so you can run a more populated simulation over several days, then go back and watch it?

In any case it looks like the sort of thing I'd spend a few hours playing around with, then leave running overnight to see how much it changes.

Edited August 17, 2018 by Jaelommiss

[Green Baron]

Nice work !

Would you tell us something about the algorithm or the fundamentals you used ?

[Eriksonn]

8 hours ago, Jaelommiss said:

Looks interesting. How much can be tweaked between simulations? Can a user start with more or less asteroids, change the range of masses they start with, their initial velocity, etc? Does it give you the option to automatically save a screenshot periodically (say, every year) so you can run a more populated simulation over several days, then go back and watch it?

The code is not that long and all the intresting stuff is at the top, so finding the values you want to change should not be that difficult. I have variables for the amount of asteroids, their average mass and the size of the cloud and a bunch more stuff you can play around with. I can link the code and you also need processing to run and edit it. (Code liked in original post)

9 hours ago, Jaelommiss said:

In any case it looks like the sort of thing I'd spend a few hours playing around with, then leave running overnight to see how much it changes.

you dont have to wait that long, becauce once the amount goes below 100 the amount you can warp increases very fast. when there is ~10 things left i get the warp to over 10k without problems at all. (~400 years per second)

 

1 hour ago, Green Baron said:

Nice work !

Would you tell us something about the algorithm or the fundamentals you used ?

I just used newtons gravity formula to calculate the forces between the bodies and then moved them according to those forces. And made them collide if they touched each other. Nothing too complex going on here.

The asteroids also dont pull on things to make the fps a bit better at the start.

[Eriksonn]

I did Another run and found something intresting. It ended up similar to the previous one but with a small object in a very inclined orbit but the same radius as the big outer planet. It is rare for those inclined ones to survive for very long.

It is in a horseshoe orbit with the big one. That basically means that it is in a slightly different orbit and alternates between inside and outside compared to the big one. The inner 2 planets are in a 5/3 resonance with each other.(5 blue orbits in the same time as 3 cyan orbits) This makes them stable and dont get messed up by the other planets.

[Green Baron]

I'd like to play around, but i have my own little mess to attend to :-)

But maybe it encourages others here who have asked about programming before. It shows what can be done with just some lines of code.

[sevenperforce]

32 minutes ago, Eriksonn said:

I did Another run and found something intresting. It ended up similar to the previous one but with a small object in a very inclined orbit but the same radius as the big outer planet. It is rare for those inclined ones to survive for very long.

It is in a horseshoe orbit with the big one. That basically means that it is in a slightly different orbit and alternates between inside and outside compared to the big one. The inner 2 planets are in a 5/3 resonance with each other.(5 blue orbits in the same time as 3 cyan orbits) This makes them stable and dont get messed up by the other planets.

It's really neat that you got a 5:3 resonance just out of the starting conditions.

[KG3]

Is it possible for planets in this simulation to capture moons?   I'm not sure how the whole moon capture works but I think it has something to do with the orbital speed of the planet being close to the escape velocity of the planet.  Also it seems that moons in a retrograde obits are considered to have been captured.  Are moons more likely to be captured into a retrograde orbit or do they just stand out more?  

Great job!  Would I have to plot motions and locations for all 800 objects before starting the simulation or have you done that already?     

[Spaceception]

1 hour ago, Eriksonn said:

I did Another run and found something intresting. It ended up similar to the previous one but with a small object in a very inclined orbit but the same radius as the big outer planet. It is rare for those inclined ones to survive for very long.

It is in a horseshoe orbit with the big one. That basically means that it is in a slightly different orbit and alternates between inside and outside compared to the big one. The inner 2 planets are in a 5/3 resonance with each other.(5 blue orbits in the same time as 3 cyan orbits) This makes them stable and dont get messed up by the other planets.

3

That looks like you recreated the outer solar system

[sevenperforce]

Question: did you start all the asteroids in circular prograde orbit?

18 minutes ago, KG3 said:

Is it possible for planets in this simulation to capture moons?   I'm not sure how the whole moon capture works but I think it has something to do with the orbital speed of the planet being close to the escape velocity of the planet.  Also it seems that moons in a retrograde obits are considered to have been captured.  Are moons more likely to be captured into a retrograde orbit or do they just stand out more?  

Suppose you place a body that is completely stationary at the very edge of Earth's gravitational sphere of influence, then give it a tiny nudge inward. Earth's gravity will pull it steadily, and it will impact the atmosphere at escape velocity. If it has a slight tangential velocity when it enters the SOI, then it will swing by perigee at well over escape velocity and then slingshot out. Think about entering orbit around a planet after a transfer: you need a burn to slow down at or near perigee, or you will not be enter orbit. In order for a moon to be captured, some portion of that orbital energy needs to be used up.

This usually happens in one of three ways:

• Collision. If the secondary hits something on its way by, then it will lose energy and be captured. It can hit an existing moon, or it can hit the planet itself. If impact with the primary is more than a glancing blow, the secondary will usually be ripped into shreds and will re-form as a new body (which is what happened when Thea struck proto-Earth, resulting in the formation of our moon). Note that this is equivalent to aerobraking/lithobraking a spacecraft.
• 3-body interaction. Another way a moon can be captured is if it imparts some of its orbital energy to something else. For example, if a dwarf planet swings through a planet's SOI and comes close to an existing moon (but doesn't collide), it can transfer some of its velocity to that moon and thus be captured. This will change the existing moon's orbit, often ejecting it from the system altogether; if not ejected, a future collision between the two is likely. You use this effect in KSP when you capture at Jool by getting a gravity assist off Tylo or Laythe. Rarely, this also happens if the incoming moon-to-be is a tight binary; one member is ejected with excess hyperbolic velocity while the other member is slowed and captured. This is what likely happened for the asteroid-moons of Mars.
• Tidal magic. If you have a sufficiently large secondary passing very close to the primary on its swingby, the second-order tidal effects can warp the shape of the secondary (and sometimes the primary as well), converting part of the orbital energy into rotation and part of it into internal heat. This can sometimes be enough for a capture.

[Eriksonn]

42 minutes ago, KG3 said:

Is it possible for planets in this simulation to capture moons?   I'm not sure how the whole moon capture works but I think it has something to do with the orbital speed of the planet being close to the escape velocity of the planet.  Also it seems that moons in a retrograde obits are considered to have been captured.  Are moons more likely to be captured into a retrograde orbit or do they just stand out more?      

I dont think that moons are possible in this simulation due to 1.The planets are very big compared to their own soi, 2. They are too close to the sun, 3.The timesteps are too large.(ie. the moon "teleports" past the planet before being able to orbit). I also think that it wont work becauce of the order in when stuff happens. When most larger planets form it is in the time when there is lots of chaos everywhere, and so not good time for moons. And in the end when things are stable there are no more tiny asteroids left to be captured anyway.

I managed to get a moon to orbit a planet by spawning them manually, but it was a very narrow range of values that could support it. Sidenote: i dont actually know the masses of the larger planets in the end of my simulations, so i might have chosen lower mass on my planet that what is possible.

42 minutes ago, KG3 said:

  Great job!  Would I have to plot motions and locations for all 800 objects before starting the simulation or have you done that already?     

All objects in the cloud are spawned randomly, but you can manually add new planets if you want to. If no velocity is given they are put on circular orbits.

Edited August 17, 2018 by Eriksonn