PROBLEMA DE LOS DOS CUERPOS

[saraou]

Holaaa!!! alguien sabe como resolver el problema de kepler de los dos cuerpos mediante c++ y empleando el metodo de integracion de euler, para una vez conocida la fuerza, poder hallar la velocidad y luego la posicion de ambos cuerpos y trazar asi su trayectoria???
llevo ya unas semanas peleandome con este proyecto pero no soy quien a sacarlo a delante!!
muchas gracias!!!

[fary]

Hola, si se  hacerlo

saludos!

[saraou]

y me podrias, por favor,  echar una mano?

[$Edu$]

Claro, mas de uno te ayudara aca, solo deja tu intento, porque tareas no se hacen acá.

[saraou]

Esto es lo que he intentado hacer, pero no hay forma. No obtengo ningun resultado.



//El cuerpo 1 se corresponde a la Tierra. El 2, a la Luna.
// Consideramos que el centro de masas del sistema Tierra-Luna se encuentra en el origen de coordenadas.
// La Tierra se encuentra a 4668 km del CM. La Luna a 379432 km.
// Para calcular la velocidad de rotacion de ambos cuerpos multiplicamos la velocidad angular del sistema (2.67e-6) por la distancia a la que se encuentra cada cuerpo.


#include <cmath>
#include <iostream>
#include <fstream>
#include <cstdlib>

using namespace std;


   
   double dt = 1;
   const int L=200;
   double d= 3.84e8;

   double X01= -4.5e3, vx01= -8.79 ,   vx1[L],  ax1[L], X1[L], fx1[L];
   double Y01= 4.5e3,  vy01= 8.79,   vy1[L],  ay1[L], Y1[L], fy1[L];
   double q1[L], q10 = -0.785;

   double X02= 2.68e5, vx02= 715.54,   vx2[L],  ax2[L], X2[L], fx2[L];
   double Y02= 2.68e5, vy02= 715.54,   vy2[L],  ay2[L], Y2[L], fy2[L];
   double q2[L], q20 = 0.785;

   double G= 6.674e-11 , m1= 5.97e24, m2=7.35e22 ;

int main()
{ 

    ofstream output("output.dat");


for(int i=1;i<=L;i++)
     {
           
             vx1[0]=vx01;
             X1[0]=X01;
             vy1[0]=vy01;
             Y1[0]=Y01;
             q1[0]=q10;
             q1=atan(Y1/X1);

             vx2[0]=vx02;
             X2[0]=X02;
             vy2[0]=vy02;
             Y2[0]=Y02;
             q2[0]=q20;
             q2=atan(Y2/X2);



           /*OX para el primer cuerpo*/

             fx1[i-1] = ((G*m1*m2)/(sqrt(pow(d,2))))*cos(q1[i-1]);
             ax1[i-1]= fx1[i-1]/m1;
             vx1=vx1[i-1]+(ax1[i-1]*dt);
             X1=X1[i-1]+(vx1[i-1]*dt);
           
            cout<<"\n Las posiciones OX1 son:"<<X1[i-1];
         

             /*OY para el primer cuerpo*/
             
             fy1[i-1] = ((G*m1*m2)/(sqrt(pow(d,2))))*sin(q1[i-1]);
             ay1[i-1]= fy1[i-1]/m1;
             vy1=vy1[i-1]+(ay1[i-1]*dt);
             Y1=Y1[i-1]+(vy1[i-1]*dt);
           
            cout<<"\n Las posiciones OY1 son:"<<Y1[i-1];



           /*OX para el segundo cuerpo*/
             fx2[i-1] = ((G*m1*m2)/(sqrt(pow(d,2))))*cos(q2[i-1]);
             ax2[i-1]= fx2[i-1]/m2;
             vx2=vx2[i-1]+(ax2[i-1]*dt);
             X2=X2[i-1]+(vx2[i-1]*dt);
           
            cout<<"\n Las posiciones OX2 son:"<<X2[i-1];
         

             /*OY para el segundo cuerpo*/
             
             fy2[i-1] = ((G*m1*m2)/(sqrt(pow(d,2))))*sin(q2[i-1]);
             ay2[i-1]= fy2[i-1]/m2;
             vy2=vy2[i-1]+(ay2[i-1]*dt);
             Y2=Y2[i-1]+(vy2[i-1]*dt);
           
            cout<<"\n Las posiciones OY2 son:"<<Y2[i-1];     

   


     output << X1[i-1] << " " << Y1[i-1] << endl;
     output << X2[i-1] << " " << Y2[i-1] << endl;
    }
   
   
return 0 ;
}





Muchas gracias a todo aquel que me ayude!!

[saraou]

He simplificqdo el problema que quiero resolver. Ahora considero que el sol esta fijo y quiero calcular la orbita de la tierra elrededor del sol. Aqui esta:

#include <cmath>
#include <iostream>
#include <fstream>
#include <cstdlib>

using namespace std;



   double dt = 1;
   const int L=365;
   double d[L];

   double X01= 1.52e8, vx01= 0 ,   vx1[L],  ax1[L], X1[L], fx1[L];
   double Y01= 0,  vy01= 2.98e4,   vy1[L],  ay1[L], Y1[L], fy1[L];
   double q1;

 

   double G= 6.674e-11 , m2=1.98e30 , m1=5.97e24 ;


int main()
{ 

    ofstream output("output.dat");


for(int i=0;i<=L;i++)
     {
           
             vx1[0]=vx01;
             X1[0]=X01;
             vy1[0]=vy01;
             Y1[0]=Y01;
             q1=atan(Y1/X1);
             d= sqrt((pow((X1),2))+(pow((Y1),2)));



             fx1= (G*m1*m2*(cos(q1)))/(pow(d,2));
             ax1= fx1/m1;
             vx1[i+1]=vx1+(ax1*dt);
             X1[i+1]=X1+(vx1*dt);
           
           // cout<<"\n Les positions OX1 sont:"<<X1;




             
             fy1=(G*m1*m2*(sin(q1)))/(pow(d,2));
             ay1= fy1/m1;
             vy1[i+1]=vy1+(ay1*dt);
             Y1[i+1]=Y1+(vy1*dt);
           
          //  cout<<"\n Les positions OY1 sont:"<<Y1;




   


     output << X1 << " " << Y1 << endl;
   

    }
   
   
return 0 ;
}



Por que obtengo una recta y no la orbita eliptica que deberia de salir?? muchas gracias!!!

[leosansan]

En las expresiones:

Código (cpp) [Seleccionar] Expandir

X1[i+1]=X1+(vx1*dt);
.........
.........
Y1[i+1]=Y1+(vy1*dt);
que te dan la posición hay aceleración, tanto en OX como en OY y veo que usas las expresiones: x=x0+vx*t .....¿ y no faltaría: +1/2ax*t^2?. Idem para la y.

[saraou]

acabo de probar y sigue dando una recta... pero gracias!!!!

[za.asi]

La correccion de leosansan no era para que de unacurva, sino para incluir la aceleración.

[cypascal]

Te comento:
1º) No estas aplicando las ecuaciones de Kepler (conservación del momento angular), sino las de gravitación universal (Newton)
2º) No vas a obtener una trayectoria elíptica de ninguna de las maneras, ya que supones la masa de la Tierra despreciable frente a la del Sol, diciendo que este último permanece quieto. Obtendras si acaso una trayectoria circular.
3º) Considerar dt=1dia es mucho considerar. Con esta consideración no se te cerrará la trayectoria ni por asomo.

Dicho todo esto, viendo el codigo, diría que falta cambiar el ángulo en cada iteración. Tampoco se muy bien que es lo que haces con las variables iniciales, que se te inicializan en cada iteración.

Igual me pongo a hacerlo de nuevo y si eso te cuento.

Un Saludo