// Maxwell_2D 를 통해서 얻은 평형 상태의 값을 불러서 계산함

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define ITINERATION 50000 // 반복 횟수
#define T 0.005           // 단위 시간
#define N 10000           // 입자의 갯수
#define L 3000            // 박스의 너비
#define R 3               // 입자 지름

int main()
{
    int i,j,k, num=0;
    double vx[N],x[N];
    double vy[N],y[N];
    double pi,tmp1,tmp2;
    double initialE;
    double finalE;
    double v[N],vsquare[N];
    double deltaE[100000]={};

    srand(time(NULL));

    FILE *vxx = fopen("vx.dat","r");
    FILE *vyy = fopen("vy.dat","r");
    FILE *xx = fopen("x.dat","r");
    FILE *yy = fopen("y.dat","r");


    fflush(stdin);
    for(i=0;i<N;i++)
    {
        fscanf(vxx,"%lf\n",&vx[i]);
        fscanf(vyy,"%lf\n",&vy[i]);
        fscanf(xx,"%lf\n",&x[i]);
        fscanf(yy,"%lf\n",&y[i]);
    }

    fclose(vxx);
    fclose(vyy);
    fclose(xx);
    fclose(yy);


    printf("\r위치 조정중");
    for(i=0;i<N-1;i++)              // 위치 중복내 확인
    {
        for(j=i+1;j<N;j++)
        {
            if(sqrt((x[i]-x[j])*(x[i]-x[j])+(y[i]-y[j])*(y[i]-y[j]))<R)
            {
                x[j]=(double)rand()/RAND_MAX*L;
                y[j]=(double)rand()/RAND_MAX*L;
                i=0;
                j=0;
            }
        }
    }

    system("cls");

    for(i=0;i<ITINERATION;i++)
    {

        for(j=0;j<N;j++)
        {
            x[j]+=vx[j]*T;            //x방향으로 T시간만큼 움직임
            y[j]+=vy[j]*T;            //y방향으로 T시간만큼 움직임

            if(x[j]>L)
                vx[j]=-fabs(vx[j]);   //0,L 경계면에서 되튐
            if(x[j]<0)
                vx[j]=fabs(vx[j]);

            if(y[j]>L)
                vy[j]=-fabs(vy[j]);   //0,L 경계면에서 되튐
            if(y[j]<0)
                vy[j]=fabs(vy[j]);
        }

        for(j=0;j<N-1;j++)
        {
            for(k=j+1;k<N;k++)
            {
                if(sqrt((x[k]-x[j])*(x[k]-x[j])+(y[k]-y[j])*(y[k]-y[j]))<R)     // 서로 R보다 작을경우(충돌)
                {
                    pi=atan((y[k]-y[j])/(x[k]-x[j]));                           // 충돌각 제시
                    initialE=(vx[j]*vx[j]+vy[j]*vy[j])/(double)2;

                    tmp1=vx[j]*cos(pi)+vy[j]*sin(pi);
                    tmp2=vx[k]*cos(pi)+vy[k]*sin(pi);
                    vx[j]=vx[j]+tmp2*cos(pi)-tmp1*cos(pi);
                    vy[j]=vy[j]+tmp2*sin(pi)-tmp1*sin(pi);
                    vx[k]=vx[k]+tmp1*cos(pi)-tmp2*cos(pi);
                    vy[k]=vy[k]+tmp1*sin(pi)-tmp2*sin(pi);                      // 탄성충돌해서 바뀌는 속도 계산

                    finalE=(vx[j]*vx[j]+vy[j]*vy[j])/(double)2;


                    //printf("%d %d %d\n",num,j,k);                             // 충돌한 공의 number 모니터 출력

                    deltaE[num]=fabs(finalE-initialE);
                    num++;
                }
            }
        }
        printf("\r%.2f%% 진행 ",(double)i/ITINERATION*100);
    }


    FILE *fp;
    fp=fopen("deltaE.dat","w");
    for(i=0;i<num;i++)
        fprintf(fp,"%f\n",deltaE[i]);                //deltaE값 파일에 출력
    fclose(fp);

    return 0;

}