Harigami
ログイン
anonymous タイトルなし
No License C
コピー
#include <stdio.h>
#include <stdlib.h>
#include<math.h>

static  const  int  NMAX = 128;
static  double  xtmp[NMAX][3];
static  double k[NMAX][3], k1[NMAX][3], k2[NMAX][3], k3[NMAX][3], k4[NMAX][3];
static  double l[NMAX][3], l1[NMAX][3], l2[NMAX][3], l3[NMAX][3], l4[NMAX][3];

void  gravity( double  (*x)[3],  double *m, double  (*a)[3],double *p, const int n){
	int i,j;
	for(i=0; i<n; i++){
		a[i][0]=a[i][1]=a[i][2]=p[i]=0.0;
		for(j=0;j<n;j++){
			if(i==j) continue;
			double dx = x[i][0] - x[j][0];
			double dy = x[i][1] - x[j][1];
			double dz = x[i][2] - x[j][2];
			double r2 = dx*dx + dy* +dz*dz;
			double rinv = 1.0 / sqrt(r2);
			double mjr3inv = m[j]*rinv*rinv*rinv;
			a[i][0] -= mjr3inv * dx;
			a[i][1] -= mjr3inv * dy;
			a[i][2] -= mjr3inv * dz;
			p[i]    -= m[j] * rinv;
		}
	}
}



double  energy( double  (*v)[3],  double *m, double *p, const  int n){

  double e = 0.0;
  int i;
  for(i=0;i<n;i++){
	e += p[i] + m[i]*(v[i][0]*v[i][0]+v[i][1]*v[i][1]+v[i][2]*v[i][2])/2;  
	}
	return e;
}


void  runge2( double  (*x)[3],  double  (*v)[3],  double *m, double  (*a)[3],
	double *p, const  double dt , const  int n){
	int i,h;
	gravity(x,m,a,p,n);
	for(i=0;i<n;i++){
		for(h=0;h<3;h++){
			k1[i][h] = v[i][h]*dt;
			l1[i][h] = a[i][h]*dt;
			xtmp[i][h] = x[i][h] + (k[i][h]/2);
		}
	}

	gravity(xtmp,m,a,p,n);
	for(i=0;i<n;i++){
		for(h=0;h<3;h++){
			k2[i][h] = (v[i][h] + l1[i][h]*dt);
			l2[i][h] = a[i][h] * dt;
		}
	}

	for(i=0;i<n;i++){
		for(h=0;h<3;h++){
			k[i][h] = (k1[i][h] + k2[i][h])/2;
			k[i][h] = (k1[i][h] + k2[i][h])/2;
			x[i][h] += k[i][h];
			v[i][h] += l[i][h];
			}
	}
}


int  main (){

	double x[NMAX ][3];
	double v[NMAX ][3];
	double a[NMAX ][3];
	double m[NMAX];
	double p[NMAX];
	int n,i,j;
	double dt,tend;
	double e0;
	double tnow =0.0;
	double e;
	FILE *file1;	
	FILE *file2;	
	FILE *file3;	
	FILE *file4;		
	
	file1 = fopen("rung2_1.dat", "w");
	file2 = fopen("rung2_2.dat", "w");
	file3 = fopen("rung2_3.dat", "w");
	file4 = fopen("e_rung2.dat", "w");


	printf("input n,dt,tend,m,x,v");
	scanf("%d,%lf,%lf",n,dt,tend);
	for(i=0;i<n; i++){
		scanf("%lf",m[i]);
		for(j=0; j<3; j++){
			scanf("%lf",x[i][j]);
		}
		for(j=0; j<3; j++){
			scanf("%lf",v[i][j]);
		}
	}
	gravity(x,m,a,p,n);
	e0 = energy(v,m,p,n);

	/*printf("%lf\t",tnow):
	for(i=0;i<n;i++){
		for(j=0;j<3;j++){
		printf("%lf\t",x[i][j]);
	}*/
	while(tnow < tend){

	runge2(x,v,m,a,p,dt,n);
	e = energy(v,m,p,n);
	tnow += dt;

	/*printf("%lf\t",tnow):
	for(i=0;i<n;i++){
		for(j=0;j<3;j++){
		printf("%lf",x[i][j]);
	}
	printf("%lf\t",e);
	printf("%lf\t",(e-e0)/e0);*/
	
		
	fprintf(file1, " %lf %lf\n",x[0][0],x[0][1]);
    	fprintf(file2, " %lf %lf\n",x[1][0],x[1][1]);
	fprintf(file3, " %lf %lf\n",x[2][0],x[2][1]);
	fprintf(file4, " %lf %lf\n",tnow,(e-e0)/e0);
	
	fclose(file1);
	fclose(file2);
	fclose(file3);
	fclose(file4);	

	}
	

}
#include <stdio.h>
#include <stdlib.h>
#include<math.h>

static  const  int  NMAX = 128;
static  double  xtmp[NMAX][3];
static  double k[NMAX][3], k1[NMAX][3], k2[NMAX][3], k3[NMAX][3], k4[NMAX][3];
static  double l[NMAX][3], l1[NMAX][3], l2[NMAX][3], l3[NMAX][3], l4[NMAX][3];

void  gravity( double  (*x)[3],  double *m, double  (*a)[3],double *p, const int n){
	int i,j;
	for(i=0; i<n; i++){
		a[i][0]=a[i][1]=a[i][2]=p[i]=0.0;
		for(j=0;j<n;j++){
			if(i==j) continue;
			double dx = x[i][0] - x[j][0];
			double dy = x[i][1] - x[j][1];
			double dz = x[i][2] - x[j][2];
			double r2 = dx*dx + dy* +dz*dz;
			double rinv = 1.0 / sqrt(r2);
			double mjr3inv = m[j]*rinv*rinv*rinv;
			a[i][0] -= mjr3inv * dx;
			a[i][1] -= mjr3inv * dy;
			a[i][2] -= mjr3inv * dz;
			p[i]    -= m[j] * rinv;
		}
	}
}



double  energy( double  (*v)[3],  double *m, double *p, const  int n){

  double e = 0.0;
  int i;
  for(i=0;i<n;i++){
	e += p[i] + m[i]*(v[i][0]*v[i][0]+v[i][1]*v[i][1]+v[i][2]*v[i][2])/2;  
	}
	return e;
}


void  runge2( double  (*x)[3],  double  (*v)[3],  double *m, double  (*a)[3],
	double *p, const  double dt , const  int n){
	int i,h;
	gravity(x,m,a,p,n);
	for(i=0;i<n;i++){
		for(h=0;h<3;h++){
			k1[i][h] = v[i][h]*dt;
			l1[i][h] = a[i][h]*dt;
			xtmp[i][h] = x[i][h] + (k[i][h]/2);
		}
	}

	gravity(xtmp,m,a,p,n);
	for(i=0;i<n;i++){
		for(h=0;h<3;h++){
			k2[i][h] = (v[i][h] + l1[i][h]*dt);
			l2[i][h] = a[i][h] * dt;
		}
	}

	for(i=0;i<n;i++){
		for(h=0;h<3;h++){
			k[i][h] = (k1[i][h] + k2[i][h])/2;
			k[i][h] = (k1[i][h] + k2[i][h])/2;
			x[i][h] += k[i][h];
			v[i][h] += l[i][h];
			}
	}
}


int  main (){

	double x[NMAX ][3];
	double v[NMAX ][3];
	double a[NMAX ][3];
	double m[NMAX];
	double p[NMAX];
	int n,i,j;
	double dt,tend;
	double e0;
	double tnow =0.0;
	double e;
	FILE *file1;	
	FILE *file2;	
	FILE *file3;	
	FILE *file4;		
	
	file1 = fopen("rung2_1.dat", "w");
	file2 = fopen("rung2_2.dat", "w");
	file3 = fopen("rung2_3.dat", "w");
	file4 = fopen("e_rung2.dat", "w");


	printf("input n,dt,tend,m,x,v");
	scanf("%d,%lf,%lf",n,dt,tend);
	for(i=0;i<n; i++){
		scanf("%lf",m[i]);
		for(j=0; j<3; j++){
			scanf("%lf",x[i][j]);
		}
		for(j=0; j<3; j++){
			scanf("%lf",v[i][j]);
		}
	}
	gravity(x,m,a,p,n);
	e0 = energy(v,m,p,n);

	/*printf("%lf\t",tnow):
	for(i=0;i<n;i++){
		for(j=0;j<3;j++){
		printf("%lf\t",x[i][j]);
	}*/
	while(tnow < tend){

	runge2(x,v,m,a,p,dt,n);
	e = energy(v,m,p,n);
	tnow += dt;

	/*printf("%lf\t",tnow):
	for(i=0;i<n;i++){
		for(j=0;j<3;j++){
		printf("%lf",x[i][j]);
	}
	printf("%lf\t",e);
	printf("%lf\t",(e-e0)/e0);*/
	
		
	fprintf(file1, " %lf %lf\n",x[0][0],x[0][1]);
    	fprintf(file2, " %lf %lf\n",x[1][0],x[1][1]);
	fprintf(file3, " %lf %lf\n",x[2][0],x[2][1]);
	fprintf(file4, " %lf %lf\n",tnow,(e-e0)/e0);
	
	fclose(file1);
	fclose(file2);
	fclose(file3);
	fclose(file4);	

	}
	

}
現在、コメントはありません。
他の人よりも先にコメントしましょう。