#include <stdio.h>
#include <math.h>
#include "mpi.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*1000.0
#define Ring 4
#define S 4


main(int argc,char* argv[]){

  int i,j,a;
  double r,v,F,tmp;
  int p,rank,nn,is,ie;

  double Fx[Ring*S];
  double Fy[Ring*S];
  double Ax[Ring*S];
  double Ay[Ring*S];
  double Vx[Ring*S];
  double Vy[Ring*S];
  double x[Ring*S];
  double y[Ring*S];

  MPI_Init(&argc,&argv);
  MPI_Comm_size(MPI_COMM_WORLD,&p);
  MPI_Comm_rank(MPI_COMM_WORLD,&rank);

  /*
  if(rank==0)
    scanf("%d\n",&S);
  MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);
  */
  nn=Ring*S/p;
  is=nn*rank;
  ie=is+nn;

  r=0.0;a=0;
  
    /*初期設定*/
  for(i=0;i<Ring;i++){
    r=r+E*60*60*24*365*10000*0.5;
    v=sqrt((G*M)/(r));
    for(j=0;j<S;j++){
      tmp=2*j*Pi/(S*1.0);
      x[a]=r*cos(tmp);
      y[a]=r*sin(tmp);
      Vx[a]=-v*sin(tmp);
      Vy[a]=v*cos(tmp);
      a++;
    }
  }

    /*力和の計算*/
    for(i=is;i<ie;i++){
tmp=x[i]*x[i]+y[i]*y[i];
      F=G*M*m/(tmp);
      Fx[i]=F*(-x[i]/sqrt(tmp));
      Fy[i]=F*(-y[i]/sqrt(tmp));
      for(j=0;j<Ring*S;j++){
	if(i!=j){
          double tmp_x,tmp_y;
          tmp_x=x[i]-x[j];
          tmp_y=y[i]-y[j];
          tmp=tmp_x*tmp_x+tmp_y*tmp_y;
          F=G*m*m/(tmp);
          Fx[i]=Fx[i]+F*(-tmp_x/sqrt(tmp));
          Fy[i]=Fy[i]+F*(-tmp_y/sqrt(tmp));
       	}
      }      
    }
    /*位置の計算*/
    for(i=is;i<ie;i++){
      Ax[i]=Fx[i]/m;
      Ay[i]=Fy[i]/m;
      Vx[i]=Vx[i]+Ax[i]*dt;
      Vy[i]=Vy[i]+Ay[i]*dt;
      x[i]=x[i]+Vx[i]*dt;
      y[i]=y[i]+Vy[i]*dt;
    }
   
    
    for(i=0;i<p;i++){
      MPI_Bcast(&x[i*nn],nn,MPI_DOUBLE,i,MPI_COMM_WORLD);
      MPI_Bcast(&y[i*nn],nn,MPI_DOUBLE,i,MPI_COMM_WORLD);
    }

if(rank==0)
for(i=0;i<Ring*S;i++)printf("%f %f\n",x[i],y[i]);
    MPI_Finalize();
}