/* ----------------------------------------*/ /* Xavier Gnata -- Laure Gonnord, oct 2004 */ /*gcc -Wall -ansi -lm -lX11 -lgraphdsu 3BM.c -o 3BM -L . -L /usr/X11R6/lib */ /* ----------------------------------------*/ #include #include #include #include "graphdsu.h" /* Macro pour calculer les coord graphiques */ #define GraphX(x) (GraphSizeX*(x)/(6e8) + GraphSizeX/2) #define GraphY(y) (GraphSizeY*(y)/(6e8) + GraphSizeY/2) /*calcul des denominateurs*/ double Dist (double X1, double X2, double Y1, double Y2) { return pow ((X1 - X2) * (X1 - X2) + (Y1 - Y2) * (Y1 - Y2), 1.5); } /*fonctions f_i*/ double f (double GM0, double GM1, double GM2, double X[3], double Y[3], double Vx[3], double Vy[3], int i) { switch (i) { case 0: return Vx[0]; case 1: return Vx[1]; case 2: return Vx[2]; case 3: return Vy[0]; case 4: return Vy[1]; case 5: return Vy[2]; case 6: return -GM1 * (X[0] - X[1]) / Dist (X[1], X[0], Y[1], Y[0]) - GM2 * (X[0] - X[2]) / Dist (X[2], X[0], Y[2], Y[0]); case 7: return -GM0 * (X[1] - X[0]) / Dist (X[1], X[0], Y[1], Y[0]) - GM2 * (X[1] - X[2]) / Dist (X[2], X[1], Y[2], Y[1]); case 8: return -GM0 * (X[2] - X[0]) / Dist (X[2], X[0], Y[2], Y[0]) -GM1 * (X[2] - X[1]) / Dist (X[2], X[1], Y[2], Y[1]); case 9: return -GM1 * (Y[0] - Y[1]) / Dist (X[1], X[0], Y[1], Y[0]) - GM2 * (Y[0] - Y[2]) / Dist (X[2], X[0], Y[2], Y[0]); case 10: return -GM0 * (Y[1] - Y[0]) / Dist (X[1], X[0], Y[1], Y[0]) - GM2 * (Y[1] - Y[2]) / Dist (X[2], X[1], Y[2], Y[1]); case 11: return -GM0 * (Y[2] - Y[0]) / Dist (X[2], X[0], Y[2], Y[0]) -GM1 * (Y[2] - Y[1]) / Dist (X[2], X[1], Y[2], Y[1]); default: printf ("Indice hors limite. f_%d n'existe pas!\n", i); exit (-1); } } /* La méthode de RK ordre 4 integation equas diffs*/ void rk (double GM0, double GM1, double GM2, double X[3], double Y[3], double Vx[3], double Vy[3], double h) { int i, i_rk; double TempX[3], TempY[3], TempVx[3], TempVy[3]; double k1[12], k2[12], k3[12], k4[12]; for (i_rk = 0; i_rk < 12; i_rk++) k1[i_rk] = h * f (GM0, GM1, GM2, X, Y, Vx, Vy, i_rk); for (i = 0; i < 3; i++) { TempX[i] = X[i] + k1[i] / 2.; TempY[i] = Y[i] + k1[i + 3] / 2.; TempVx[i] = Vx[i] + k1[i + 6] / 2.; TempVy[i] = Vy[i] + k1[i + 9] / 2.; } for (i_rk = 0; i_rk < 12; i_rk++) k2[i_rk] = h * f (GM0, GM1, GM2, TempX, TempY, TempVx, TempVy, i_rk); for (i = 0; i < 3; i++) { TempX[i] = X[i] + k2[i] / 2.; TempY[i] = Y[i] + k2[i + 3] / 2.; TempVx[i] = Vx[i] + k2[i + 6] / 2.; TempVy[i] = Vy[i] + k2[i + 9] / 2.; } for (i_rk = 0; i_rk < 12; i_rk++) k3[i_rk] = h * f (GM0, GM1, GM2, TempX, TempY, TempVx, TempVy, i_rk); for (i = 0; i < 3; i++) { TempX[i] = X[i] + k3[i]; TempY[i] = Y[i] + k3[i + 3]; TempVx[i] = Vx[i] + k3[i + 6]; TempVy[i] = Vy[i] + k3[i + 9]; } for (i_rk = 0; i_rk < 12; i_rk++) k4[i_rk] = h * f (GM0, GM1, GM2, TempX, TempY, TempVx, TempVy, i_rk); for (i = 0; i < 3; i++) { X[i] = X[i] + (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) / 6.; Y[i] = Y[i] + (k1[i + 3] + 2 * k2[i + 3] + 2 * k3[i + 3] + k4[i + 3]) / 6.; Vx[i] = Vx[i] + (k1[i + 6] + 2 * k2[i + 6] + 2 * k3[i + 6] + k4[i + 6]) / 6.; Vy[i] = Vy[i] + (k1[i + 9] + 2 * k2[i + 9] + 2 * k3[i + 9] + k4[i + 9]) / 6.; } } int main (int argc, char *argv[]) { /* Fichier pour stocker les résultats */ FILE *file; /* Variable de travail */ int i; /* Taille de la fenetre graphique */ int GraphSizeX = 500; int GraphSizeY = 500; /* Positions */ double X[3]; double Y[3]; /* Vitesses projetées sur X et Y */ double Vx[3]; double Vy[3]; /* Pour tracer les tajectoires, il nous faut le point n-1 et le point n. De même pour les vitesses */ double OldX[3]; double OldY[3]; /* Pas de temps en seconde */ int h = 100; /* Nombre de pas de temps écoulés */ int NbDeltaT = 0; /* Nombre de pas de temps MAX */ int MaxNbDeltaT; /* CI : Mars et deux de ses satellites. Les CI sont telles que les deux satellites sur des orbites circulaires G=6.67259e-11 et MO = 6.4185e23 : on va pas coder le produit... */ double GM0 = 6.67259 * 6.4185e12; double GM1 = 6.67259 * 1.27e5; double GM2 = 6.67259 * 1.8e4; X[0] = 0; X[1] = 3.752e7; X[2] = -2.346e8; Y[0] = 0; Y[1] = 0; Y[2] = 0; Vx[0] = 0; Vx[1] = 0; Vx[2] = 0; Vy[0] = 0; Vy[1] = 950; /* 1069 si on veut un cercle */ Vy[2] = -380; /* 427 si on veut un cercle */ /* Ouverture du fichier resultat */ if ((file = fopen ("3BMResults.txt", "wr")) == NULL) { printf ("Error opening file\n"); exit (1); } /* On suppose que l'utilisateur ne va pas entrer n'importe quoi... */ printf ("Combien de pas de simulation? (%d seconde par pas) : ", h); scanf ("%d", &MaxNbDeltaT); fflush (stdin); /* Affichage utilisant la lib graphdsu */ Initialiser (GraphSizeX, GraphSizeY); /* GraphInverse(); n'a pas l'air de marcher donc bourrin pour avoir un fond noir. C'est joli un fond noir. */ /* RectanglePlein (0, 0, GraphSizeX, GraphSizeY);*/ ChangeCouleur (2); /* Mars au centre. On n'est pas à l'échelle */ CerclePlein (GraphSizeX / 2, GraphSizeY / 2, 5); /* On stocke les CI pour pouvoir afficher le premier point */ for (i = 0; i < 3; i++) { OldX[i] = X[i]; OldY[i] = Y[i]; } /* On intègre et on affiche les trajectoires */ for (NbDeltaT = 1; NbDeltaT <= MaxNbDeltaT; NbDeltaT++) { printf ("%d secondes écoulées depuis le début de la simulation.\n", h * NbDeltaT); rk (GM0, GM1, GM2, X, Y, Vx, Vy, h); fprintf (file, "%g\t%g\t%g\t%g\t%g\t%g\n", X[0], Y[0], X[1], Y[1], X[2], Y[2]); ChangeCouleur (4); Ligne (GraphX (OldX[1]), GraphY (OldY[1]), GraphX (X[1]), GraphY (Y[1])); ChangeCouleur (3); Ligne (GraphX (OldX[2]), GraphY (OldY[2]), GraphX (X[2]), GraphY (Y[2])); for (i = 0; i < 3; i++) { OldX[i] = X[i]; OldY[i] = Y[i]; } } /* This is the end... */ while (!TestClavier ()); Clore (); fclose (file); return 0; }