Harigami
ログイン
まんだげり タイトルなし
No License 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);  
  }
#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);  
  }
現在、コメントはありません。
最初のコメンターになりませんか?