Upload New File

src/jupiter.c

+#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);
+}
+