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

#define G 6.672e-11
#define M 2.058915e41
#define m 1.989e30
#define E 2.998e8
#define Pi 3.142593
#define dt 60*60*24.0*365
#define Ring 10
#define S 20

typedef struct star{
  double Fx;
  double Fy;
  double Ax;
  double Ay;
  double Vx;
  double Vy;
  double x;
  double y;
}Star;

main(){
  Star star[200];
  int i,j,k,a;
  double r,v,F,tmp;

  r=0.0;a=0;

  FILE *pp,*fp;

  pp=popen("gnuplot -geometry +100+100","w");
  if(pp==NULL){
    perror("Cannot execute gnuplot");
    exit(1);
  }
  fprintf(pp,"set yrange [-5:5]\n");
  fprintf(pp,"set xrange [-5:5]\n");

  fp=fopen("tmp.dat","w+");
  if(fp==NULL){
    perror("Cannot open temorary file");
    exit(1);
  }
  
    /*初期設定*/
  for(i=0;i<Ring;i++){
    r=r+E*60*60*24*365*0.5;
    v=sqrt((G*M)/(r));
    for(j=0;j<S;j++){
      tmp=2*j*Pi/(S*1.0);
      star[a].x=r*cos(tmp);
      star[a].y=r*sin(tmp);
      star[a].Vx=-v*sin(tmp);
      star[a].Vy=v*cos(tmp);
      a++;
    }
  }

  for(k=0;k<1000;k++){
    /*力和の計算*/
    for(i=0;i<Ring*S;i++){
      tmp=star[i].x*star[i].x+star[i].y*star[i].y;
      F=G*M*m/(tmp);
      star[i].Fx=F*(-star[i].x/sqrt(tmp));
      star[i].Fy=F*(-star[i].y/sqrt(tmp));
      for(j=0;j<Ring*S;j++){
	if(i!=j){
          double tmp_x,tmp_y;
          tmp_x=star[i].x-star[j].x;
          tmp_y=star[i].y-star[j].y;
          tmp=tmp_x*tmp_x+tmp_y*tmp_y;
          F=G*m*m/(tmp);
          star[i].Fx=star[i].Fx+F*(-tmp_x/sqrt(tmp));
          star[i].Fy=star[i].Fy+F*(-tmp_y/sqrt(tmp));
       	}
      }
    }
    /*位置の計算*/
    for(i=0;i<Ring*S;i++){
      star[i].Ax=star[i].Fx/m;
      star[i].Ay=star[i].Fy/m;
      star[i].Vx=star[i].Vx+star[i].Ax*dt;
      star[i].Vy=star[i].Vy+star[i].Ay*dt;
      star[i].x=star[i].x+star[i].Vx*dt;
      star[i].y=star[i].y+star[i].Vy*dt;
    }
    /*gnuplotに出力*/
    if(k%10==0){
      double x,y;
      fseek(fp,0,SEEK_SET);
      for(i=0;i<Ring*S;i++){
	x=star[i].x/(r/5.0);
	y=star[i].y/(r/5.0);
	fprintf(fp,"%10.6f %10.6f\n",x,y);
      }
      fflush(fp);
      fprintf(pp,"plot 'tmp.dat' using 1:2 title '' with points ps 3\n");
      fflush(pp);
      sleep(1);
    }
  }
}