francosauvisky/README.md

## README.md

      
    Raw
  

              README.md
            
          
    Arbitrary Particle Simulator

What is this?

A small pair of programs to simulate physical point-like particles and plot them. The number of dimensions, forces and particles is arbitrary, limited only by memory and processing requirements.
The sim.c is responsible for generating the data from the initial conditions and the myplot.m is responsible for normalizing the data in time (using splines) and plotting it in real-time.
How to use?

The C code in this project comes with a 4-particle system. Taking that as example, it shouldn't be hard to adapt it to a different system. New forces can be easily created and added to a system. Note that you can add internal and external forces (every force is calculated with every pair, including itself, that it can be used as an external force).
The modifiers are generalized forces which can alter directly the particle properties as speed, position, etc. That's useful for collisions and similar processes which are easily described analytically but aren't very easy to implement only with forces.
If a particular force or modifier requires some particle property (for example, charge), one can easily add it to the Particle struct and include it on the system.
The myplot.m is a Octave script that creates an animated real-time plot of the system, so it's important to keep the time scale on a visible range (not too small neither too large), although it can be easily adapted.
Running

To run the simulator, just put every files on the same folder and run octave -qf ./myplot.m there. It will compile the code, run it and plot it automatically. To print GIF you need Imagemagick and Gifsicle installed on your system.

  
## myplot.m
% Copyright (c) 2016, Franco Sauvisky
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions
% are met:
% 1. Redistributions of source code must retain the above copyright
%    notice, this list of conditions and the following disclaimer.
% 2. Redistributions in binary form must reproduce the above copyright
%    notice, this list of conditions and the following disclaimer in the
%    documentation and/or other materials provided with the distribution.
% 3. The name of the author may not be used to endorse or promote products
%    derived from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
% IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
% OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
% INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
% NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
% THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
% THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

disp("Calculating..."); fflush(stdout);

fflag = system("gcc sim.c -O3 -lm -o sim && ./sim > orbit1.dat");
if fflag != 0
	return;
end

disp("Done. Processing..."); fflush(stdout);

fps = 30;

fh = fopen("./orbit1.dat");
raw = fscanf(fh, "%f");
fclose(fh);

n = raw(1); % Number of particles
raw = raw(2:end);

data = reshape(raw, [2*n+1, length(raw)/(2*n+1)])';

tmax = max(data(:,1));
normaldata = [(0:1/fps:tmax)', spline(data(:,1), data(:,2:end),0:1/fps:tmax)'];
datalength = size(normaldata, 1);

baseaxis = [min(min(data(:,[2:2:end]))),max(max(data(:,[2:2:end]))),...
min(min(data(:,[3:2:end]))),max(max(data(:,[3:2:end])))];
mycmap = colormap(lines);

cgif = 0;
while cgif != 3
	cgif = menu("Done.", "Plot real-time", "Create GIF file", "Exit");

	if cgif == 1
		if isempty(findall(0,'Type','Figure'))
			figure()
		end

		clf;
		disp("Click on figure to plot..."); fflush(stdout);
		waitforbuttonpress();

		i = 1;
		tic;
		while i < datalength
			scatter(normaldata(i,[2:2:end]),normaldata(i,[3:2:end]),5,mycmap(1:n,:),"filled");
			grid on
			axis(baseaxis,"equal");
		title(["t = " num2str(normaldata(i,1),"%2.2f")]);
		pause(0.001);
		i = round(fps*toc);
		end

		hold on
		for i = 1:n
			plot(data(:,2*i),data(:,2*i+1),"color", mycmap(i,:), "linewidth", 1);
		end
		hold off
	end

	if cgif == 2
		disp("Printing frames..."); fflush(stdout);
		set(0, 'defaultfigurevisible', 'off');
		mkdir output;
		for t = 1:size(normaldata,1)
			scatter(normaldata(t,[2:2:end]),normaldata(t,[3:2:end]),5,mycmap(1:n,:),"filled");
			grid on
			axis(baseaxis,"equal");
			title(["t = " num2str(normaldata(t,1),"%2.2f")]);
			filename=sprintf('./output/frame%d.png',t);
			print("-r0", filename);
		end
		set(0, 'defaultfigurevisible', 'on');
		disp("Done. Converting and creating gif..."); fflush(stdout);
		system("for i in $(ls ./output/*.png); do convert $i $i.gif; rm $i; done");
		system(["gifsicle --loop --delay=" num2str(round(100/fps)) " $(ls -v1 ./output/*.gif) > output.gif"]);
		system("rm -r ./output");
	end
end

disp("Done!")

## sim.c
// Copyright (c) 2016, Franco Sauvisky
// All rights reserved.

// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions
// are met:
// 1. Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
// 2. Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
// 3. The name of the author may not be used to endorse or promote products
//    derived from this software without specific prior written permission.

// THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
// OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
// IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
// NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
// THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#include <stdio.h>
#include <string.h>
#include <stdbool.h>
#include <math.h>

#define dim 2
#define maxParticles 10
#define maxForces 5
#define maxModifiers 5

#define maxTimeStep 1e-2
#define minTimeStep 1e-6
#define dt0 1e-4

#include "sim_internal.c"

#define tmax 10

void grav_force(double *, struct Particle, struct Particle);
void eletr_force(double *, struct Particle, struct Particle);
void mag_ext_force(double *, struct Particle, struct Particle);
void eletr_ext_force(double *, struct Particle, struct Particle);
void drag_force(double *, struct Particle, struct Particle);
void square_wall(struct Particle *, struct Particle *);
void sphere_collision(struct Particle *, struct Particle *);

int
main(void)
{
	struct point_system psystem = {
		.Bodies = {
			{
				.pos = {0,0},
				.vel = {0,0},
				.mass = 10,
				.radius = 0.5,
				.charge = 10,
				.fixed = false
			},
			{
				.pos = {3.9,4.1},
				.vel = {-8,-2},
				.mass = 10,
				.radius = 0.1,
				.charge = -10,
				.fixed = false
			},
			{
				.pos = {3,-3},
				.vel = {0,4},
				.mass = 10,
				.radius = 0.3,
				.charge = -10,
				.fixed = false
			},
			{
				.pos = {-3,3},
				.vel = {0,-4},
				.mass = 20,
				.radius = 0.3,
				.charge = -0,
				.fixed = false
			},
		},
		.Forces = {
			grav_force,
			eletr_force,
			mag_ext_force,
			eletr_ext_force,
			drag_force,
		},
		.Modifiers = {
			square_wall,
			sphere_collision,
		},
	};

	config_system(&psystem);

	printf("%d ", psystem.particleN);
	print_system(psystem);

	for(; psystem.time < tmax;)
	{
		double dt = dynamic_dt(psystem);
		calc_system(&psystem, dt);
		print_system(psystem);
	}
}

#define grav_c 1
void
grav_force(double *force, struct Particle particle_0, struct Particle particle_1)
{
	if(particle_0.id == particle_1.id || particle_0.fixed) return;

	double r_01[dim];

	sum_vec(r_01, particle_1.pos, particle_0.pos, -1);

	double dist = r_abs(r_01);

	for(int i = 0; i < dim; i++)
	{
		force[i] += grav_c * particle_0.mass * particle_1.mass * r_01[i] / pow(dist, 3);
	}
}

#define eletr_c 2
void
eletr_force(double *force, struct Particle particle_0, struct Particle particle_1)
{
	if(particle_0.id == particle_1.id || particle_0.fixed) return;

	double r_01[dim];

	sum_vec(r_01, particle_1.pos, particle_0.pos, -1);

	double dist = r_abs(r_01);

	for(int i = 0; i < dim; i++)
	{
		force[i] -= eletr_c * particle_0.charge * particle_1.charge * r_01[i] / pow(dist, 3);
	}
}

#define mag_ext_c 1
void
mag_ext_force(double *force, struct Particle particle_0, struct Particle particle_1)
{
	if(particle_0.id != particle_1.id || particle_0.fixed) return;

	force[0] += mag_ext_c * particle_1.charge * particle_0.vel[1];
	force[1] -= mag_ext_c * particle_1.charge * particle_0.vel[0];
}

#define eletr_ext_c 20
void
eletr_ext_force(double *force, struct Particle particle_0, struct Particle particle_1)
{
	if(particle_0.id != particle_1.id || particle_0.fixed) return;

	force[0] += eletr_ext_c * particle_1.charge;
}

#define drag_c 2
void
drag_force(double *force, struct Particle particle_0, struct Particle particle_1)
{
	if(particle_0.id != particle_1.id || particle_0.fixed) return;

	for(int i = 0; i < dim; i++)
	{
		force[i] -= drag_c*particle_0.vel[i];
	}
}

#define topw 5
#define bottomw -5
#define rightw 5
#define leftw -5
void
square_wall(struct Particle *particle_0, struct Particle *particle_1)
{
	if((*particle_0).id != (*particle_1).id) return;

	if((*particle_0).pos[1] + (*particle_0).radius > topw && (*particle_0).vel[1] > 0)
		(*particle_0).vel[1] *= -1;

	if((*particle_0).pos[0] + (*particle_0).radius > rightw && (*particle_0).vel[0] > 0)
		(*particle_0).vel[0] *= -1;

	if((*particle_0).pos[1] - (*particle_0).radius < bottomw && (*particle_0).vel[1] < 0)
		(*particle_0).vel[1] *= -1;

	if((*particle_0).pos[0] - (*particle_0).radius < leftw && (*particle_0).vel[0] < 0)
		(*particle_0).vel[0] *= -1;
}

void
sphere_collision(struct Particle *particle_0, struct Particle *particle_1)
{
	if((*particle_0).id >= (*particle_1).id) return;

	double r_01[dim];

	sum_vec(r_01, (*particle_1).pos, (*particle_0).pos, -1);
	double dist = r_abs(r_01);

	if(dist < (*particle_0).radius + (*particle_1).radius)
	{
		r_norm(r_01, r_01);

		double proj0, proj1, vcm;

		proj0 = dot_product((*particle_0).vel, r_01);
		proj1 = dot_product((*particle_1).vel, r_01);

		if(proj0 < proj1) return;

		if(!(*particle_0).fixed && !(*particle_1).fixed)
		{
			vcm = (proj0 * (*particle_0).mass + proj1 * (*particle_1).mass)/
			((*particle_0).mass + (*particle_1).mass);
		}
		else
		{
			vcm = 0;
		}

		sum_vec((*particle_0).vel, (*particle_0).vel, r_01, 2*vcm - 2*proj0);
		sum_vec((*particle_1).vel, (*particle_1).vel, r_01, 2*vcm - 2*proj1);
	}
}

## sim_internal.c
// Copyright (c) 2016, Franco Sauvisky
// All rights reserved.

// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions
// are met:
// 1. Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
// 2. Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
// 3. The name of the author may not be used to endorse or promote products
//    derived from this software without specific prior written permission.

// THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
// OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
// IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
// NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
// THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

double zero_vec[dim] = {0};

struct Particle
{
	int id;
	double pos[dim];
	double vel[dim];
	double mass;
	double radius;
	double charge;
	bool fixed;
};

struct point_system
{
	int particleN;
	int forceN;
	int modifierN;
	double time;
	struct Particle Bodies[maxParticles];
	void (*Forces[maxForces])(double*, struct Particle, struct Particle);
	void (*Modifiers[maxModifiers])(struct Particle*, struct Particle*);
};

void
config_system(struct point_system *system)
{
	(*system).time = 0;

	(*system).particleN = 0;
	for (int i = 0; i < maxParticles; ++i)
	{
		if((*system).Bodies[i].mass != 0)
			++(*system).particleN;
	}

	(*system).forceN = 0;
	for (int i = 0; i < maxForces; ++i)
	{
		if((*system).Forces[i] != NULL)
			++(*system).forceN;
	}

	(*system).modifierN = 0;
	for (int i = 0; i < maxModifiers; ++i)
	{
		if((*system).Modifiers[i] != NULL)
			++(*system).modifierN;
	}

	for (int i = 0; i < (*system).particleN; i++)
	{
		(*system).Bodies[i].id = i+1;
	}
}

void
sum_vec(double *c, double *a, double *b, double multiplier_constant)
{
	for(int i = 0; i < dim; i++)
	{
		c[i] = a[i] + multiplier_constant*b[i];
	}
}

double
dot_product(double *a, double *b)
{
	double cnt = 0;
	for(int i = 0; i < dim; i++)
	{
		cnt += a[i]*b[i];
	}
	return cnt;
}

double
r_abs(double *radius)
{
	return sqrt(dot_product(radius, radius));
}

void
r_norm(double *b, double *a)
{
	double norm = r_abs(a);
	for(int i = 0; i < dim; i++)
	{
		b[i] = a[i]/norm;
	}
}

void
calc_system(struct point_system *mysystem, double dt)
{
	for(int i = 0; i < (*mysystem).particleN; i++)
	{
		double force[dim] = {0};

		for(int j = 0; j < (*mysystem).particleN; j++)
		{
			for(int k = 0; k < (*mysystem).forceN; k++)
			{
				(*mysystem).Forces[k](force, (*mysystem).Bodies[i], (*mysystem).Bodies[j]);
			}

			for(int k = 0; k < (*mysystem).modifierN; k++)
			{
				(*mysystem).Modifiers[k](&(*mysystem).Bodies[i], &(*mysystem).Bodies[j]);
			}
		}

		if(!(*mysystem).Bodies[i].fixed)
		{
			sum_vec((*mysystem).Bodies[i].vel, (*mysystem).Bodies[i].vel, force, dt/(*mysystem).Bodies[i].mass);
			sum_vec((*mysystem).Bodies[i].pos, (*mysystem).Bodies[i].pos, (*mysystem).Bodies[i].vel, dt);
		}
	}

	(*mysystem).time += dt;
}

void
print_system(struct point_system system)
{
	printf("%lf", system.time);

	for(int i = 0; i < system.particleN; i++)
	{
		for(int j = 0; j < dim; j++)
		{
			printf(" %lf", system.Bodies[i].pos[j]);
		}
	}
	printf("\n");
}

// dt = dt0 * min(dist/vel) where dt0 << 1, so that vel*dt << dist
double
dynamic_dt(struct point_system system)
{
	double dt = maxTimeStep;

	for (int i = 0; i < system.particleN; i++)
	{
		if(system.Bodies[i].fixed == false)
		{
			double vel = r_abs(system.Bodies[i].vel);
			for (int j = 0; j < system.particleN; j++)
			{
				if(i != j)
				{
					double rad[dim], dist;
					sum_vec(rad, system.Bodies[i].pos, system.Bodies[j].pos, -1);
					dist = r_abs(rad);
					if(dt0*dist/vel < dt)
					{
						dt = dt0*dist/vel;
					}
				}

			}
		}
	}
	return fmax(dt, minTimeStep);
}
	% Copyright (c) 2016, Franco Sauvisky
	% All rights reserved.

	% Redistribution and use in source and binary forms, with or without
	% modification, are permitted provided that the following conditions
	% are met:
	% 1. Redistributions of source code must retain the above copyright
	% notice, this list of conditions and the following disclaimer.
	% 2. Redistributions in binary form must reproduce the above copyright
	% notice, this list of conditions and the following disclaimer in the
	% documentation and/or other materials provided with the distribution.
	% 3. The name of the author may not be used to endorse or promote products
	% derived from this software without specific prior written permission.

	% THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
	% IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
	% OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
	% IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
	% INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
	% NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
	% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
	% THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
	% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
	% THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

	disp("Calculating..."); fflush(stdout);

	fflag = system("gcc sim.c -O3 -lm -o sim && ./sim > orbit1.dat");
	if fflag != 0
	return;
	end

	disp("Done. Processing..."); fflush(stdout);

	fps = 30;

	fh = fopen("./orbit1.dat");
	raw = fscanf(fh, "%f");
	fclose(fh);

	n = raw(1); % Number of particles
	raw = raw(2:end);

	data = reshape(raw, [2n+1, length(raw)/(2n+1)])';

	tmax = max(data(:,1));
	normaldata = [(0:1/fps:tmax)', spline(data(:,1), data(:,2:end),0:1/fps:tmax)'];
	datalength = size(normaldata, 1);

	baseaxis = [min(min(data(:,[2:2:end]))),max(max(data(:,[2:2:end]))),...
	min(min(data(:,[3:2:end]))),max(max(data(:,[3:2:end])))];
	mycmap = colormap(lines);

	cgif = 0;
	while cgif != 3
	cgif = menu("Done.", "Plot real-time", "Create GIF file", "Exit");

	if cgif == 1
	if isempty(findall(0,'Type','Figure'))
	figure()
	end

	clf;
	disp("Click on figure to plot..."); fflush(stdout);
	waitforbuttonpress();

	i = 1;
	tic;
	while i < datalength
	scatter(normaldata(i,[2:2:end]),normaldata(i,[3:2:end]),5,mycmap(1:n,:),"filled");
	grid on
	axis(baseaxis,"equal");
	title(["t = " num2str(normaldata(i,1),"%2.2f")]);
	pause(0.001);
	i = round(fps*toc);
	end

	hold on
	for i = 1:n
	plot(data(:,2i),data(:,2i+1),"color", mycmap(i,:), "linewidth", 1);
	end
	hold off
	end

	if cgif == 2
	disp("Printing frames..."); fflush(stdout);
	set(0, 'defaultfigurevisible', 'off');
	mkdir output;
	for t = 1:size(normaldata,1)
	scatter(normaldata(t,[2:2:end]),normaldata(t,[3:2:end]),5,mycmap(1:n,:),"filled");
	grid on
	axis(baseaxis,"equal");
	title(["t = " num2str(normaldata(t,1),"%2.2f")]);
	filename=sprintf('./output/frame%d.png',t);
	print("-r0", filename);
	end
	set(0, 'defaultfigurevisible', 'on');
	disp("Done. Converting and creating gif..."); fflush(stdout);
	system("for i in $(ls ./output/*.png); do convert $i $i.gif; rm $i; done");
	system(["gifsicle --loop --delay=" num2str(round(100/fps)) " $(ls -v1 ./output/*.gif) > output.gif"]);
	system("rm -r ./output");
	end
	end

	disp("Done!")
	// Copyright (c) 2016, Franco Sauvisky
	// All rights reserved.

	// Redistribution and use in source and binary forms, with or without
	// modification, are permitted provided that the following conditions
	// are met:
	// 1. Redistributions of source code must retain the above copyright
	// notice, this list of conditions and the following disclaimer.
	// 2. Redistributions in binary form must reproduce the above copyright
	// notice, this list of conditions and the following disclaimer in the
	// documentation and/or other materials provided with the distribution.
	// 3. The name of the author may not be used to endorse or promote products
	// derived from this software without specific prior written permission.

	// THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
	// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
	// OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
	// IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
	// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
	// NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
	// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
	// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
	// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
	// THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

	#include <stdio.h>
	#include <string.h>
	#include <stdbool.h>
	#include <math.h>

	#define dim 2
	#define maxParticles 10
	#define maxForces 5
	#define maxModifiers 5

	#define maxTimeStep 1e-2
	#define minTimeStep 1e-6
	#define dt0 1e-4

	#include "sim_internal.c"

	#define tmax 10

	void grav_force(double *, struct Particle, struct Particle);
	void eletr_force(double *, struct Particle, struct Particle);
	void mag_ext_force(double *, struct Particle, struct Particle);
	void eletr_ext_force(double *, struct Particle, struct Particle);
	void drag_force(double *, struct Particle, struct Particle);
	void square_wall(struct Particle , struct Particle );
	void sphere_collision(struct Particle , struct Particle );

	int
	main(void)
	{
	struct point_system psystem = {
	.Bodies = {
	{
	.pos = {0,0},
	.vel = {0,0},
	.mass = 10,
	.radius = 0.5,
	.charge = 10,
	.fixed = false
	},
	{
	.pos = {3.9,4.1},
	.vel = {-8,-2},
	.mass = 10,
	.radius = 0.1,
	.charge = -10,
	.fixed = false
	},
	{
	.pos = {3,-3},
	.vel = {0,4},
	.mass = 10,
	.radius = 0.3,
	.charge = -10,
	.fixed = false
	},
	{
	.pos = {-3,3},
	.vel = {0,-4},
	.mass = 20,
	.radius = 0.3,
	.charge = -0,
	.fixed = false
	},
	},
	.Forces = {
	grav_force,
	eletr_force,
	mag_ext_force,
	eletr_ext_force,
	drag_force,
	},
	.Modifiers = {
	square_wall,
	sphere_collision,
	},
	};

	config_system(&psystem);

	printf("%d ", psystem.particleN);
	print_system(psystem);

	for(; psystem.time < tmax;)
	{
	double dt = dynamic_dt(psystem);
	calc_system(&psystem, dt);
	print_system(psystem);
	}
	}

	#define grav_c 1
	void
	grav_force(double *force, struct Particle particle_0, struct Particle particle_1)
	{
	if(particle_0.id == particle_1.id \|\| particle_0.fixed) return;

	double r_01[dim];

	sum_vec(r_01, particle_1.pos, particle_0.pos, -1);

	double dist = r_abs(r_01);

	for(int i = 0; i < dim; i++)
	{
	force[i] += grav_c * particle_0.mass * particle_1.mass * r_01[i] / pow(dist, 3);
	}
	}

	#define eletr_c 2
	void
	eletr_force(double *force, struct Particle particle_0, struct Particle particle_1)
	{
	if(particle_0.id == particle_1.id \|\| particle_0.fixed) return;

	double r_01[dim];

	sum_vec(r_01, particle_1.pos, particle_0.pos, -1);

	double dist = r_abs(r_01);

	for(int i = 0; i < dim; i++)
	{
	force[i] -= eletr_c * particle_0.charge * particle_1.charge * r_01[i] / pow(dist, 3);
	}
	}

	#define mag_ext_c 1
	void
	mag_ext_force(double *force, struct Particle particle_0, struct Particle particle_1)
	{
	if(particle_0.id != particle_1.id \|\| particle_0.fixed) return;

	force[0] += mag_ext_c * particle_1.charge * particle_0.vel[1];
	force[1] -= mag_ext_c * particle_1.charge * particle_0.vel[0];
	}

	#define eletr_ext_c 20
	void
	eletr_ext_force(double *force, struct Particle particle_0, struct Particle particle_1)
	{
	if(particle_0.id != particle_1.id \|\| particle_0.fixed) return;

	force[0] += eletr_ext_c * particle_1.charge;
	}

	#define drag_c 2
	void
	drag_force(double *force, struct Particle particle_0, struct Particle particle_1)
	{
	if(particle_0.id != particle_1.id \|\| particle_0.fixed) return;

	for(int i = 0; i < dim; i++)
	{
	force[i] -= drag_c*particle_0.vel[i];
	}
	}

	#define topw 5
	#define bottomw -5
	#define rightw 5
	#define leftw -5
	void
	square_wall(struct Particle particle_0, struct Particle particle_1)
	{
	if((particle_0).id != (particle_1).id) return;

	if((particle_0).pos[1] + (particle_0).radius > topw && (*particle_0).vel[1] > 0)
	(particle_0).vel[1] = -1;

	if((particle_0).pos[0] + (particle_0).radius > rightw && (*particle_0).vel[0] > 0)
	(particle_0).vel[0] = -1;

	if((particle_0).pos[1] - (particle_0).radius < bottomw && (*particle_0).vel[1] < 0)
	(particle_0).vel[1] = -1;

	if((particle_0).pos[0] - (particle_0).radius < leftw && (*particle_0).vel[0] < 0)
	(particle_0).vel[0] = -1;
	}

	void
	sphere_collision(struct Particle particle_0, struct Particle particle_1)
	{
	if((particle_0).id >= (particle_1).id) return;

	double r_01[dim];

	sum_vec(r_01, (particle_1).pos, (particle_0).pos, -1);
	double dist = r_abs(r_01);

	if(dist < (particle_0).radius + (particle_1).radius)
	{
	r_norm(r_01, r_01);

	double proj0, proj1, vcm;

	proj0 = dot_product((*particle_0).vel, r_01);
	proj1 = dot_product((*particle_1).vel, r_01);

	if(proj0 < proj1) return;

	if(!(particle_0).fixed && !(particle_1).fixed)
	{
	vcm = (proj0 * (particle_0).mass + proj1 (*particle_1).mass)/
	((particle_0).mass + (particle_1).mass);
	}
	else
	{
	vcm = 0;
	}

	sum_vec((particle_0).vel, (particle_0).vel, r_01, 2vcm - 2proj0);
	sum_vec((particle_1).vel, (particle_1).vel, r_01, 2vcm - 2proj1);
	}
	}