# Simulating a closed particle system in Java

I have this program that simulates a closed particle system. Closed in this context means that the sum of all energies is constant. My primary concern is code itself, yet I would like to hear comments regarding physics as well. Here is my code:

Configuration.java

package net.coderodde.simulation;

public final class Configuration {

/**
* Defines the drawing scale. A distance of one unit length corresponds to
* the length of 100 pixels.
*/
public static final int PIXELS_PER_UNIT_LENGTH = 10;

/**
* The rejection force constant.
*/
public static final double FORCE_CONSTANT = 1000.0;
}

Particle.java

package net.coderodde.simulation;

import java.awt.Color;
import java.awt.Graphics;
import java.util.Objects;
import static net.coderodde.simulation.Configuration.FORCE_CONSTANT;
import static net.coderodde.simulation.Configuration.PIXELS_PER_UNIT_LENGTH;
import static net.coderodde.simulation.Utils.checkNonInfinite;
import static net.coderodde.simulation.Utils.checkNonNaN;
import static net.coderodde.simulation.Utils.checkNonNegative;

/**
* This class defines a particle in the simulation. The entire weight of a
* particle is considered to be fully focused in the center of this particle.
*
* @author Rodion "rodde" Efremov
* @version 1.6 (Sep 2, 2017)
*/
public final class Particle {

/**
* The mass of this particle.
*/
private final double mass;

/**
* The radius of the graphical representation of this particle.
*/

/**
* The color of the graphical representation of this particle.
*/
private final Color color;

/**
* The current x-coordinate of this particle.
*/
private double x;

/**
* The current y-coordinate of this particle.
*/
private double y;

/**
* The current velocity to the right. May be negative when the particle
* moves to the left.
*/
private double velocityX;

/**
* The current velocity downwards. May be negative when the particle moves
* upwards.
*/
private double velocityY;

/**
* Constructs a new particle.
*
* @param mass the weight of the new particle.
* @param color  the color of the new particle.
*/
public Particle(double mass, int radius, Color color) {
this.mass = checkMass(mass);
this.color = Objects.requireNonNull(color,
"The particle color is null.");
}

/**
* Copy-constructs a new particle.
*
* @param other the other particle to copy.
*/
public Particle(Particle other) {
this.mass      = other.mass;
this.color     = other.color;
this.x         = other.x;
this.y         = other.y;
this.velocityX = other.velocityX;
this.velocityY = other.velocityY;
}

public Vector getVelocityVector() {
return new Vector(velocityX, velocityY);
}

public double getX() {
return x;
}

public double getY() {
return y;
}

public double getVelocityX() {
return velocityX;
}

public double getVelocityY() {
return velocityY;
}

public double getMass() {
return mass;
}

public void setX(double x) {
this.x = checkX(x);
}

public void setY(double y) {
this.y = checkY(y);
}

public void setVelocityX(double velocityX) {
this.velocityX = checkVelocityX(velocityX);
}

public void setVelocityY(double velocityY) {
this.velocityY = checkVelocityY(velocityY);
}

/**
* Returns the current speed of this particle.
*
* @return the current speed.
*/
public double getSpeed() {
double vxSquared = velocityX * velocityX;
double vySquared = velocityY * velocityY;
return Math.sqrt(vxSquared + vySquared);
}

/**
* Returns the distance between this particle and {@code other}.
*
* @param other the other particle.
* @return the distance between two particles.
*/
public double getDistance(Particle other) {
double dx = x - other.x;
double dy = y - other.y;
return Math.sqrt(dx * dx + dy * dy);
}

/**
* Computes the kinetic energy of this particle.
*
* @return the kinetic energy.
*/
public double getKineticEnergy() {
double speed = getSpeed();
return 0.5 * mass * speed * speed;
}

/**
* Computes the rejection force between this and {@code other} particles.
*
* @param other the other particle.
* @return the rejection force.
*/
public double getRejectionForce(Particle other) {
double distance = getDistance(other);
return FORCE_CONSTANT * mass * other.getMass() / (distance * distance);
}

/**
* Computes the potential energy between this and {@code other} particle.
*
* @param other the other particle.
* @return potential energy.
*/
public double getPotentialEnergy(Particle other) {
return FORCE_CONSTANT * mass * other.getMass() / getDistance(other);
}

/**
* Draws this particle on a canvas.
*
* @param g the graphics context.
*/
public void draw(Graphics g) {
int effectiveX = (int)(x * PIXELS_PER_UNIT_LENGTH);
int effectiveY = (int)(y * PIXELS_PER_UNIT_LENGTH);

g.setColor(color);
}

@Override
public String toString() {
return "[x=" + x + ", y=" + y + ", velocityX=" + velocityX +
", velocityY=" + velocityY + "]";
}

private double checkMass(double mass) {
checkNonNaN(mass, "The particle mass is NaN.");
checkNonNegative(mass, "The particle mass is non-positive.");
checkNonInfinite(mass, "The particle mass is infinite.");
return mass;
}

throw new IllegalArgumentException(
}

}

private double checkCoordinate(double coordinate,
String errorMessageNaN,
String errorMessageInfinite) {
checkNonNaN(coordinate, errorMessageNaN);
checkNonInfinite(coordinate, errorMessageInfinite);
return coordinate;
}

private double checkX(double x) {
checkCoordinate(x,
"The x-coordinate is NaN.",
"The x-coordinate is infinite.");
return x;
}

private double checkY(double y) {
checkCoordinate(y,
"The y-coordinate is NaN.",
"The y-coordinate is infinite.");
return y;
}

private double checkVelocityX(double velocityX) {
checkCoordinate(velocityX,
"The x-velocity is NaN.",
"The x-velocity is infinite.");
return velocityX;
}

private double checkVelocityY(double velocityY) {
checkCoordinate(velocityY,
"The y-velocity is NaN.",
"The y-velocity is infinite.");
return velocityY;
}
}

SimulationApp.java

package net.coderodde.simulation;

import java.awt.Color;
import java.awt.Dimension;
import java.awt.Toolkit;
import java.util.ArrayList;
import java.util.List;
import java.util.Random;
import static net.coderodde.simulation.Configuration.PIXELS_PER_UNIT_LENGTH;

/**
* This class implements the entire simulation program.
*
* @author Rodion "rodde" Efremov
* @version 1.6 (Sep 2, 2017)
*/
public final class SimulationApp {

/**
* The minimum particle mass.
*/
private static final double MINIMUM_MASS = 15.0;

/**
* The maximum particle mass.
*/
private static final double MAXIMUM_MASS = 30.0;

/**
* Reserve the number of pixels for the title bar.
*/
private static final int TITLE_BAR_RESERVED_HEIGHT = 50;

/**
* The default number of particles in the simulation.
*/
private static final int DEFAULT_PARTICLES = 6;

/**
* The time step.
*/
private static final double TIME_STEP = 0.01;

/**
* The number of milliseconds spent between two consecutive time quants.
*/
private static final int SLEEP_TIME = 20;

/**
* Used for randomly generating the color components.
*/
private static final int COLOR_CHANNEL_MAX = 256;

/**
* The maximum initial velocity horizontally and/or vertically.
*/
private static final double MAX_INITIAL_VELOCITY = 40.0;

/**
* Defines the entry point of the program.
*
* @param args the command line arguments.
*/
public static void main(String[] args) {
Dimension screenDimension = Toolkit.getDefaultToolkit().getScreenSize();
screenDimension.height -= TITLE_BAR_RESERVED_HEIGHT;

double worldWidth = (1.0 * screenDimension.width)
/ PIXELS_PER_UNIT_LENGTH;

double worldHeight = (1.0 * screenDimension.height)
/ PIXELS_PER_UNIT_LENGTH;

long seed = System.currentTimeMillis();
Random random = new Random(seed);

System.out.println("Seed = " + seed);

List<Particle> particles = getParticles(DEFAULT_PARTICLES,
worldWidth,
worldHeight,
random);

SimulationCanvas simulationCanvas = new SimulationCanvas();

Simulator simulator = new Simulator(particles,
simulationCanvas,
worldWidth,
worldHeight,
TIME_STEP,
SLEEP_TIME);

simulationCanvas.setSimulator(simulator);

SimulationFrame simulationFrame =
new SimulationFrame(simulationCanvas,
screenDimension.width,
screenDimension.height);

SimulationFrameKeyListener keyListener =
new SimulationFrameKeyListener(simulator);

simulator.run();
}

private static List<Particle> getParticles(int particles,
double worldWidth,
double worldHeight,
Random random) {
List<Particle> particleList = new ArrayList<>(particles);

for (int i = 0; i < particles; ++i) {
worldWidth,
worldHeight));
}

return particleList;
}

private static Particle createRandomParticle(Random random,
double worldWidth,
double worldHeight) {
double mass = MINIMUM_MASS +
(MAXIMUM_MASS - MINIMUM_MASS) * random.nextDouble();
Color color = new Color(random.nextInt(COLOR_CHANNEL_MAX),
random.nextInt(COLOR_CHANNEL_MAX),
random.nextInt(COLOR_CHANNEL_MAX));

Particle particle = new Particle(mass, radius, color);

particle.setX(worldWidth * random.nextDouble());
particle.setY(worldHeight * random.nextDouble());

particle.setVelocityX(MAX_INITIAL_VELOCITY * random.nextDouble());
particle.setVelocityY(MAX_INITIAL_VELOCITY * random.nextDouble());

return particle;
}
}

SimulationCanvas.java

package net.coderodde.simulation;

import java.awt.Canvas;
import java.awt.Color;
import java.awt.Graphics;
import java.util.List;
import java.util.Objects;

/**
* This class implements a simple canvas for drawing the simulated particle
* system.
*
* @author Rodion "rodde" Efremov
* @version 1.6 (Sep 2, 2017)
*/
public final class SimulationCanvas extends Canvas {

/**
* The list of particles.
*/
private List<Particle> particles;

/**
* The simulation engine.
*/
private Simulator simulator;

@Override
public void paint(Graphics g) {
update(g);
}

@Override
public void update(Graphics g) {
double totalEnergy = simulator.computeTotalEnergy();
String totalEnergyString = "Total energy: " + totalEnergy;

g.setColor(getBackground());
g.clearRect(0, 0, getWidth(), getHeight());

for (Particle particle : particles) {
particle.draw(g);
}

g.setColor(Color.WHITE);
g.drawChars(totalEnergyString.toCharArray(),
0,
totalEnergyString.length(),
0,
20);
}

void setParticles(List<Particle> particles) {
this.particles = Objects.requireNonNull(
particles,
"The particle list is null.");
}

void setSimulator(Simulator simulator) {
this.simulator = simulator;
}
}

SimulationFrame.java

package net.coderodde.simulation;

import java.awt.Color;
import java.util.Objects;
import javax.swing.JFrame;

public final class SimulationFrame extends JFrame {

private static final String FRAME_TITLE = "Closed system simulation";

public SimulationFrame(SimulationCanvas simulationCanvas,
int width,
int height) {
super(FRAME_TITLE);
Objects.requireNonNull(simulationCanvas, "The input canvas is null.");
setSize(width, height);
simulationCanvas.setSize(width, height);
setLocation(0, 0);
setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
simulationCanvas.setBackground(Color.BLACK);
setResizable(false);
setVisible(true);
}
}

SimulationFrameKeyListener.java

package net.coderodde.simulation;

import java.awt.event.KeyEvent;
import java.awt.event.KeyListener;
import java.util.Objects;

public final class SimulationFrameKeyListener implements KeyListener {

private final Simulator simulator;

public SimulationFrameKeyListener(Simulator simulator) {
this.simulator = Objects.requireNonNull(simulator,
"The simulator is null.");
}

@Override
public void keyTyped(KeyEvent e) {
simulator.togglePause();
}

@Override
public void keyPressed(KeyEvent e) {

}

@Override
public void keyReleased(KeyEvent e) {

}
}

Simulator.java

package net.coderodde.simulation;

import java.util.List;
import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
import java.util.Map;
import java.util.Objects;
import static net.coderodde.simulation.Utils.checkNonInfinite;
import static net.coderodde.simulation.Utils.checkNonNaN;
import static net.coderodde.simulation.Utils.checkNonNegative;

/**
* This class implements the actual simulator.
*
* @author Rodion "rodde" Efremov
* @version 1.6 (Sep 2, 2017)
*/
public final class Simulator {

/**
* The list of particles.
*/
private final List<Particle> particles = new ArrayList<>();

/**
* Holds the canvas for drawing the system.
*/
private final SimulationCanvas simulationCanvas;

/**
* The time quant.
*/
private final double timeStep;

/**
* The total energy of the simulated system.
*/
private final double totalEnergy;

/**
* The width of the system.
*/
private final double worldWidth;

/**
* The height of the system.
*/
private final double worldHeight;

/**
* Number of milliseconds between two time quants.
*/
private final int sleepTime;

/**
* The exit flag
*/
private volatile boolean exit = false;

/**
* The pause flag.
*/
private volatile boolean pause = true;

/**
* Used for mapping the particles to their respective force vectors.
*/
private final Map<Particle, Vector> particleToForceVectorMap =
new HashMap<>();

public Simulator(List<Particle> particles,
SimulationCanvas simulationCanvas,
double worldWidth,
double worldHeight,
double timeStep,
int sleepTime) {
Objects.requireNonNull(particles, "The particle list is null.");

this.simulationCanvas =
Objects.requireNonNull(
simulationCanvas,
"The simulation canvas is null.");

checkNotEmpty(particles);
copy(particles);
checkParticlesDoNotOverlap();
this.worldWidth = checkWorldWidth(worldWidth);
this.worldHeight = checkWorldHeight(worldHeight);
this.timeStep = checkTimeStep(timeStep);
this.sleepTime = checkSleepTime(sleepTime);
totalEnergy = computeTotalEnergy();
simulationCanvas.setParticles(this.particles);
}

public void togglePause() {
pause = !pause;
}

public void run() {
while (!exit) {
if (!pause) {
performStep();
simulationCanvas.repaint();
}

sleep(sleepTime);
}
}

List<Particle> getParticles() {
return Collections.<Particle>unmodifiableList(particles);
}

/**
* Checks that the particle list is not empty.
*
* @param particles the particles list.
*/
private void checkNotEmpty(List<Particle> particles) {
if (particles.isEmpty()) {
throw new IllegalArgumentException("No particles given.");
}
}

/**
* Makes internal copies of all the particles so that client programmer
* cannot interfere.
*
* @param particles the particle list.
*/
private void copy(List<Particle> particles) {
for (Particle particle : particles) {
}
}

/**
* Performs one simulation step.
*/
private void performStep() {
// Compute the force vectors of all partices:
computeForceVectors();
updateParticleVelocities();
moveParticles();
resolveWorldBorderCollisions();
normalizeVelocityVectors();
particleToForceVectorMap.clear();
}

/**
* Computes all the repelling force vectors for each particle.
*/
private void computeForceVectors() {
for (Particle particle : particles) {
Vector vector = new Vector();

for (Particle other : particles) {
if (particle == other) {
// Do not compute the force from and to itself.
continue;
}

Vector aux = computeForceVector(particle, other);
vector = vector.plus(aux);
}

particleToForceVectorMap.put(particle, vector);
}
}

/**
* Computes a repelling force vector from {@code other} to {@code target}.
*
* @param target the target particle.
* @param other  the particle exerting repelling force towards
*               {@code target}.
* @return the force vector.
*/
private Vector computeForceVector(Particle target, Particle other) {
double vectorLength = target.getRejectionForce(other);
double dx = target.getX() - other.getX();
double dy = target.getY() - other.getY();
double angle = Math.atan2(dy, dx);
double xComponent = vectorLength * Math.cos(angle);
double yComponent = vectorLength * Math.sin(angle);
return new Vector(xComponent, yComponent);
}

/**
* Updates the velocities of each particle.
*/
private void updateParticleVelocities() {
for (Map.Entry<Particle, Vector> e
: particleToForceVectorMap.entrySet()) {
Particle particle = e.getKey();
Vector vector = e.getValue();
// Make the force 'vector' a acceleration vector:
vector = vector.multiply(1.0 / particle.getMass());

// Update the velocity components:
particle.setVelocityX(
particle.getVelocityX() + vector.getX() * timeStep);

particle.setVelocityY(
particle.getVelocityY() + vector.getY() * timeStep);
}
}

/**
* Moves all the particles.
*/
private void moveParticles() {
for (Particle particle : particles) {
particle.setX(particle.getX() + particle.getVelocityX() * timeStep);
particle.setY(particle.getY() + particle.getVelocityY() * timeStep);
}
}

/**
* Resolves all the border collisions.
*/
private void resolveWorldBorderCollisions() {
for (Particle particle : particles) {
if (particle.getY() <= 0.0 || particle.getY() >= worldHeight) {
particle.setVelocityY(-particle.getVelocityY());
}

if (particle.getX() <= 0.0 || particle.getX() >= worldWidth) {
particle.setVelocityX(-particle.getVelocityX());
}
}
}

/**
* Normalizes the current velocity vectors such that the total energy of the
* system remains constant.
*/
private void normalizeVelocityVectors() {
double totalEnergyDelta = computeTotalEnergyDelta();
double factor = getNormalizationConstant(totalEnergyDelta);

for (Particle particle : particles) {
particle.setVelocityX(factor * particle.getVelocityX());
particle.setVelocityY(factor * particle.getVelocityY());
}
}

/**
* Computes the difference between initial total energy and current total
* energy.
*
* @return the total energy difference.
*/
private double computeTotalEnergyDelta() {
double currentTotalEnergy = computeTotalEnergy();
double totalEnergyDelta = totalEnergy - currentTotalEnergy;
}

/**
* Computes such a velocity normalization constant, that the total energy of
* the system remains constant.
*
* @param totalEnergyDelta the difference of initial and current total
*                         energies.
* @return the velocity normalization constant.
*/
private double getNormalizationConstant(double totalEnergyDelta) {
double aux = totalEnergyDelta / computeTotalKineticEnergy() + 1;

if (aux < 0.0) {
return 1.0;
}

return Math.sqrt(aux);
}

/**
* Computes the sum of kinetic energies of all the particles.
*
* @return the sum of kinetic energies.
*/
private double computeTotalKineticEnergy() {
double kineticEnergy = 0.0;

for (Particle particle : particles) {
kineticEnergy += particle.getKineticEnergy();
}

return kineticEnergy;
}

/**
* Computes the current total energy.
*
* @return the current total energy.
*/
public double computeTotalEnergy() {
double totalEnergy = 0.0;

for (Particle particle : particles) {
totalEnergy += particle.getKineticEnergy();
}

for (int i = 0; i < particles.size(); ++i) {
Particle particle1 = particles.get(i);

for (int j = i + 1; j < particles.size(); ++j) {
Particle particle2 = particles.get(j);
totalEnergy += particle1.getPotentialEnergy(particle2);
}
}

}

/**
* Checks that there is no two different particles on the same spot.
*/
private void checkParticlesDoNotOverlap() {
for (int i = 0; i < particles.size(); ++i) {
Particle particle1 = particles.get(i);

for (int j = i + 1; j < particles.size(); ++j) {
Particle particle2 = particles.get(j);

if (particle1.getX() == particle2.getX()
&& particle1.getY() == particle2.getY()) {
throw new IllegalStateException(
"Two particles occupy the same spot.");
}
}
}
}

private double checkTimeStep(double timeStep) {
checkNonNaN(timeStep, "The time step is NaN.");
checkNonNegative(timeStep,
"The time step is non-positive: " + timeStep + ".");
checkNonInfinite(timeStep, "The time step is infinite.");
return timeStep;
}

private double checkWorldWidth(double worldWidth) {
return checkWorldDimension(
worldWidth,
"The world width is NaN.",
"The world width is non-positive: " + worldWidth,
"The world width is infinite.");
}

private double checkWorldHeight(double worldHeight) {
return checkWorldDimension(
worldHeight,
"The world height is NaN.",
"The world height is non-positive: " + worldHeight,
"The world height is infinite.");
}

private double checkWorldDimension(double dimension,
String errorMessageNaN,
String errorMessageNonPositive,
String errorMessageInfinite) {
checkNonNaN(dimension, errorMessageNaN);
checkNonNegative(dimension, errorMessageNonPositive);
checkNonInfinite(dimension, errorMessageInfinite);
return dimension;
}

private int checkSleepTime(int sleepTime) {
if (sleepTime < 1) {
throw new IllegalArgumentException(
"The sleep time is non-positive: " + sleepTime + ".");
}

return sleepTime;
}

private static void sleep(int milliseconds) {
try {
} catch (InterruptedException ex) {

}
}
}

Utils.java

package net.coderodde.simulation;

public final class Utils {

private Utils() {}

static void checkNonNaN(double value, String errorMessage) {
if (Double.isNaN(value)) {
throw new IllegalArgumentException(errorMessage);
}
}

static void checkNonNegative(double value, String errorMessage) {
if (value < 0.0) {
throw new IllegalArgumentException(errorMessage);
}
}

static void checkNonInfinite(double value, String errorMessage) {
if (Double.isInfinite(value)) {
throw new IllegalArgumentException(errorMessage);
}
}
}

Vector.java

package net.coderodde.simulation;

import static net.coderodde.simulation.Utils.checkNonInfinite;
import static net.coderodde.simulation.Utils.checkNonNaN;

/**
* This class implements a two-dimensional vector.
*
* @author Rodion "rodde" Efremov
* @version 1.6 (Sep 2, 2017)
*/
public final class Vector {

/**
* The x-component of this vector.
*/
private final double x;

/**
* The y-component of this vector.
*/
private final double y;

public Vector(double x, double y) {
this.x = checkX(x);
this.y = checkY(y);
}

public Vector() {
this(0.0, 0.0);
}

public double getX() {
return x;
}

public double getY() {
return y;
}

public Vector plus(Vector other) {
return new Vector(x + other.x, y + other.y);
}

public Vector multiply(double factor) {
return new Vector(x * factor, y * factor);
}

public double dotProduct(Vector other) {
return x * other.x + y * other.y;
}

@Override
public String toString() {
return "(x=" + x + ", y=" + y + ")";
}

private double checkX(double x) {
return check(x,
"The x-component is NaN.",
"The x-component is infinite.");
}

private double checkY(double y) {
return check(y,
"The y-component is NaN.",
"The y-component is infinite.");
}

private double check(double value,
String errorMessageNaN,
String errorMessageInfinite) {
checkNonNaN(value, errorMessageNaN);
checkNonInfinite(value, errorMessageInfinite);
return value;
}
}

Pressing any key toggles a pause on or off. If it does not work, first click the window so that it gets the focus.

Windows If you suffer from the flickering, check this repo.

I'll try to comment more on the physics side of things, since others have already done a good job of reviewing your code structure.

It seems that your overall physics looks pretty much exactly like an n-body simulation of Newtonian gravity, except that in your simulation the particles repel rather than attract each other. (Equivalently, you could say that your code is simulating Newtonian gravity, but with a negative gravitational constant.)

The reason I mention this is that there's a huge body of literature, code examples and even dedicated software libraries for efficient n-body Newtonian gravity simulation. While it's perfectly OK to "reinvent the wheel" and write your own as a learning exercise (or simply because you don't need anything fancy), if you ever get stuck or wish to explore new techniques for improving your simulation, it's useful to know that this huge body of existing research exists.

Anyway, let's move on to some specific criticism:

There are several places in your physics code where you first multiply something with a constant and later divide it by the same constant, or take the square root of something only to later square it again. Not only do such redundant calculations needlessly slow down your simulation, but they also reduce its accuracy. In particular:

• Instead of calling Math.sqrt() in Particle.getDistance() and then squaring the result in Particle.getRejectionForce(), you should give your Particle class a getDistanceSquared(Particle other) method that directly returns the squared distance.

• In your simulator class, you first compute the total force acting on each particle by summing up the repulsive forces between each pair of particles (which are proportional to their masses) and then divide this total force by the particle's mass to get the acceleration experienced by the particle.

It would be more efficient to rename your Particle.getRejectionForce() method into e.g. getRejectionAcceleration(), remove the multiplication by mass from it, and also remove the vector multiplication by 1.0 / particle.getMass() from Simulator.updateParticleVelocities(). (Of course, you should also replace any references to "force" with "acceleration" in your simulator code.)

In particular, this change would let your code properly handle zero-mass "test particles" that are repelled by the other particles but do not themselves exert any measurable repulsion forces on them. (If you wanted to have lots of such massless particles, you'd probably want to also add some extra optimizations for them so that you don't needlessly waste time computing the zero forces they exert on other particles. This would also ensure that you don't end up trying to calculate 0/0 = NaN if two test massless particles ever get pushed to the exact same position.)

Speaking of inefficient math, your Simulator.computeForceVector() method (which, BTW, really seems like it belongs in the Particle class) needlessly uses trig functions for the simple task of normalizing a vector. It would probably be more efficient to rewrite it something like this:

private Vector computeForceVector(Particle target, Particle other) {
double dx = target.getX() - other.getX();
double dy = target.getY() - other.getY();
Vector direction = new Vector(dx, dy).normalize();
return direction.multiply(target.getRejectionForce(other));
}

where the Vector.normalize() method would look something like this:

public Vector normalize() {
double factor = 1.0 / Math.sqrt(x * x + y * y);
return new Vector(x * factor, y * factor);
}

Alternatively, if you wished to minimize the number of temporary Vector objects constructed, you could divide dx and dy by the distance between the particles (and multiply them by the force / acceleration) before constructing the vector, e.g. like this:

private Vector computeForceVector(Particle target, Particle other) {
double dx = target.getX() - other.getX();
double dy = target.getY() - other.getY();
double factor = target.getRejectionForce(other) / Math.sqrt(dx*dx + dy*dy);
return new Vector(dx * factor, dy * factor);
}

You could even just merge this method with Particle.getRejectionForce(), which would save you from having to redundantly calculate dx, dy and the sum of their squares twice.

In any case, you might want to benchmark these different implementations and see how their performance compares, since it's not always obvious due to things like JIT optimizations.

You could also combine the force (or acceleration) calculation step and the particle velocity update step in your simulation code, since the forces / accelerations experienced by the particles don't depend on their velocities. (Even if you did later want to introduce velocity-dependent forces like friction, these would presumably only depend on each particle's own velocity.) Doing so would let you get rid of the particleToForceVectorMap entirely.

(You do still need to keep the position update step separate, since the forces experienced by a particle do depend on the positions of other particles.)

On the other hand, saving the acceleration of each particle until the position update step would let you replace the simple Euler integration with (velocity) Verlet integration. Compared to Euler integration, the Verlet method is more accurate and less sensitive to changes in the time step size. It's also time-reversible, such that reversing the velocities of all particles makes them retrace their path backwards exactly (up to the limits of floating-point accuracy, anyway), which tends to be a desirable feature for physics simulations.

If you decide to save the accelerations vectors, however, it would IMO make more sense to add the x and y acceleration components directly as properties of the Particle class, rather than storing them in a separate map. You could then also move the total force / acceleration calculation into the Particle class, letting each particle compute the total repulsive force it experiences, given a list of all the (other) particles in the system. Since your particles already handle the pairwise force calculations anyway, this would seem natural to me.

Finally, note that your current simulation algorithm takes O(n²) time to simulate a system with n particles, since it needs to compute the forces between each pair of particles. If you wish to simulate a large number of particles, you should consider using something like the Barnes–Hut algorithm, which can simulate n particles in O(n log n) time by approximating the forces using a spatial subdivision scheme.

• Barnes-Hut is a nice to know. Accepted! Oct 4, 2017 at 13:45

I don't like

:P - no, really:

/**
* The mass of this particle.
*/
private final double mass;

That should be really obvious.

And there's also those ones:

/**
* The time quant.
*/
private final double timeStep;

What's a quant? And is it the same as step?

For me, excessive comments are a code smell.

### Particle : separation of conerns

The particle itself has state for its presentation: radius and color. It also draws itself to a Graphics object. That's two jobs the Particle does, better split it up.

### Particle.check*() methods

The check methods do return a value, you then assign those return values to the member variables. That's a pattern I haven't seen very often. Usually, you first do your checks, and then assign the values.

(small detail: you can static import the requiredNonNull method)

### Particles.getKineticEnergy()

In the simulator, you call a method computeTotalEnergy, here it's get. I like compute. Compute is a good prefix, since you tell the user, that something is getting computed and it's not a state of the object. But at the same time, it is implementation detail: The user shouldn't care if it's state or if it's computed.

### Simulation.main

You can cast screenDimension.width to double, instead of multiplying it by 1.0.

There's no need to create a Random with the system time, default constructor does something similiar. Except you want to rerun the simulation of the same seed.

You got a bidirectional dependency for the SimulationCanvas and the Simulator. One of the main object oriented principles is lose coupling, having circular dependencies leads to the opposite.

In general, the main method is way too long, a few parts can be extracted to separate method, for instance the creation of the List<Particle> particles or the creation of the Simulator

### Simulation

DEFAULT_PARTICLES: Better: DEFAULT_NUMBER_OF_PARTICLES. Further down in the main method, the getParticles call would make more sense. It can be misinterpreted as a predefined list of particles.

getParticles: -> createParticles, or better: createRandomParticles. The particles parameter should be named amountOfParticles

### SimulationFrameKeyListener

you can extend KeyAdapter, so you don't have to overwrite all the methods.

### Simulator

You pass a list of Particles using the constructor, those are the copied (to avoid "messing stuff up"). At the same time, you provide a public method, which returns those particles. The list itself may be unmodifiable, but the Particles itself still are.

computeForceVectors: Instead of using continue, you can negate the if condition and move the calculation into the if block. To be honest: I don't like the use continue.

I don't understand, why the vectors are stored with the particles in a Map. Shouldn't the vector be a state of a Particle?

sleep: Even though in most of the time, the information, that a Thread has been interrupted during a sleep is not needed, I wouldn't leave the catch block blank.

The Simulator class it self is also too long for my liking. The simulator's job is mainly to calculate stuff. So, I think the threading can be decoupled from the Simulator, the threading is a separate concern anyway, at least in my opinion. Then, maybe introduce a package protected SimulatorCheckUtils class. Without those two things, the class should be a bit tidied up and only focuses on the main job.

Other than that: The code is quite straightforward.

When initializing the velocities, you currently generate only velocities with $v_x \ge 0, v_y \ge 0$. It should also be possible to have initial velocities that go in the other directions.