ORBITA DE LA TIERRA EN TORNO AL SOL

El siguiente programa muestra la simulación de la trayectoria de la tierra en torno al Sol con ordenador. Se dice que está simplificado ya que se considera al sol un objeto inmóvil y se aplican ecuaciones de gravitación de Newton.
La trayectoria se da en forma de coordenadas XY, así como la distancia del sistema TIERRA-SOL que debe mantenerse constante.

El programa esta implementado para un dt de un segundo, lo cual hace que el programa se ejecute lento.

La diferencia entre la posición inicial y final (que deberían ser iguales ya que el programa calcula la trayectoria de todo un año) se debe al propio error del calculo numérico  así como los errores en las magnitudes de las masas del sol y de la tierra, errores también en la velocidad inicial de la tierra así como en la distancia inicial y la constante  de gravitación universal G. De todos modos no queda muy mal, y se ve más o menos cual es la trayectoria que describe. Espero que les sirva:

#include <stdio.h>
#include <math.h>


const int L=31536000;/*SEGS. EN UN AÑO*/
const double dt=1;/*DT ES UN SEG.*/
const double m1=5.97e24, m2=1.99e30, G= 6.674e-11;

int i;
double F,coseno,seno,d;

struct estado{
    double X;
    double Y;
    double Vx;
    double Vy;
    double Ax;
    double Ay;
}Solido1[L];

//Situacion inicial;   


int main()
{
    Solido1[0].X=1.5e11;
    Solido1[0].Y=0;
    Solido1[0].Vx=0;
    Solido1[0].Vy=2.98e4;
   
    d=sqrt(pow(Solido1[0].X,2)+pow(Solido1[0].X,2));
    F=(G*m1*m2)/(pow(d,2));
   
    Solido1[0].Ax=-F/m1;
    Solido1[0].Ay=0;
   
    printf("\n\n     ORBITA DE LA TIERRA EN TORNO AL SOL\n\n\n");
   
    for (i=1;i<=L;i++)
    {
        Solido1[i].X=Solido1[i-1].X+Solido1[i-1].Vx*dt;
        Solido1[i].Y=Solido1[i-1].Y+Solido1[i-1].Vy*dt;
       
        d=sqrt(pow(Solido1[i].X,2)+pow(Solido1[i].Y,2));
        coseno=Solido1[i].X/d;
        seno=Solido1[i].Y/d;
       
        Solido1[i].Vx=Solido1[i-1].Vx+Solido1[i-1].Ax*dt;
        Solido1[i].Vy=Solido1[i-1].Vy+Solido1[i-1].Ay*dt;
       
        F=(G*m1*m2)/(pow(d,2));
        Solido1[i].Ax=-F*coseno/m1;
        Solido1[i].Ay=-F*seno/m1;
       
        if ((i%86400)==1){/*PARA NO IMPRIMIR DEMASIDADOS DATOS*/
            printf(" DIA %4d)  X=%10.4f MKm",i/86400,Solido1[i].X/1000000000);
            printf("   Y=%10.4f MKm",Solido1[i].Y/1000000000);
            printf("   d=%10.4f MKm\n",d/1000000000);
        }   
    }
   
    printf("\n\n MKm significa MegaKilometro o Millones de kilometros.");
    return 0;
}
   

Un Saludo.