Need help with N-body simulation

[megapiano]

TL;DR Summary

I need help modifying my N-body simulation's force calculation formula so that galaxy simulations give rise to a spiral-arm aesthetic.

Many years ago, for my high school senior project, I wrote a solver for the N-body problem that is performance-optimized using the Barnes-Hut algorithm. (the optimization algorithm is not relevant to my question.)

In one particular simulation, I simulated a spiral galaxy. The simulation is not meant to be scientifically accurate, but is instead meant to demonstrate the simulator. You can view a video of this simulation below:


Now, as you can see in the video, the galaxy develops pseudo-spiral-arms. I say "pseudo" because the arms are not developed according to density wave theory. Instead, I incorporated a trigonometric function, such as sin(), in the gravity-force calculation of the particles. Without this addition, "vanilla" N-body simulations of galaxies do not develop the spiral arm structure you see in the video.

Here are a couple screenshots from the video that show that sin() is being used in the force calculation as evidenced by the wave-like streams of particles annotated by the red lines (note that there are N >2 arms in these screenshots. In the video, if you keep watching, the number of these is reduced and the existing arms grow in size.)

The problem: Now, many years later, I have recovered the source-code of the simulator that is an earlier version of what is used in the video. And in the force calculation of this version of the code, there is no sin() implementation in the force calculation, and thus, I cannot reproduce the spiral-galaxy aesthetic that is seen in the video: the pseudo-spiral-arms do not form.

After many attempts I am unable or recreate the argument that I passed to the sin() function or where I even put it in the force calculation.

Here is the force calculation as derived from Newton's law of gravitation. The smoothing variable is used to prevent singularities from forming and also to prevent wild sling-shot effects:

forceGravity = (G * particleMass1 * particleMass2) / (distanceBetweenParticles^2 + smoothing)

forceGravity is a quantity and I calculate the 3d vector to apply to both of the particle's acceleration vectors according to their relative positions.

So, somewhere in there was a sin() function that gave rise to the aesthetic of spiral arms.

Can anyone figure out/know how to re-implement sin() in the force calculation? I think it was used as a coefficient or something, but I'm unsure what I passed as the parameter to the function. Candidates for variables used in the formula passed as the parameter may include: distance to galactic center, expected galactic radius, something to do with the position of the particles relative to the galactic center, or others that I can't think of.

Any help is greatly appreciated.

 

[jedishrfu]

This is very hard to figure out with what little you've given us.

I did find a galaxyapp.java source file from the OpenSourcePhysics collection written in java. While it likely doesn't use the same algorithm, it might jog your memory and you'll remember how you did it.

To run the code you'll need java on your machine and the OSP jar file. More details on using OSP can be found at their website:

https://www.compadre.org/osp

The draw() method uses sin/cos for x,y positioning of the stars and computes the angle for these functions:

  public void draw(DrawingPanel panel, Graphics g) {
    g.setColor(Color.black);
    g.fillRect(0, 0, panel.getWidth(), panel.getHeight()); // color background black
    g.setColor(Color.yellow);
    for(int label = 0;label<numberOfCells;label++) {
      if(starLifeTime[label]>0) {
        int r = cellR[label];
        int a = cellA[label];
        double angle = twoPiOver6*a/r+(v*t)/r; // angle for cell at time t
        double x = r*Math.cos(angle);
        double y = r*Math.sin(angle);
        double ds = starLifeTime[label]/15.0;
        int leftPixels = panel.xToPix(x);
        int topPixels = panel.yToPix(y);
        int heightPixels = (int) (panel.getYPixPerUnit()*ds);
        int widthPixels = (int) (panel.getXPixPerUnit()*ds);
        g.fillOval(leftPixels, topPixels, widthPixels, heightPixels);
        starLifeTime[label]--;                 // decrease star cluster lifetime
      }
    }
  }

Here's the complete source:

/*
 * Open Source Physics software is free software as described near the bottom of this code file.
 *
 * For additional information and documentation on Open Source Physics please see: 
 * <http://www.opensourcephysics.org/>
 */

package org.opensourcephysics.sip.ch19;
import org.opensourcephysics.controls.*;
import org.opensourcephysics.frames.*;
import org.opensourcephysics.display.*;
import java.awt.*;

/**

 * GalaxyApp models evolution of a galaxy using percolation ideas

 * Model proposed by Schulman and Seiden

 * @version 1.0

 * @author J. Tobochnik  3/25/05

 *

 */
public class GalaxyApp extends AbstractSimulation implements Drawable {
  DisplayFrame frame = new DisplayFrame("Galaxy");
  final static double twoPi = 2.0*Math.PI;
  final static double twoPiOver6 = twoPi/6;
  double p, v, dt, t;
  int numberOfRings, numberOfActiveCells, numberOfCells;
  int[] starLifeTime, cellR, cellA, activeCellLabel;

  public GalaxyApp() {
    frame.addDrawable(this);
  }

  public void initialize() {
    t = 0;
    numberOfRings = control.getInt("Number of rings");
    numberOfActiveCells = control.getInt("Initial number of active cells");
    p = control.getDouble("star formation probability");
    v = control.getDouble("cell velocity");
    dt = control.getDouble("time step");
    frame.setPreferredMinMax(-numberOfRings-1.0, numberOfRings+1.0, -numberOfRings-1.0, numberOfRings+1.0);
    numberOfCells = 3*numberOfRings*(numberOfRings+1);
    cellR = new int[numberOfCells];
    cellA = new int[numberOfCells];
    int cellLabel = 0;
    // initial values of r and a for each cell label
    for(int r = 1;r<=numberOfRings;r++) {
      for(int a = 0;a<r*6;a++) {
        cellR[cellLabel] = r;
        cellA[cellLabel] = a;
        cellLabel++;
      }
    }
    starLifeTime = new int[numberOfCells];
    activeCellLabel = new int[numberOfCells];
    initializeGalaxy();
  }

  public void initializeGalaxy() {
    int i = 0;
    while(i<numberOfActiveCells) {
      int label = (int) (Math.random()*numberOfCells);
      if(starLifeTime[label]!=15) {
        starLifeTime[label] = 15; // activate region for 15 time steps
        activeCellLabel[i] = label;
        i++;
      }
    }
  }

  public void formNewStars() {
    // copy list of active cells into another array
    int[] currentActiveCellLabel = activeCellLabel.clone();
    int currentNumberOfActiveCells = numberOfActiveCells;
    numberOfActiveCells = 0;
    for(int i = 0;i<currentNumberOfActiveCells;++i) {
      int cellLabel = currentActiveCellLabel[i];
      int r = cellR[cellLabel];
      int a = cellA[cellLabel];
      // activate neighboring cells in same ring
      createStars(r, pbc(a+1, r));
      createStars(r, pbc(a-1, r));
      // activate cells in next larger ring
      if(r<numberOfRings-1) {
        int ap = aForOtherRadius(a, r, r+1);
        createStars(r+1, pbc(ap, r+1));
        createStars(r+1, pbc(ap+1, r+1));
      }
      // activate cells in next smaller ring
      if(r>1) {
        int am = aForOtherRadius(a, r, r-1);
        createStars(r-1, pbc(am, r-1));
        createStars(r-1, pbc(am+1, r-1));
      }
    }
  }

  public int pbc(int a, int r) {
    return(a+6*r)%(6*r);
  }

  public int aForOtherRadius(int a, int r, int rOther) {
    // angle corresponding to a at time t
    double angle = twoPiOver6*a/r+((v*t)/r);
    angle -= twoPi*(int) (angle/twoPi);
    // change in angle for cell at other r
    double angleChange = ((v*t)/rOther);
    angleChange -= twoPi*(int) (angleChange/twoPi);
    // return value of a for cell at other r
    return(int) ((rOther/twoPiOver6)*(angle-angleChange));
  }

  public void createStars(int r, int a) {
    int label = a+3*r*(r-1);
    if((Math.random()<p)&&(starLifeTime[label]!=15)) {
      activeCellLabel[numberOfActiveCells] = label;
      numberOfActiveCells++;
      starLifeTime[label] = 15;
    }
  }

  public void doStep() {
    t += dt;
    formNewStars();
    frame.setMessage("t = "+decimalFormat.format(t)+" #active = "+numberOfActiveCells);
  }

  public void reset() {
    control.setValue("Number of rings", 50);
    control.setValue("Initial number of active cells", 200);
    control.setValue("star formation probability", 0.18);
    control.setValue("cell velocity", 1.0);
    control.setValue("time step", 10.0);
  }

  public void draw(DrawingPanel panel, Graphics g) {
    g.setColor(Color.black);
    g.fillRect(0, 0, panel.getWidth(), panel.getHeight()); // color background black
    g.setColor(Color.yellow);
    for(int label = 0;label<numberOfCells;label++) {
      if(starLifeTime[label]>0) {
        int r = cellR[label];
        int a = cellA[label];
        double angle = twoPiOver6*a/r+(v*t)/r; // angle for cell at time t
        double x = r*Math.cos(angle);
        double y = r*Math.sin(angle);
        double ds = starLifeTime[label]/15.0;
        int leftPixels = panel.xToPix(x);
        int topPixels = panel.yToPix(y);
        int heightPixels = (int) (panel.getYPixPerUnit()*ds);
        int widthPixels = (int) (panel.getXPixPerUnit()*ds);
        g.fillOval(leftPixels, topPixels, widthPixels, heightPixels);
        starLifeTime[label]--;                 // decrease star cluster lifetime
      }
    }
  }

  public static void main(String[] args) {
    SimulationControl.createApp(new GalaxyApp());
  }
}

/* 
 * Open Source Physics software is free software; you can redistribute
 * it and/or modify it under the terms of the GNU General Public License (GPL) as
 * published by the Free Software Foundation; either version 2 of the License,
 * or(at your option) any later version.

 * Code that uses any portion of the code in the org.opensourcephysics package
 * or any subpackage (subdirectory) of this package must must also be be released
 * under the GNU GPL license.
 *
 * This software is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
 * or view the license online at http://www.gnu.org/copyleft/gpl.html
 *
 * Copyright (c) 2007  The Open Source Physics project
 *                     http://www.opensourcephysics.org
 */

 

[megapiano]

This is very hard to figure out with what little you've given us.

I did find a galaxyapp.java source file from the OpenSourcePhysics collection written in java. While it likely doesn't use the same algorithm, it might jog your memory and you'll remember how you did it.

Can you post a link to that project? I can't find it on the OpenSourcePhysics website. Or maybe just a sceenshot of the simulation?

 

[jedishrfu]

It’s part of examples embedded in the source jar file. It’s related to a chapter in their book on computer simulations.

 

[pbuk]

What is the point in creating animations of "galaxies" with artificially created spiral arms by changing the laws of gravitation? You may as well use a paint brush: we cannot learn anything from these and they are not as beautiful as the real thing.

 

[megapiano]

What is the point in creating animations of "galaxies" with artificially created spiral arms by changing the laws of gravitation? You may as well use a paint brush: we cannot learn anything from these and they are not as beautiful as the real thing.

I was just a high school student when I developed this project. The purpose was to create an N-body solver, and then reduce the complexity from O(N^2) to O(nlog(n)) by implementing the Barnes-Hut recursive tree-building optimization algorithm. This was a computer science project that interested me, not a PhD guided research program. I had to present the project to a couple science classes, and thought a reasonably aesthetic visualization would be good.

These days, I would like to make it more scientifically accurate, but I would also like to recover the spiral-arm aesthetic first.

@jedishrfu
That is a very interesting algorithm that is used. Unfortunately, it is not an N-body solver, and I need to spend more time studying it to see if I can adapt anything from it.

However, I might write an implementation of the that algorithm in C#. Or maybe re-implement it in Java in a different way that will make it easier for me to tinker with and learn from.

If you need more info on how my N-body solver works, it's pretty straightforward if we ignore the optimization algorithm (which we can do because it is not functionally relevant to the physics except during force calculation)

Here is some pseudo-code of the essentials (the actual code is in C#):

//use this force calculation if two particles are attracting
forceGravity = (G * particleMass1 * particleMass2) / (distanceBetweenParticles^2 + smoothing)
//use this force calculation if particle is attracting to a tree-node
forceGravity = (G * particleMass1 * nodeMass) / (distanceFromNodeCenterOfGravityToParticle^2 + smoothing)
//
vec3Gravity = particlePosition1 - particlePosition2
//
particleNetForce += vec3Gravity.normalize() * forceGravity
//
particleAcceleration = (particleNetForce / particleMass)
//
particleVelocity += particleAcceleration * deltaT
//
particlePosition += particleVelcocity * deltaT

 

[pbuk]

I was just a high school student when I developed this project. The purpose was to create an N-body solver, and then reduce the complexity from O(N^2) to O(nlog(n)) by implementing the Barnes-Hut recursive tree-building optimization algorithm. This was a computer science project that interested me, not a PhD guided research program. I had to present the project to a couple science classes, and thought a reasonably aesthetic visualization would be good.

That is a challenging project for a high school student, and it looks like you executed it well from that point of view.

These days, I would like to make it more scientifically accurate

Yes, as an adult with that experience behind you you should be aiming higher.

but I would also like to recover the spiral-arm aesthetic first.

That will act in exactly the opposite direction of scientific accuracy.

Here is some pseudo-code of the essentials (the actual code is in C#):

That code uses Euler's method and also adds a constant acceleration to each step, both of which make the method non-conservative. I suggest your first move should be to address this by using a symplectic method and dealing with close approaches in a different way (what is the problem with close approaches anyway: stars being ejected from a galaxy is a real thing? And in 3D close approaches are a very rare event: you are working in 3D?).

The first test of the improved accuracy of your method should be ensuring it maintains a stable 2-body orbit for many rotations.

If you are comfortable in C# then Java would be a step backwards: I cannot see any good reason to go in this direction.

 

[megapiano]

That code uses Euler's method and also adds a constant acceleration to each step, both of which make the method non-conservative.

Sorry, how do you mean? The acceleration vector is re-calculated every step according to the netforce vector, which is reset and then cumulates gravitational forces at every step.

Or do you mean that the acceleration vector doesn't have a deltaTime coefficient? In the actual code, it does, but I removed it in that post because it didn't seem right; I'm not super sure why I did that, because I remember testing out simple scenarios when I originally wrote the solver and the computation was correct according to the math I tested it against

what is the problem with close approaches anyway: stars being ejected from a galaxy is a real thing? And in 3D close approaches are a very rare event: you are working in 3D?

Yes, the simulator is in 3d. The smoothing value is added to prevent the sling-shot effect in all simulations, not just that particular one, which may or may not suffer from the effect. Whether or not the effect occurs in a simulation is entirely dependent on the initial conditions and parameters (G, unit scaling, etc.) of whatever the current simulation is. Adding a smoothing value is pretty standard for n-body solvers.

The first test of the improved accuracy of your method should be ensuring it maintains a stable 2-body orbit for many rotations.

Yes, it can do this. And, as you know, the smaller the timestep, the better.

If I recall correctly, smoothed particle hydrodynamics (SPH) is a technique used to help account for the effects of dark matter in simulations of large enough scale, not super sure if I'm remembering this correctly. I was an undergraduate research assistant in the astrophysics department working on computational models of asteroids, and SPH is a technique used there.

Also, note that that simulation doesn't follow the appropriate galactic rotation curve, no n-body simulations do so without actually specifically implementing it in the code.

It might be useful to fork the code and create a modified simulator specifically for galaxy simulations.

 

[pbuk]

Sorry, how do you mean? The acceleration vector is re-calculated every step according to the netforce vector, which is reset and then cumulates gravitational forces at every step.

Yes, this is commonly known as Euler's method.

An example of a symplectic method and why it is more suitable: https://en.wikipedia.org/wiki/Verlet_integration

Adding a smoothing value is pretty standard for n-body solvers.

Yes (often called softening rather than smoothing), however not necessarily the method you are using (which is effectively Plummer softening) which has two critical disadvantages: firstly that it affects the calculation whether it is needed or not, and secondly (see below)...

Yes, it can do this. And, as you know, the smaller the timestep, the better.

No, your method adds the softening constant at every step so the smaller the timestep the LESS accurate it gets (Edit: but you won't normally see this in a 2-body stability test because there is no close approach and so $r^2 \gg \varepsilon ^ 2$ and the softening constant is lost in roundoff error).

 

[megapiano]

No, your method adds the softening constant at every step so the smaller the timestep the LESS accurate it gets

Hmm, I hadn't considered that. The softening constant is adjusted depending on the simulation, in attempt to make the error negligible, but I hadn't considered that the errors add up every step.

Here is what I have now:

// euler
netforce = sum(gravity)
acceleration = gravity / mass
velocity += acceleration * dt
position += velocity * dt

// verlet
position += velocity * dt + acceleration * dt * dt * 0.5
netforce = sum(gravity)
acceleration = netforce / mass
velocity += (accelerationOld + acceleration) * dt * 0.5
accelerationOld = acceleration

Is it necessary in the Verlet method to update the position before the other vectors?

 

[jedishrfu]

Thats the thing about simulations and why different integration methods are used because error is introduced at each step.

For montonically increasing or decreasing simulations you would choose a method that adds (or subtracts) small amounts of error like Euler.

https://en.wikipedia.org/wiki/Euler_method

For a periodic simulation like the three-body problem, you'd choose a method that introduces error based on the period sometime more sometimes less so that on average it evens out. If you were to use Euler, your planets would either collide or fly away.

When I took a computer sim course in grad school, the instructor noted that error relates to external energy added to or removed from the system under simulation to which some of us responded oh so we're living in a simulation now ala The Matrix movie. We all had a good laugh but then puzzled about it later.

 

[pbuk]

Is it necessary in the Verlet method to update the position before the other vectors?

Yes. For a derivation see https://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet or https://www.physics.udel.edu/~bnikolic/teaching/phys660/numerical_ode/node5.html.

Note that the algorithm you have quoted is specifically the velocity Verlet algorithm, not just Verlet (which is sometimes called Störmer-Verlet to avoid confusion).

 

[megapiano]

@jedishrfu
@pbuk
Thanks for the help.

Well, I was able to recover the part of the source-code I was looking for.

The force calculation is:
d = distance between particle1 and particle2 c = some constant forceGravity = (mass1 * mass2) * (1.0 - sin(min((d * d) / c, 1) * (pi / 2)));

Is anyone familiar with this method? I did not come up with it on my own, I remember adapting it from somewhere, but I can't remember exactly where it came from.

 

[jedishrfu]

It somewhat resembles Newtons equation: $f_r = G m_1 m_2 / r^2$

where G = 1 and the expression $(1 - sin ( \theta ) )$

with $\theta = min ( d^2 / c , 1 ) \pi/2$

so when $d^2 / c$ is between 0 and 1 it is used as the angle factor otherwise 1 is used

the min() insures the angle is between 0 and 90 degrees

and so the $(1 - sin ( \theta ) )$ varies between 0 and 1

 

[megapiano]

Yes, it also seems to make the gravitational field work in a finite distance.

I still cannot find the paper that implemented this specific algorithm.

I also think this method is only appropriate for n-body simulations of specifically galaxies, as it gives rise to spiral-arms and even the bar at the galactic center. I think I may use this method again in more simulations, but also implement a correct galactic rotation curve, which the simulation in the original post does not do.

 

[jedishrfu]

I found this reference detailing how various particle sim schemes work perhaps you'll find something in the links in the article.

https://www.cs.cmu.edu/afs/cs/academic/class/15850c-s96/www/nbody.html

and here's an article on deriving $1-sin(x) = sin^2(x) + cos^2(x) - 2 sin(x/2) cos(x/2)$ which might give you a clue of where the approximation came from.

https://www.quora.com/What-is-the-formula-for-1-sinx

 

[jedishrfu]

Rosetta code has an n-body task with many implementations:

https://rosettacode.org/wiki/N-body_problem

 

[pbuk]

when $d^2 / c$ is between 0 and 1 it is used as the angle factor otherwise 1 is used

... which limits the range of effect of the softening factor thus mitigating the adverse effect on conservation from Plummer softening mentioned above.

 

[pbuk]

Yes, it also seems to make the gravitational field work in a finite distance.

Not really, rather it acts so that the gravitational potential tends to 0 as $r \to 0$.

 

[pbuk]

here's an article on deriving $1-sin(x) = sin^2(x) + cos^2(x) - 2 sin(x/2) cos(x/2)$ which might give you a clue of where the approximation came from.

https://www.quora.com/What-is-the-formula-for-1-sinx

Really? I thought softening factors are usually derived from some assumption about mass distribution (with Plummer softening assuming a particle's mass is distributed between $0 \ \mathtt{ and } \ \infty$).