nbody/main.cpp at master · shreyExp/nbody · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
#include <iostream>
#include "math.h"
#include <armadillo>
#include <string>
#include "getini.h"
#include <sstream>
#include "simout.h"
#include "ProgressBar.hpp"
using namespace arma; //To  use anmadillo objects
using namespace std;
void minmax(double**, int, int);
inline Col<double> acc(Col<double> , Col<double> , double );
int main(int argc, char** argv){
	string f;
	if (argc>1){
		f = argv[1];
	}else{
		f = "parameters.json";
	}
	const char* fi = f.c_str();
	GetIni ini(fi);
	double *mass,  **x, **y;
	const int nbodies = ini.nOfBodies();
	const int iter = ini.getIteration();
	const double h = ini.getIterStepSize();
	Col<double> xa;
	ini.getXprojection(xa);// to be written
	Col<double> ya;
	ini.getYprojection(ya);// to be edited
	ini.getMass(mass); //to be written
	Col<double> *R;
	Col<double> *v;
	Col<double> CM;
	Col<double> a[nbodies];

	ini.getPositions(R);
	ini.getVelocities(v);
	bool cmflag = 0;
	if(ini.getCMFlag())
		cmflag = 1;
	x = new double*[iter];
	y = new double*[iter];
	for(int l = 0; l < iter; l++){
		x[l] = new double[nbodies];
		y[l] = new double[nbodies];
	}

	ProgressBar progressBar(iter, 70, '#', '-');

	cout<<"Iterating....."<<endl;
	double totalMass = 0;
	if(cmflag){
		for(int i = 0; i < nbodies; i++)
			totalMass += mass[i];
	}

	for(int k = 0; k < iter; k++){
		if(cmflag)
			CM<<0<<0<<0;
		for(int j = 0; j < nbodies; j++){//This can be parallelized
			a[j]<<0<<0<<0; //To stop accumulation : a should be new for new calculation
			for(int i = 0; i < nbodies; i++){
				if(i == j){
					continue;
				}
				a[j] += acc(R[j],R[i], mass[i]);
			}
			v[j] = v[j] + h*a[j];
			R[j] = R[j] + h*v[j];
			if(cmflag)
				CM = CM + mass[j]*R[j];
			x[k][j] = dot(R[j], xa);
			y[k][j] = dot(R[j], ya);
			//cout<<"x"<<"["<<k<<"]"<<"["<<j<<"]"<<"= "<<x[k][j]<<endl;
			//cout<<"y"<<"["<<k<<"]"<<"["<<j<<"]"<<"= "<<y[k][j]<<endl;
		}
		if(cmflag){
			CM = CM/totalMass;
			for(int j = 0; j < nbodies; j++){
				x[k][j] = x[k][j] - dot(CM, xa);
				y[k][j] = y[k][j] - dot(CM, ya);
			}
		}
		++progressBar;
		if(k % 10 == 0)
			progressBar.display();
	}
	//for(int i = 0; i < 10; i++)
	//	for(int j = 0; j < nbodies; j++)
	//		cout<<"x["<<i<<"]["<<j<<"] = "<<x[i][j]<<endl;

	cout<<"Iterating Complete.";
	string filename = ini.getVideoFileName();
	int duration = ini.getVideoDuration();
	//minmax(x, nbodies, iter);
	//Simout vid;
	//vid.x = x;
	//vid.y = y;
	//vid.nbodies = nbodies;
	//vid.iter = iter;
	//vid.duration = duration;
	//vid.filename = filename;
	//vid.videoGen();
	//cout<<" Number of bodies in the ,ain program: "<<nbodies<<endl;
	generateVideo(x, y, iter, nbodies, duration, filename);
}
inline Col<double> acc(Col<double> R1, Col<double> R2, double m){
	double G = 6.67408e-20;
	//double G = 6.67408e-7;
	Col<double> a = (G*m/pow((norm(R2-R1)),3))*(R2-R1);
	return a;
}
void minmax(double **x, int nbodies, int iter){
	double xmax = x[0][0];
	double xmin = x[0][0];
	for(int i = 0; i < iter; i++){
		for(int j = 0; j < nbodies; j++){
			if(xmax < x[i][j]){
				xmax = x[i][j];
			}
			if(xmin > x[i][j]){
				xmin = x[i][j];
			}
		}
	}
	cout<<"Max is"<<xmax<<endl;
	cout<<"Min is"<<xmin<<endl;

}