Skip to content
Snippets Groups Projects
Commit fe61d58d authored by Johannes Lukas Bartunek's avatar Johannes Lukas Bartunek
Browse files

Delete main.c

parent 255b0a25
Branches
No related tags found
No related merge requests found
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//compile options gcc -O2 main.c -Wall -lm
typedef struct {
double x[3];
double v[3];
double a[3];
double mass;
} Planet;
Planet Planets[] = {
{
{0,0,0},
{0,0,0},
{0,0,0},
1.98847*pow(10,30)
},
{
{149597870700,0,0},
{0,30080,0},
{0,0,0},
5.9722*pow(10,24)
}
};
//TODO add softening factor
static inline double distance( const double x1[3], const double x2[3]) {
const double x = x1[0] - x2[0];
const double y = x1[1] - x2[1];
const double z = x1[2] - x2[2];
return sqrt(x*x + y*y + z*z);
}
double energy(int n, double G){
double e;
int i=0, j=0, k=0;
e = 0.0;
for(i=0; i<n; i++){
for(k=0; k<3; k++){
e += 0.5 * Planets[i].mass * Planets[i].v[k] * Planets[i].v[k];
}
for(j=i+1; j < n; j++) {
double dist = distance(Planets[i].x, Planets[j].x);
e -= (G * Planets[i].mass * Planets[j].mass) / dist;
}
}
return e;
}
int main(int argc, char * argv[]) {
const long steps = atoi(argv[1]);
/*if (argc == 1) {
const long steps = atoi(argv[1]);
} else {
const long steps = 10000;
}*/
//FILE *output = fopen("output.xyz", "w+");
FILE *output = fopen("test", "w+");
//defer fclose(output);
if (output == NULL) {
printf("HELP\n COULDNT OPEN FILE\nSIGSEVV"); //TODO stop here
} else {
printf("Output File opened");
}
//this could also be done via the preprocessor
const int noOfBodies = sizeof(Planets)/sizeof(Planet);
printf("size: %d", noOfBodies);
const double G = 6.6743015*pow(10,-11);
const double dt = 0.1;
printf("energy before: %.lf\n", energy(noOfBodies, G));
//double vx, vy, vz;
long run = 0;
//double acc_buff[3];
while (run < steps) {
for (int l = 0; l < noOfBodies; l++) {
Planets[l].v[0] += 0.5 * Planets[l].a[0]* dt;
Planets[l].v[1] += 0.5 * Planets[l].a[1]* dt;
Planets[l].v[2] += 0.5 * Planets[l].a[2]* dt;
if (run %(long int) pow(10,7)== 0) {
//fprintf(output, "%d %lf %lf %lf\n", index, Planets[index].x[0], Planets[index].x[1], Planets[index].x[2]);
fprintf(output, "%lf,%lf,%lf,", Planets[l].x[0], Planets[l].x[1], Planets[l].x[2]);
}
}
for (int k = 0; k < noOfBodies; k++) {
Planets[k].x[0] += dt * Planets[k].v[0];
Planets[k].x[1] += dt * Planets[k].v[1];
Planets[k].x[2] += dt * Planets[k].v[2];
Planets[k].a[0] = 0;
Planets[k].a[1] = 0;
Planets[k].a[2] = 0;
//printf("%d @ %d; ---x: %lf, y: %lf, z: %lf\n\n",run, k, Planets[k].x[0], Planets[k].x[1], Planets[k].x[2]);
}
for (int i = 0; i < noOfBodies; i++) {
for (int j = i+1; j < noOfBodies; j++) {
double dist = distance(Planets[i].x, Planets[j].x);
double a = G/pow(dist, 3);
//this might be the wrong way.......
Planets[i].a[0] += a * Planets[j].mass * ( Planets[j].x[0] - Planets[i].x[0] );
Planets[i].a[1] += a * Planets[j].mass * ( Planets[j].x[1] - Planets[i].x[1] );
Planets[i].a[2] += a * Planets[j].mass * ( Planets[j].x[2] - Planets[i].x[2] );
Planets[j].a[0] += a * Planets[i].mass * ( Planets[i].x[0] - Planets[j].x[0] );
Planets[j].a[1] += a * Planets[i].mass * ( Planets[i].x[1] - Planets[j].x[1] );
Planets[j].a[2] += a * Planets[i].mass * ( Planets[i].x[2] - Planets[j].x[2] );
//printf("%lf\n", Planets[i].a[1]);
}
}
if (run %(long int) pow(10,7)== 0) {
//fprintf(output, "%d %lf %lf %lf\n", index, Planets[index].x[0], Planets[index].x[1], Planets[index].x[2]);
fprintf(output, "\n" );
}
for (int m = 0; m < noOfBodies; m++) {
Planets[m].v[0] += 0.5 * Planets[m].a[0]* dt;
Planets[m].v[1] += 0.5 * Planets[m].a[1]* dt;
Planets[m].v[2] += 0.5 * Planets[m].a[2]* dt;
}
run++;
}
printf("energy before: %.lf\n", energy(noOfBodies, G));
printf("distance: %lf\n", distance(Planets[0].x, Planets[1].x));
fclose(output);
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment