Xzin
Messages postés2Date d'inscriptiondimanche 13 mars 2011StatutMembreDernière intervention13 mars 2011
-
13 mars 2011 à 13:10
Xzin
Messages postés2Date d'inscriptiondimanche 13 mars 2011StatutMembreDernière intervention13 mars 2011
-
13 mars 2011 à 13:16
Bonjour, j'aurai besoin d'un peu d'aide pour mon projet d'informatique.
Mon thème est la capture d'une planete par une etoile qui s'approche. Je dois résoudre ce problème à 3 corps complet avec RK4. Et j'affiche le résultat sous gnuplot.
J'ai dejà programmé en C++ la partie différentielle, affichage gnuplot mais mon probleme est sur les Conditions initiales de l'approche de l'étoile.
Voici mon programme :
#include
#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;}
}
