SOLUTION GRAPHIQUE APPROCHÉE AU PROBLÈME DES N CORPS

/* Solution graphique approchée du problème des N corps */
/* Metaldwarf 2005 */
/* Testé sous Linux avec la libSDL 1.2 */
/* gcc ncorps.c -o ncorps -lm -Wall `sdl-config --cflags --libs` */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <SDL/SDL.h>

#define SDL_VIDEO_FLAGS (SDL_HWSURFACE | SDL_ANYFORMAT)

/* Le nombre de corps */
#define N	5

/************************ Partie Algorithmique ************************/

/* La constante universelle de gravitation */
#define G 6.67e-11

/* vecteur du plan */
typedef struct _vect
{
	double x;
	double y;
} vect;

vect vecteur_nul = {0,0};

typedef struct _corps
{
	double m;	/* masse */
	vect pos;	/* position */
	vect v;		/* vitesse */
	vect a;		/* acceleration */
} corps;	

/* norme d'un vecteur */
double norme2(vect u)
{
	return sqrt((u.x * u.x) + (u.y * u.y));
}

/* res = u1 + u2 */
vect vect_add(vect u1, vect u2)
{
	vect res;
	res.x = u1.x + u2.x;
	res.y = u1.y + u2.y;
	return res;
}

/* res = u1 - u2 */
vect vect_sub(vect u1, vect u2)
{
	vect res;
	res.x = u1.x - u2.x;
	res.y = u1.y - u2.y;
	return res;
}

/* res = a.u  (a réel) */
vect mul_scal(double a, vect u)
{
	vect res;
	res.x = a * u.x;
	res.y = a * u.y;
	return res;
}

/* force gravitationnelle exercée par p1 sur p2 */
vect force_grav(corps p1, corps p2)
{
	double d;	/* distance entre les deux corps */
	vect u;		/* vecteur unitaire dirigé de p2 vers p1 */
	
	d = norme2(vect_sub(p1.pos, p2.pos));
	u = mul_scal((1 / d),vect_sub(p1.pos,p2.pos));
	return mul_scal(G * (p1.m * p2.m) * (1/(d*d)), u);
}

/* application du principe fondamental de la dynamique au point p, d'index i dans un tableau de taille p */
/* modifie son accélération */
void pfd(corps *tab, int i, int n)
{
	int k;
	vect f = vecteur_nul;	/* vecteur force */
	corps p = tab[i];
	
	for(k=0;k<n;k++)
	{
		if(k != i)
			f = vect_add(f, force_grav(tab[k], p));
	}
	tab[i].a = mul_scal((1 / p.m),f);
}

/* optimisation : F(1->2) = -F(2->1) donc on peut éviter la moitie des calculs */
void pfd2(corps *tab, int i, int n)
{
	int k;
	corps p = tab[i];
	vect v = vecteur_nul;
	
	/* premier corps : on reinitialise l'acceleration de tous les corps */
	if(i==0)
		for(k=0;k<n;k++)
			tab[k].a = vecteur_nul;
	/* Dernier corps : toutes les contributions ont deja ete calculee */
	if(i==(n-1)) return;
	
	for(k=i+1;k<n;k++)
	{
		v = force_grav(tab[k],p);
		tab[i].a = vect_add(tab[i].a,mul_scal((1/p.m),v));
		tab[k].a = vect_add(tab[k].a,mul_scal((-1/tab[k].m),v));
	}
 }

/* on a : dv/dt = a    et   dpos/dt = v 
   donc : dv = a * dt  et   dpos = v * dt */
/* modifie la position d'un corps (on considère que l'accélération a déja été calculée */
void change_pos(corps *tab, int i, double dt)
{
	tab[i].v = vect_add(tab[i].v , mul_scal(dt,tab[i].a));
	tab[i].pos = vect_add(tab[i].pos , mul_scal(dt,tab[i].v));
}

/* met à jour le tableau de corps en avancant de dt dans le temps */
/* NB : Les deux boucles doivent etre laisees comme ceci sinon l'acceleration des corps n'est pas calculee correctement */
/* En effet certains corps on deja bouges avant de calculer l'acceleration des suivants */
void maj(corps *tab, int n, double dt)
{
	int i;
	for(i=0;i<n;i++)	/* on change les accélérations des corps */
	{
		pfd2(tab,i,n);
	}
	for(i=0;i<n;i++)	/* puis on met à jour les positions */
	{
		change_pos(tab,i,dt);
	}
}

/************************  Partie Graphique  ************************/ 

#define NB_COLOR	5
#define XRES		1024
#define YRES		768
Uint32 color[NB_COLOR] = {0x00ff0000, 0x0000ff00, 0x000000ff, 0x00ff00ff, 0x00ffff00};

/* repris de Linux Magazine 65 (p.80) */
void PutPixel(SDL_Surface *surface, Uint16 x, Uint16 y, Uint32 color)
{
	Uint8 bpp = surface->format->BytesPerPixel;
	Uint8 *p = ((Uint8 *) surface->pixels) + y * surface->pitch + x * bpp;
	
	switch(bpp)
	{
		case 1:
			*p = (Uint8) color;
			break;
		case 2:
			*(Uint16 *)p = (Uint16) color;
			break;
		case 3:
			if(SDL_BYTEORDER == SDL_BIG_ENDIAN)
			{
				*(Uint16 *)p = ((color >> 8) & 0xff00) | ((color >> 8) & 0xff);
				*(p+2) = color & 0xff;
			}
			else
			{
				*(Uint16 *)p = color & 0xffff;
				*(p+2) = (color >> 16) & 0xff;
			}
			break;
		case 4:
			*(Uint32 *)p = color;
			break;
	}
}

/* Anciennes positions des corps dans le systeme de coordonnees de l'ecran */
Uint16 old_pos_x[N];
Uint16 old_pos_y[N];

void animer(corps *tab, int n, double dt, SDL_Surface *screen)
{
	int i;
	int need_update = 0;
	while(1)
	{
		maj(tab,n,dt);

		for(i=0;i<n;i++)
		{	
			Uint16 new_x = (Uint16)tab[i].pos.x;
			Uint16 new_y = (Uint16)tab[i].pos.y;
			if((new_x <= 0) || (new_x >= XRES)
			 ||(new_y <= 0) || (new_y >= YRES))
			{
			 	printf("Un corps a quitté l'écran. Fin du programme\n");
				return;
			}
			/* On ne change l'affichage que si le corps a effectivement bouge a l'ecran */
			if((new_x != old_pos_x[i]) || (new_y != old_pos_y[i]))
			{
				old_pos_x[i] = new_x;
				old_pos_y[i] = new_y;
				SDL_LockSurface(screen);
				PutPixel(screen,new_x, new_y, color[i % NB_COLOR]);
				SDL_UnlockSurface(screen);
				need_update++;
			}
		}
		/* Pour pallier a la lenteur de l'affichage on ne le rafraichit que toutes les 10 modifications */
		if(need_update >= 10)
		{
			need_update = 0;
			SDL_Flip(screen);
		}	
	}
}

int main(int argc, char **argv)
{
	SDL_Surface *screen;
	corps tab[]= {	{1000000000 ,{512,384},	vecteur_nul	,vecteur_nul},		/* le "soleil" */
			{1000000    ,{350,384},	{0,0.01}	,vecteur_nul},		/* les planetes */
			{3000000    ,{512,450},	{0.02,0} 	,vecteur_nul},
			{1000500    ,{300,300},	{0.007,-0.001}	,vecteur_nul},
			{2000000    ,{650,450},	{-0.003,0.003}	,vecteur_nul}};

	if(SDL_Init(SDL_INIT_VIDEO) == -1)
	{
		printf("Initialisation de SDL impossible : %s\n",SDL_GetError());
		return -1;
	}
	screen = SDL_SetVideoMode(XRES,YRES,24,SDL_VIDEO_FLAGS);
	SDL_FillRect(screen, NULL, SDL_MapRGB(screen->format, 0x00,0x00,0x00));	// écran noir
	SDL_ShowCursor(SDL_DISABLE);
	SDL_Flip(screen);
	
	animer(tab,N,1,screen);
	
	SDL_Quit();	
	return 0;
}