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



//define una estructura de parametros para usar durante todo el programa
struct Par{
  double gamma, mu,l;
} aa;


// Ecuaciones del sistema a resolver
void ecuaciones(int n, double v[], double dv[], double t){
  double r,sigma,tita;
  r=v[0];sigma=v[1], tita = v[2];

  dv[0]= sigma; //de la parte radial
  dv[1]= -aa.gamma/(aa.mu*r*r)+aa.l*aa.l/(aa.mu*aa.mu*r*r*r); //de la parte radial
  dv[2] = aa.l/(aa.mu*r*r); //de la parte angular

  return;
}


//Programa principal
int main(){

  int i,j;
  FILE *ptr;
  double v[4],t,dt,t_pre,t_max,x,y;

  //Archivo de salida
  ptr=fopen("orb_k.dat","w");

//tiempo de integración y paso de integración
  dt=0.01;
  t_max=1000;

 
//valores de los parámetros

  aa.mu = 0.8;
  aa.l = 1600;   //
  aa.gamma= 100000;


 //Condiciones Iniciales
  v[0]=aa.l*aa.l/(aa.mu*aa.gamma);  // esta es la condicion inicial r_0=l^2/mu*gamma, para la orbita circular
  v[1]=1.;   // para que haga orbita circular la c.i. debe ser 0, despues recordar E = 1/2mu((d(r_0)/dt)^2+r_0^2*theta_punto^2)-gamma/r_0) y se conserva!!



//empieza a integrar
  t=0.;
  while(t<t_max){
    //Integra las ecuaciones utilizando el metodo de Runge Kutta
    rk4(ecuaciones,v,4,t,dt);

  //define las variables en cartesianas
  x = v[0]*cos(v[2]); 
  y = v[0]*sin(v[2]); 

   //Imprime la integracion en un vector

   fprintf(ptr,"%lg\t%lg\t%lg\t%lg\n",t,x,y,v[0]);

    t+=dt;
  }
  fclose(ptr);
  return(0);
}