相互作用を表したいのに一つの質点からの力のみになってしまう。

[ちむ]

コード:

#include <stdlib.h>
#include<stdio.h>
#include <GL/glut.h>
#include <math.h>

#define h 10.
#define N 9
#define a 0.001
#define g  9.8;


int width = 300, height = 600;

/* 初期値 */
float y[N] = {10.,10,10,0,0,0,-10,-10,-10};
float Vy[N] ={ 0.,0.,0.,0.,0.,0.,0.,0.,0.};
float t =1;
float dt =1;
float tmax=10.;
float x[N] = {-10,0.,10,-10,0,10,-10,0,10};
float Vx[N] = {0.,0.,0.,0.,0.,0.,0.,0.,0};
float a1=4/(3*powf(h,8.));
float a2=-24/(3*powf(h,8.));
float p=1.  ;

/*距離の計算*/
void distance(float x[],float y[], float dis[][N]){
	int i,j;
	for(i=0; i<N; i++){
		for(j=0; j<N;j++){
			dis[i][j]=sqrtf( powf(x[j]-x[i],2)+powf(y[j]-y[i],2)); 
		}
	}
}
/*位置ベクトルの差*/
void residual_vector(float x[], float y[], float res_x[][N],float res_y[][N]){
	int i,j;
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
		res_x[i][j]=x[i]-x[j];
		res_y[i][j]=y[i]-y[j];
		}
	}
}
/*運動方程式*/
float  *Fx(float x[] , float y[]){
	float dis[N][N],res_x[N][N],res_y[N][N];
	static float Fx[N];
	int j,i;
	distance(x,y, dis);
	residual_vector(x,y,res_x,res_y);
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			if (i==j)
				Fx[i]=0;
			else
				Fx[i]=(a/dis[i][j])*res_x[i][j];
		}
	}
	return Fx ;
}
float * Fy(float x[] , float y[]){
	float dis[N][N],res_x[N][N],res_y[N][N];
	static float Fy[N];
	int j,i;
	distance(x,y, dis);
	residual_vector(x,y,res_x,res_y);
	for(i=0; i<N; i++){
		for(j=0; j<N; j++){
			if (i==j) 
				Fy[i]=0;
			else
			    Fy[i]=(a/dis[i][j])*res_y[i][j];
		}
	}
	return Fy ;
}


void Calculate(  )
{    
    float kx2[N], kx1[N], ky1[N],ky2[N],kx3[N],ky3[N],kx4[N],ky4[N];
	int i;
	for(t=0; t<tmax; t+=dt){
		for(i=0;i<N;i++) {
                kx1[i]=dt**(Fx(x,y)+i);
                ky1[i]=dt**(Fy(x,y)+i);
				
				x[i]=x[i]+kx1[i]/2.0;
				y[i]=y[i]+ky1[i]/2.0;
                kx2[i]=dt**(Fx(x,y)+i);
				ky2[i]=dt**(Fy(x,y)+i);

				x[i]=x[i]+kx2[i]/2.0-ky1[i]/2.0;
				y[i]=y[i]+ky2[i]/2.0-ky1[i]/2.0;
                kx3[i]=dt**(Fx(x,y)+i);
                ky3[i]=dt**(Fy(x,y)+i);

				x[i]=x[i]+kx3[i]/2.0-kx2[i]/2.0;
				y[i]=y[i]+ky3[i]/2.0-kx2[i]/2.0;

                kx4[i]=dt**(Fx(x,y)+i);
                ky4[i]=dt**(Fy(x,y)+i);

                Vx[i]=Vx[i]+(kx1[i]+2.0*kx2[i]+2.0*kx3[i]+kx4[i])/6.0;
				Vy[i]=Vy[i]+(ky1[i]+2.0*ky2[i]+2.0*ky3[i]+ky4[i])/6.0;
				x[i]=x[i]-kx3[i]/2.0-kx2[i]/2.0;
				y[i]=y[i]-ky3[i]/2.0-kx2[i]/2.0;
		}
   }
}

void FreeFall(void){ 
	int i;
   /* axis */
  glColor3f(.2, .3, .5);
 
  glEnd();

  /* draw points */
  glColor3f( .9, .2, .2 );
  glPointSize(5.);
  glBegin( GL_POINTS );
  for(i=0;i<N;i++)
    glVertex2f( x[i], y[i]);
  glEnd();

 
}

void display(void){
  glClearColor( .0, .0, .0, .0);
  glClear( GL_COLOR_BUFFER_BIT );
  FreeFall();
  glutSwapBuffers();
}

void idle(void){
  float dis[N][N],res_x[N][N],res_y[N][N];
  int i;
    Calculate();
  for(i=0; i<N; i++){
	y[i] += Vy[i]*dt;
	x[i] += Vx[i]*dt;

	distance(x,y,dis);
	residual_vector(x,y,res_x,res_y);
	

	printf("%f,%f,%f\n",dis[i][0], dis[i][1],res_x[i][0]);
  
 
  
	/* 衝突 */
	if(x[i]>=100.){
	Vx[i]=-Vx[i];
	}
	if(y[i]<-100.){
	  Vy[i]=-Vy[i];
	}
	if(x[i]<-100.){
	  Vx[i]=-Vx[i];
	}
	if(y[i]>=100.){
	  Vy[i]=-Vy[i];
	}
	glutPostRedisplay();
  }
}

int main(int argc, char** argv){
  /* init gult */
  glutInit( &argc, argv );
  glutInitDisplayMode( GLUT_DOUBLE | GLUT_RGB );
  glutInitWindowPosition( 0, 0 );
  glutInitWindowSize( width, height );
  glutCreateWindow( "FreeFall" );

  /* call-back */
  glutDisplayFunc( display );
  glutIdleFunc( idle );

  /* 2D-setup */
  glMatrixMode( GL_PROJECTION );
  glLoadIdentity();
  gluOrtho2D( -100., 100., -100., 100.);

  glutMainLoop();
  
  return 0;
}

互いの相互作用ではなくx[9]番目とほかの質点との相互作用になってしまう　当方プログラミングは独学でむちゃくちゃです　恐らくcalclate関数の中身が問題なのだと思います。

[usao]

あるタイムステップにおいて
ある質点iに対して加わる力は，他質点群との間の粒子間力の総和となるべきです．
問題の現象は，関数Fx(),Fy()の計算結果が上書き代入になっていることが原因と思われます．
(line57等)