#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fstream.h>
#include <iostream.h>

/*-------------------構造點的類型-------------------*/
struct Point{
        double x;
        double y;
};

/*--------------------全局變量----------------------*/
double G;               //萬有引力常量  
int N;                  //天體的數目
Point* P;               //每個天體的位置
Point* v;               //每個天體的速度
Point* f;               //每個天體受到的力
double* m;              //每個天體的質量
Point size;             //空間的大小
float C;                //一常數
double dt;              //很小的時間周期

/*---------------初始化各天體的情況-----------------*/
void Initial(char* filename){
	ifstream f1(filename);
	if(!f1){
		cout<<"open file fail!"<<endl;
		return;
	}
	f1>>N;         //從文件讀入天體的數目      
	P = new Point[N];
	v = new Point[N];
	f = new Point[N];
	m = new double[N];
	f1>>G>>size.x>>size.y>>dt>>C;      //從文件讀入G，size，dt
	for(int i=0;i<N;i++){
		f1>>m[i]>>P[i].x>>P[i].y>>v[i].x>>v[i].y;     //從文件讀入各天體的質量，位置和速度
	}
	for(int j=0;j<N;j++) { f[j].x = 0.0; f[j].y = 0.0; }   //初始化每個天體的受力為0
	f1.close();
}

/*-----------計算dt時間內每個天體受到的力-----------*/
void calculateForces(int i){
	double distance,magnitude;
	Point direction;
	for(int j=0;j<N;j++){
		if(j!=i){
			distance = (float)sqrt((P[i].x-P[j].x)*(P[i].x-P[j].x)+(P[i].y-P[j].y)*(P[i].y-P[j].y));    //計算兩點間距離
			if(distance!=0){
				magnitude = (G*m[i]*m[j])/(distance*distance+C*C);   //計算萬有引力
				direction.x = P[j].x-P[i].x;                         //記錄力的x方向
				direction.y = P[j].y-P[i].y;                         //記錄力的y方向
				f[i].x = f[i].x+magnitude*direction.x/distance;      //把收到的所有萬有引力累加
				f[i].y = f[i].y+magnitude*direction.y/distance;
			}
		}
	}
}

/*----------計算dt時間后，天體的位置和速度----------*/
void moveBodies(int i){
	Point dv;
	Point dp;
	Point tempPoint;
	dv.x = f[i].x/m[i]*dt;       dv.y = f[i].y/m[i]*dt;         //計算x，y方向的dv＝a*dt,a=F/m
	dp.x = (v[i].x+dv.x/2)*dt;   dp.y = (v[i].y+dv.y/2)*dt;     //計算x，y偏移量dp＝(v+dv/2)dt
	tempPoint.x=P[i].x+dp.x;     tempPoint.y = P[i].y+dp.y;     //對位置的偏移量進行累加
	if(tempPoint.x<0||tempPoint.x>size.x||tempPoint.y<0||tempPoint.y>size.y){  //判斷下一位置是否越界
		P[i].x = P[i].x;     //如果越界，天體保留原來位置，速度方向取反
		P[i].y = P[i].y;
		v[i].x = -v[i].x;
		v[i].y = -v[i].y;
	}else{                   //記錄新的位置和方向
		P[i].x = tempPoint.x;
		P[i].y = tempPoint.y;
		v[i].x = v[i].x + dv.x;
		v[i].y = v[i].y + dv.y;
	}
	f[i].x = 0.0; f[i].y = 0.0;
}

/*-------------------輸出各天體的情況---------------*/
void End(char* filename){
	ofstream f2(filename);
	if(!f2){
		cout<<"open file fail!"<<endl;
		return;
	}
	f2<<N<<endl<<G<<endl<<size.x<<" "<<size.y<<endl<<dt<<endl<<C<<endl;        //輸出n，G，size，dt到結果文件
	for(int i=0;i<N;i++){
		f2<<m[i]<<" "<<P[i].x<<" "<<P[i].y<<" "<<v[i].x<<" "<<v[i].y<<endl;    //輸出各天體的新情況到結果文件
	}
	delete[] P;
	delete[] v;
	delete[] f;
	delete[] m;
	f2.close();
}




/*--------------------------主函數------------------*/
int main(int argc,char* argv[]){
	int myid,numprocs;
    double starttime,endtime;               //記錄運行的開始時間和結束時間

    MPI_Init(&argc,&argv);                  //啟動MPI計算   
    MPI_Comm_size(MPI_COMM_WORLD,&numprocs);//確定進程數
    MPI_Comm_rank(MPI_COMM_WORLD,&myid);    //確定自己的進程標識
    starttime=MPI_Wtime();                  //記錄當前時間
    Initial("sample_input.in");             //讀取數據
	double* PP = new double[2*N];           //記錄天體的x，y坐標
	double* vv = new double[2*N];           //記錄天體的x，y方向
	int d=N/numprocs;                       //記錄每個進程分配到的天體的數目
	for(int time=0; time<1000;time+=1){     //進行循環計算1000次
		for(int i=myid*d;i<(myid+1)*d;i++){ //每個進程計算d個天體
			calculateForces(i);             //計算第i個天體的力
			moveBodies(i);                  //計算第i個天體新的位置和力
			PP[2*i]=P[i].x;                 //記錄天體的坐標，方向
			PP[2*i+1]=P[i].y;
			vv[2*i]=v[i].x;
			vv[2*i+1]=v[i].y;
		}
		if(myid==(numprocs-1)){             //如果是最后一個進程，把剩下的沒分配的天體進行計算
			for(int j=(myid+1)*d;j<N;j++){
				calculateForces(j);
				moveBodies(j);
				PP[2*j]=P[j].x;
				PP[2*j+1]=P[j].y;
				vv[2*j]=v[j].x;
				vv[2*j+1]=v[j].y;
			}
		}
		MPI_Barrier(MPI_COMM_WORLD);        //進行各進程同步
		for(int j=0;j<numprocs;j++){        //各進程依次廣播自己計算的位置
			MPI_Barrier(MPI_COMM_WORLD);    
			if (j!=(numprocs-1))	MPI_Bcast(&PP[j*2*d],2*d,MPI_DOUBLE,j,MPI_COMM_WORLD);    
			else MPI_Bcast(&PP[j*2*d],(N-j*d)*2,MPI_DOUBLE,j,MPI_COMM_WORLD);
			MPI_Barrier(MPI_COMM_WORLD);    //廣播完各進程同步
			for(int k=j*d;k<(j+1)*d;k++){   //收到廣播信息后，其他進程更新自己的值
				P[k].x=PP[2*k];
				P[k].y=PP[2*k+1];
			}
			MPI_Barrier(MPI_COMM_WORLD);    //更新完成后同步
		}
	}
	MPI_Barrier(MPI_COMM_WORLD);            //計算完1000后同步
	for(int j=0;j<numprocs;j++){            //各進程依次廣播自己的速度
		MPI_Barrier(MPI_COMM_WORLD);        //發送速度前同步
		if (j!=(numprocs-1))	MPI_Bcast(&vv[j*2*d],2*d,MPI_DOUBLE,j,MPI_COMM_WORLD);   
		else MPI_Bcast(&vv[j*2*d],(N-j*d)*2,MPI_DOUBLE,j,MPI_COMM_WORLD);
		MPI_Barrier(MPI_COMM_WORLD);        //廣播完同步
		for(int k=j*d;k<(j+1)*d;k++){       //收到廣播信息后，其他進程更新自己的值
			v[k].x=vv[2*k];
			v[k].y=vv[2*k+1];
		}
	}
	MPI_Barrier(MPI_COMM_WORLD);            //計算完成后同步
	if(myid==0)     End("result1000.data"); //如果是0進程，輸出運行結果到文件中
	endtime=MPI_Wtime();                    //計算每個進程整個程序的運行時間
	if(myid==0) cout<<"總的運行時間："<<endtime-starttime<<"s"<<endl;     //輸出每個進程的運行時間 
	delete[] PP;
	delete[] vv; 
	MPI_Finalize();                         //結束MPI運算 
	return 0;
}