Capture d une etoile par une planete

#include<iostream>
#include<fstream>
#include<math.h>
#include<cstdlib>
#define PM 4
#define N 12
using namespace std;

void masse(double m[]);
void sys(double q[],double t, double qp[],int n);
double stabilite(double q[]);
void determination_CI(double q[], int np, double dt);
void val_init(double q[], double alpha, double Di, double Vi);
void calcul_bonne_trajectoire(double q[], int np, double dt);
void tracer_gnuplot();
void rk4(void(*sd)(double *,double,double *,int),double q[],double t,double dt,int n);


int main()
{
  int np;
  double q[N], dt, tfin;
  
  np=100001; // Nombre de point
  tfin=100000; // temps de fin de la simulation
  dt=tfin/(np-1); // Intervalle de temps
  
  determination_CI(q,np,dt);
  //val_init(q,0.1,150,0.0002);
  calcul_bonne_trajectoire(q,np,dt); // On calcul la bonne trajectoire et l'enregistre
  
  tracer_gnuplot(); // on trace
  
  return 0;
} 

// Déclaration masses 
double masse(int a)
{
  //En masse solaire
  double m[3];
  m[0]=1.00000597682; // Masse etoile 1 
  m[1]=0.00000354786104043; // Masse planete 1
  m[2]=1.00000597682; //Masse etoile 2
  return m[a];
}

//Systeme differentiel
void sys(double q[],double t, double qp[],int n)
{
  double G=2.96E-4;
  static double c=3./2.;
  double r1,r2,r3;
  double m[3];

  m[0]=masse(0);
  m[1]=masse(1);
  m[2]=masse(2);

  //Distance Etoile 1-Planete 1
  r1=pow((q[0]-q[8])*(q[0]-q[8])+(q[1]-q[9])*(q[1]-q[9]),c);

  //Distance Planete 1-Etoile 2
  r2=pow((q[0]-q[4])*(q[0]-q[4])+(q[1]-q[5])*(q[1]-q[5]),c);
  
  //Distance Etoile 1-Etoile 2
  r3=pow((q[4]-q[8])*(q[4]-q[8])+(q[5]-q[9])*(q[5]-q[9]),c);

  //Systeme Planete 1
  qp[0]=q[2];
  qp[1]=q[3];
  qp[2]=-G*m[0]*(q[0])/r1-G*m[2]*(q[0]-q[4])/r2;
  qp[3]=-G*m[0]*(q[1])/r1-G*m[2]*(q[1]-q[5])/r2;

  //Systeme Etoile 1
  qp[8]=q[10];
  qp[9]=q[11];
  qp[10]=-G*m[2]*(q[8]-q[4])/r3-G*m[1]*(q[8]-q[0])/r1;
  qp[11]=-G*m[2]*(q[9]-q[5])/r3-G*m[1]*(q[9]-q[1])/r1;

  //Systeme Etoile 2
  qp[4]=q[6];
  qp[5]=q[7];
  qp[6]=-G*m[0]*(q[4])/r3-G*m[1]*(q[4]-q[0])/r2;
  qp[7]=-G*m[0]*(q[5])/r3-G*m[1]*(q[5]-q[1])/r2;
  
}

void val_init(double q[],double alpha, double Di, double Vi)
{
  //Valeur en UA.
  //Valeurs initiales Planete 1
  q[0]=-30.0023653; //x1
  q[1]=7.1169847; //y1
  q[2]=-0.000865429; //vx1
  q[3]=-0.00072490; //vy1

  //Valeur initiales Etoile 1
  q[8]=0;
  q[9]=0;
  q[10]=0;
  q[11]=0;

  //Valeur initiales Etoile 2
  q[4]=Di/sqrt(2);
  q[5]=Di/sqrt(2);
  q[6]=-Vi*cos(alpha);
  q[7]=-Vi*sin(alpha);

}

//Détermination des conditions initiales pour que ça marche

void determination_CI(double q[], int np, double dt)
{
  int l; //Variable boucle
  int i=0,j=0,k=0;
  double t=0,alpha,alpha_max,Vi,Vimax,Di,Dimax,r1,r2;
  
  // Test des différentes conditions initiales 
  //On fait varier la vitesse et l'angle de l'étoile et la position initiale
  
  Di=120;
  Dimax=200;
  // Variation de la position 
  
  while (Di<=Dimax && i==0)
    {
      
      // Variation de la vitesse
      Vi=0.001;
      Vimax=0.003;
      while (Vi <= Vimax && i==0)
	{
	  // Variation de l'angle 
	  alpha=0;
	  alpha_max=M_PI/3;
	  while (alpha<=alpha_max && i==0)
	    {
	      t=0;
	      val_init(q,alpha,Di,Vi); // initialisation conditions

	      for(l=1;l<=np;l++)
		{
		  rk4(sys,q,t,dt,N);
		  t=t+dt;
		} 

	      r2=sqrt((q[0]-q[4])*(q[0]-q[4])+(q[1]-q[5])*(q[1]-q[5])); //Distance Etoile 2 - Planete
	      r1=sqrt((q[0]-q[8])*(q[0]-q[8])+(q[1]-q[9])*(q[1]-q[9])); //Distance Etoile 1 - Planete

	      if (r2<500) {j=1;}
	      k=stabilite(q);
	      i=k*j;
	      alpha=alpha+alpha_max/30;

	      cout << "Distance Etoile 1 - Planète : " << r1 << endl;
	      cout << "Distance Etoile 2 - Planète : " << r2 << endl;
	      cout << "Angle initial Etoile 2 : " << alpha << endl;
	      cout << "Vitesse initiale Etoile 2 : " << Vi << endl;
	      cout << "Distance Initiale Etoile 2 : " << Di << endl;
	      
	      if(i==1)
		{
		  cout << "Orbite stable" << endl;
		}
	      else
		{
		  cout << "Orbite instable" << endl;
		}
	    }
	  Vi=Vi+0.0005;
	}
      Di=Di+10;
    }
  
  val_init(q,alpha,Di,Vi); // On met les bonnes CI		
}

// Critère de stabilité 
double stabilite(double q[])
{ 
  int i=0; //i=0 => orbite instable, i=1 => orbite stable
  double G=2.96E-4;
  static double c=3./2.;

  double f1p, f2p, r1, r2;
  double rapport,vlib,vplan;

  double m[3];
  m[0]=masse(0);
  m[1]=masse(1);
  m[2]=masse(2);
  
  //Distance Etoile 1-Planete 1
  r1=pow((q[0]-q[8])*(q[0]-q[8])+(q[1]-q[9])*(q[1]-q[9]),c);
  
  //Distance Planete 1-Etoile 2
  r2=pow((q[0]-q[4])*(q[0]-q[4])+(q[1]-q[5])*(q[1]-q[5]),c);
  
  f1p=G*m[0]*m[1]/r1; //force etoile 1 sur planete
  f2p=G*m[1]*m[2]/r2; // force etoile 2 sur planete
  
  rapport=f1p/f2p; // Si rapport<<1 alors la force de l'étoile 1 est négligeable
  vlib=sqrt(2*G*m[2]/r2); // Vitesse libération 
  vplan=sqrt(q[2]*q[2]+q[3]*q[3]); // Vitesse de la planete 
  
  if(rapport<0.01 && vplan<vlib) {i=1;}

  cout << endl;
  cout << "Rapport :" << rapport << endl;
  cout << "Vitesse liberation :" <<  vlib << endl;
  cout << "Vitesse planete :" << vplan << endl;

  return (i);
}


//Calcul Trajectoire
void calcul_bonne_trajectoire(double q[], int np, double dt)
{
  int i;
  double t=0;
  fstream traj("traj2.dat",ios::out);
  fstream origine("origine.dat",ios::out);

  for(i=1;i<=np;i++)
    {
      traj << t << " " << q[0] << " " << q[1] << " " << q[4] << " " << q[5] << endl;
      origine << t << " " << q[8] << " " << q[9] << endl;
      rk4(sys,q,t,dt,N);
      t=t+dt;
    }
  traj.close();
  origine.close();
}


// Commandes gnuplot 
void tracer_gnuplot()
{
  fstream gnu("test.gnu",ios::out);
  gnu << "set xlabel 'Abscisse'" << endl;
  gnu << "set ylabel 'Ordonnée'" << endl;
  gnu << "plot 'traj2.dat' every 5 using 2:3 title 'Planète' w l,'traj2.dat' every 5 using 4:5 title 'Etoile' w l, 'origine.dat' every 5 using 2:3 title 'Soleil' w l" << endl;
  gnu << "pause -1" << endl;
  gnu.close();
  
  system("gnuplot test.gnu");
  
}


/* Runge Kunta */
void rk4(void(*sd)(double *,double,double *,int),double q[],double t,double dt,int n)
{
  /* Declarations et initialisations */
  int i,k,p;
  double *y[PM+1],*z[PM];
  static const double a[PM][PM]={{1./2,0,0,0},{0,1./2,0,0},{0,0,1,0},{1./6,1./3,1./3,1./6}};
  static const double b[PM]={0,1./2,1./2,1};
  /* Allocations */
  for(i=0;i<PM+1;i++) y[i]=(double *)malloc(n*sizeof(double));
  for(i=0;i<PM;i++) z[i]=(double *)malloc(n*sizeof(double));
  /* Calcul */
  for(i=0;i<n;i++) y[0][i]=q[i];
  for(p=1;p<=PM;p++)
    {
      sd(y[p-1],t+b[p-1]*dt,z[p-1],n);
      for(i=0;i<n;i++) y[p][i]=q[i];
      for(k=0;k<p;k++) {for(i=0;i<n;i++) y[p][i]=y[p][i]+dt*a[p-1][k]*z[k][i];}
    }
  for(i=0;i<n;i++) q[i]=y[PM][i];
  /* Desallocations */
  for(i=0;i<PM+1;i++) {free(y[i]); y[i]=NULL;}
  for(i=0;i<PM;i++) {free(z[i]); z[i]=NULL;}
}