Skip to content

Instantly share code, notes, and snippets.

@earthnuker

earthnuker/nbody.cpp

Created October 25, 2016 20:58
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save earthnuker/8c14b0c2f365f1a8c0b4aa5c9e4bb29f to your computer and use it in GitHub Desktop.
Save earthnuker/8c14b0c2f365f1a8c0b4aa5c9e4bb29f to your computer and use it in GitHub Desktop.
Download ZIP
N-Body Simulation
Raw
nbody.cpp
#include <iostream>
#include <cmath>
#include <vector>
#include <random>
#include <chrono>
#include <algorithm>
#include <functional>
#include <fstream>
#include <omp.h>
#include <parallel/algo.h>
#define G 6.6738480
using namespace std;
class Body {
public:
long double x,y;
long double mass;
long double vx,vy;
Body(long double px,long double py) {
this->x=px;
this->y=py;
this->mass=10;
auto M=1.0/(0.1+sqrt((this->x)*(this->x)+(this->y)*(this->y)));
M*=0.01;
this->vx=-py*M*this->mass;
this->vy=px*M*this->mass;
}
Body(long double px,long double py,long double mass_) {
this->x=px;
this->y=py;
this->mass=mass_;
auto M=1.0/(0.1+sqrt((this->x)*(this->x)+(this->y)*(this->y)));
M*=0.01;
this->vx=-py*M*this->mass;
this->vy=px*M*this->mass;
}
long double dist(const Body& other) {
long double dx=(this->x)-(other.x);
long double dy=(this->y)-(other.y);
return sqrt(dx*dx+dy*dy);
}
bool operator==(const Body &other) {
return ((this->x==other.x)&&(this->y==other.y));
}
bool operator!=(const Body &other) {
return !(*this==other);
}
void print() {
cout<<"x="<<this->x<<" ";
cout<<"y="<<this->y<<" ";
cout<<"M="<<this->mass<<endl;
}
};
void progbar(string label,uintmax_t progress,uintmax_t total,uintmax_t width) {
cerr<<"\r"<<label<< " [";
uintmax_t pos = width * ((long double)progress/(long double)total);
for (uintmax_t i = 0; i < width; ++i) {
if (i < pos) cerr << "#";
else cerr << " ";
}
cerr << "] " << progress<<"/"<<total<<" (" << int(((long double)progress/(long double)total) * 100.0) << "%)";
cerr.flush();
}
int main(int argc,char *argv[]) {
omp_set_num_threads(omp_get_max_threads());
cout.setf(ios::fixed, ios::floatfield);
cout.setf(ios::showpoint);
auto dt=0.0001;
vector<Body> bodies,bodies_modified;
uintmax_t nb;
uintmax_t itcnt;
uintmax_t substeps;
if (argc<2) {
nb=10;
} else {
nb=atoi(argv[1]);
}
if (argc<3) {
itcnt=512;
} else {
itcnt=atoi(argv[2]);
}
if (argc<4) {
substeps=10;
} else {
substeps=atoi(argv[3]);
}
dt=dt/substeps;
const auto S=1.0;
long double c_time=0;
auto seed=chrono::system_clock::now().time_since_epoch().count();
mt19937_64 generator(seed);
uniform_real_distribution<long double> dist(-10,10);
auto rng=bind(dist,generator);
auto mass=rng()*10;
for (uintmax_t n=0;n<nb;n++) {
while (abs(mass)<1) mass=rng()*10;
bodies.emplace_back(rng()*5,rng()*5,50);
}
cout<<"IT,ID,PX,PY,VX,VY,M"<<endl;
for (auto b=bodies.begin(); b<bodies.end(); b++) {
auto i=b-bodies.begin();
cout<<0<<","<<i<<","<<0<<","<<b->x<<","<<b->y<<","<<b->vx<<","<<b->vy<<","<<b->mass<<endl;
}
cerr<<"Simulating "<<bodies.size()<<" Bodies for "<<itcnt<<" Iterations with "<<substeps<<" Substeps, Timestep: "<<dt<<endl;
for (uintmax_t itc=1; itc<=itcnt; itc++) {
string status="";
status="Calculating Iteration "+to_string(itc)+" of "+to_string(itcnt);
bodies_modified=bodies;
uintmax_t sstp;
for (sstp=1; sstp<=substeps; sstp++) {
progbar(status,sstp-1,substeps,25);
cerr<<" ["<<bodies.size()<<"]";
#pragma omp parallel for schedule(dynamic,256)
for (auto b1=bodies.begin(); b1<bodies.end(); b1++) {
for (auto b2=bodies.begin(); b2<bodies.end(); b2++) {
if (b1!=b2) {
auto d = b1->dist(*b2);
auto dx = b2->x-b1->x;
auto dy = b2->y-b1->y;
auto f = G*((b1->mass*b2->mass)/(d*d+S*S));
auto acc = 0.01*(f/abs(b1->mass))*dt;
auto i = b1-bodies.begin();
bodies_modified[i].vx+= acc*dx;
bodies_modified[i].vy+= acc*dy;
bodies_modified[i].x+= bodies_modified[i].vx*dt;
bodies_modified[i].y+= bodies_modified[i].vy*dt;
}
}
}
bodies=bodies_modified;
}
c_time+=dt;
for (auto b=bodies.begin(); b!=bodies.end(); b++) {
auto i=b-bodies.begin();
cout<<itc<<","<<i<<","<<c_time<<","<<b->x<<","<<b->y<<","<<b->vx<<","<<b->vy<<","<<b->mass<<endl;
}
progbar(status,sstp-1,substeps,25);
cerr<<" ["<<bodies.size()<<"] T:"<<c_time;
cerr<<endl;
/*
auto mass=rng()*10;
for (uintmax_t n=0;n<(nb);n++) {
while (abs(mass)<1) mass=rng()*10;
bodies.emplace_back(rng()*5,rng()*5,50);
}
*/
}
}
Raw
render.py
import numpy as np
import pylab as plt
import sys
data_format=[
('it','intc'),
('id','intc'),
('T','float64'),
('px','float64'),
('py','float64'),
('vx','float64'),
('vy','float64'),
('m','float64')
]
def gpot(x,y):
return lambda X,Y:((x-X)**2+(y-Y)**2)**0.9
def plot(local_data):
fnum=local_data['it'][0]
T=local_data['T'][0]
zoom=1
cx=np.average(local_data['px'])
cy=np.average(local_data['py'])
varx=np.max(local_data['px'])-np.min(local_data['px'])
vary=np.max(local_data['py'])-np.min(local_data['py'])
varm=np.max([varx,vary])/zoom
varm=150
sx,ex=(cx-varm,cx+varm)
sy,ey=(cy-varm,cy+varm)
#sx,ex=np.min(local_data['px']),np.max(local_data['px'])
#sy,ey=np.min(local_data['py']),np.max(local_data['py'])
plt.title("Time: {:f}".format(T))
plt.hist2d(local_data['px'],local_data['py'],bins=1000,range=[[sx,ex],[sy,ey]])
plt.savefig('Render_Out/h_{}.png'.format(fnum),dpi=200)
plt.cla()
if len(sys.argv)>1:
datafile=sys.argv[1]
else:
datafile=sys.stdin
head=next(datafile)
def getframe(handle):
pframe=None
for line in handle:
frame=int(line.split(",")[0])
if pframe is None:
pframe=frame
if frame!=pframe:
return
pframe=frame
yield line
if __name__=="__main__":
while 1:
try:
frame_data=np.loadtxt(getframe(datafile),dtype=data_format,skiprows=1,delimiter=",")
except StopIteration:
break
plot(frame_data)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment