diff --git a/src/jupiter.c b/src/jupiter.c
new file mode 100644
index 0000000000000000000000000000000000000000..7cbe00aa5c19e9012139a9f4c1806224eac3d08e
--- /dev/null
+++ b/src/jupiter.c
@@ -0,0 +1,166 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#define pi 3.141592653589793
+#define solar_mass 1
+#define G (4 * pi * pi)
+#define days_per_year 365.24
+#define JD 5.2
+#define VJ0 2.75535903023
+
+
+typedef struct {
+ double x[3];
+ double v[3];
+ double a[3];
+ double mass;
+} Planet;
+
+//Unit scaling: Mass in solar mass, distance in AU, velocity in AU/year, and constant G in AU^3/solar mass*year^2
+
+Planet Planets[] = {
+ {
+ {0,0,0},
+ {0,0,0},
+ {0,0,0},
+ solar_mass
+ },
+ {
+ {JD, 0, 0},
+ {0, VJ0, 0},
+ {0,0,0},
+ 9.54791938424326609e-04 * solar_mass
+ },
+ {
+ {-JD, 0, 0},
+ {0, -VJ0, 0},
+ {0, 0, 0},
+ 2.85885980666130812e-04 * solar_mass
+ }
+};
+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 + 0.001);
+
+}
+
+void update(int steps){
+ long run = 0;
+ const int noOfBodies = sizeof(Planets)/sizeof(Planet);
+ const double dt = 0.01;
+
+ 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;
+ }
+ 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;
+ }
+
+ 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);
+
+ 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] );
+ }
+ }
+
+ 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++;}
+}
+
+int main() {
+ FILE *starting = fopen("starting.txt", "w+");
+ if (starting == NULL) {
+ printf("HELP\n COULDNT OPEN FILE\nSIGSEVV");
+ } else {
+ printf("Output File opened\n");
+ }
+ const int noOfBodies = sizeof(Planets)/sizeof(Planet);
+ double s =0;
+ double ds=0.001;
+
+ for(s=1.5; s<=3; s+=ds){
+
+ Planets[2].x[0]=-s*JD;
+ Planets[2].x[1]=0;
+ Planets[2].x[2]=0;
+ Planets[2].v[0]=0;
+ Planets[2].v[1]=-sqrt(1/s)*VJ0;
+ Planets[2].v[2]=0;
+ for(int i=0; i<3;i++){
+ for(int k=0; k<3; k++){
+ Planets[k].a[i]=0;
+ }
+ }
+ for(int i=0; i<3; i++){
+ Planets[0].x[i]=0;
+ Planets[0].v[i]=0;
+ }
+ Planets[1].x[0]=JD;
+ Planets[1].x[1]=0;
+ Planets[1].x[2]=0;
+ Planets[1].v[0]=0;
+ Planets[1].v[1]=VJ0;
+ Planets[1].v[2]=0;
+ update(10000);
+ double r10= sqrt(Planets[2].x[1]*Planets[2].x[1]+Planets[2].x[0]*Planets[2].x[0]);
+ update(10000);
+ double r1= sqrt(Planets[2].x[1]*Planets[2].x[1]+Planets[2].x[0]*Planets[2].x[0]);
+
+ Planets[2].x[0]=(-s*JD)-0.01;
+ Planets[2].x[1]=0;
+ Planets[2].x[2]=0;
+ Planets[2].v[0]=0;
+ Planets[2].v[1]=-sqrt(G/((s*JD)-0.01));
+ Planets[2].v[2]=0;
+ for(int i=0; i<3;i++){
+ for(int k=0; k<3; k++){
+ Planets[k].a[i]=0;
+ }
+ }
+ for(int i=0; i<3; i++){
+ Planets[0].x[i]=0;
+ Planets[0].v[i]=0;
+ }
+ Planets[1].x[0]=JD;
+ Planets[1].x[1]=0;
+ Planets[1].x[2]=0;
+ Planets[1].v[0]=0;
+ Planets[1].v[1]=VJ0;
+ Planets[1].v[2]=0;
+ update(10000);
+ double r20= sqrt(Planets[2].x[1]*Planets[2].x[1]+Planets[2].x[0]*Planets[2].x[0]);
+ update(10000);
+ double r2=sqrt(Planets[2].x[1]*Planets[2].x[1]+Planets[2].x[0]*Planets[2].x[0]);
+ double d0=fabs(r10-r20);
+ double d=fabs(r2-r1);
+ double l= log(d/d0);
+ fprintf(starting, "%lf %lf\n", r10, l);
+
+ }
+ fclose(starting);
+}
+