ТЕМА:

#include <stdlib.h>
#include <math.h>
#include <GL/glut.h>
#define NMAX 1000
#define PMAX 10000
int N=10;
int K=3;
int T=1;
int P=1000;
int Width=1000, Height=800;
double K1,K2,K3,V1;
double R[NMAX], VX[NMAX], VY[NMAX], XX[NMAX], YY[NMAX];
double X[PMAX], Y[PMAX]; 
int L[NMAX], Z[NMAX];
int T2;

void nbodyinit( ) {
static int i;
	K1=1.0*K;
	T2=T*K*K;
	V1=1/sqrt(K1);	
	R[0]=3.0; 
	VX[0]=0;
	VY[0]=0;		
	XX[0]=0;
	YY[0]=0;
	L[0]=0;	
	Z[0]=1;	
	for(i=1;i<N;i++) {
		R[i]=(rand()%1000+1)*0.001; 
		VX[i]=(rand()%1000+1)*0.0003*V1;
		VY[i]=(rand()%1000+1)*0.0003*V1;		
		XX[i]=(rand()%400-200)*K1;
		YY[i]=(rand()%400-200)*K1;
		if(XX[i]<0) VY[i]=-VY[i];	
		if(YY[i]>0) VX[i]=-VX[i];
		L[i]=i; 		
		Z[i]=1;		
	}
}

void nbodystep( ) {
static double x,y,ai,aj,r1,r2,m,mi,mj;	
static int i, j;
	for(i=0;i<(N-1);i++) for(j=i+1;j<N;j++) {	
		x=XX[j]-XX[i];
		y=YY[j]-YY[i];		
		r1=R[i]+R[j];		
		r2=x*x+y*y;
		if(Z[i] && Z[j]) {
			if( r2 > r1*r1 ) {
				r1=1/sqrt(r2);
				ai=R[j]*r1; ai=ai*ai*ai;
				VX[i]+=ai*x;
				VY[i]+=ai*y;
				aj=R[i]*r1; aj=aj*aj*aj; 
				VX[j]-=aj*x;
				VY[j]-=aj*y;				
			} else { 
				L[j]=L[i]; 			
				Z[j]=2; 
			}
		}	
	}
	for(i=0;i<N;i++) if(Z[i]) { 
		XX[i]+=VX[i];
		YY[i]+=VY[i];	
	}
	for(i=N-1;i>=0;i--) if(Z[i]==2) { 
		j=L[i];
		mi=R[i]; mi=mi*mi*mi;
		mj=R[j]; mj=mj*mj*mj;
		m=mi+mj;
		VX[j]=(VX[i]*mi+VX[j]*mj)/m;
		VY[j]=(VY[i]*mi+VY[j]*mj)/m;			
		XX[j]=(XX[i]*mi+XX[j]*mj)/m;
		YY[j]=(YY[i]*mi+YY[j]*mj)/m;		
		R[j]=pow(m,1.0/3.0);		
		Z[i]=0;				
	}
}

void idle( void ) {
	nbodystep( );	
	glutPostRedisplay( );
}

void Display( void ) {
	static int i, t, p;	
	static double k,x,y;
	glClearColor( 0, 0, 0, 1 );
	glClear( GL_COLOR_BUFFER_BIT );
	glEnable(GL_POINT_SMOOTH);
	glColor3f( 0.5, 0.5, 0.5 );
	glPointSize( 1 );	
	for(i=0;i<P;i++) {
		glBegin( GL_POINTS );
		glVertex2f( X[i], Y[i] );
		glEnd( );			
	}
	t=(t+1)%T2;
	k=1/K1;
	glColor3f( 1, 1, 0 );
	for(i=0;i<N;i++) {
		if(i == 1) glColor3f( 0, 1, 0 );	
		if(Z[i]) { 
			glPointSize( 2.0+2.0*R[i] );
			glBegin( GL_POINTS );
			x=XX[i]*k+(Width/2); y=YY[i]*k+(Height/2);
			glVertex2f( x, y );
			if((t==0)&&(x>0)&&(x<Width)&&(y>0)&&(y<Height))		
				{ p=(p+1)%P; X[p]=x; Y[p]=y; }
			glEnd( );		
		}
	}			
	glFlush( );
}

void Reshape( int w, int h ) {
	static int i;
	Width = w; Height = h;
	glViewport( 0, 0, w, h );
	glMatrixMode( GL_PROJECTION );
	glLoadIdentity( );
	glOrtho( 0, w, 0, h, -1, 1 );
	for(i=0;i<P;i++) { X[i]=0; Y[i]=0; }
}

void Keyboard( unsigned char key, int x, int y ) {
	const char ESCAPE = '\033';
	if( key == ESCAPE ) exit( 0 ) ;
}

void main( int argc, char* argv[ ] ) {
	static int i,m;
	i=argc; while(i--) {
		if(argv[i][0]=='n') sscanf(argv[i],"n%d", &N);
		if(argv[i][0]=='k') sscanf(argv[i],"k%d", &K);
		if(argv[i][0]=='t') sscanf(argv[i],"t%d", &T);
		if(argv[i][0]=='p') sscanf(argv[i],"p%d", &P);
	}	
	printf("N=%d\n",N);
	printf("K=%d\n",K);	
	printf("T=%d\n",T);
	printf("P=%d\n",P);
	nbodyinit( );
	glutInit( &argc, argv );
	glutInitDisplayMode( GLUT_RGB );
	glutInitWindowSize( Width , Height );
	glutInitWindowPosition( 50, 50 );		
	glutCreateWindow( "Test" );
	glutDisplayFunc( Display );
	glutReshapeFunc( Reshape );
	glutKeyboardFunc( Keyboard );
	glutIdleFunc( idle );
	glutMainLoop( );
}