Harigami
タイトルなし
C
#include <stdio.h>
#include <stdlib.h>
#include<math.h>
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*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,j;
int nmax = 128;
double xtmp[nmax][3];
double k[nmax][3], k1[nmax][3], k2[nmax][3], k3[nmax][3], k4[nmax][3];
double l[nmax][3], l1[nmax][3], l2[nmax][3], l3[nmax][3], l4[nmax][3];
gravity(x,m,a,p,n);
for(i=0;i<n;i++){
for(j=0;j<n;j++){
k1[i][j] = v[i][j]*dt;
l1[i][j] = a[i][j]*dt;
xtmp[i][j] = x[i][j] + (k1[i][j]/2);
}
}
gravity(xtmp,m,a,p,n);
for(i=0;i<n;i++){
for(j=0;j<n;j++){
k2[i][j] = (v[i][j] + l1[i][j]*dt);
l2[i][j] = a[i][j] * dt;
}
}
for(i=0;i<n;i++){
for(j=0;j<n;j++){
k3[i][j] = (v[i][j] + l2[i][j]*dt);
l3[i][j] = a[i][j] * dt;
}
}
for(j=0;j<n;j++){
k4[i][j] = (v[i][j] + l3[i][j]*dt);
l4[i][j] = a[i][j] * dt;
k[i][j] = (k1[i][j] + 2*k2[i][j]+2*k3[i][j]+k4[i][j])/2;
l[i][j] = (l1[i][j] + 2*l2[i][j]+2*l3[i][j]+l4[i][j])/2;
x[i][j] += k[i][j];
v[i][j] += l[i][j];
}
}
int main (){
int nmax =128;
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("rung4_1.dat", "w");
file2 = fopen("rung4_2.dat", "w");
file3 = fopen("rung4_3.dat", "w");
file4 = fopen("e_rung4.dat", "w");
printf("input n,dt,tend,m,x,v\n");
scanf("%d %lf %lf",&n,&dt,&tend);
for(i=0;i<3;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);
}
タイトルなし
タイトルなし