// Maxwell_2D 를 통해서 얻은 평형 상태의 값을 불러서 계산함

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define ITER 50
#define T 0.005   // 단위 시간
#define N 10000   // 입자의 갯수
#define L 3000    // 박스의 너비
#define R 3       // 입자 지름

int main()
{
    int i,j,k,l;
    double vx[N],x[N];
    double vy[N],y[N];
    double pi,tmp1,tmp2;
    double initialE;
    double finalE;
    double v[N],vsquare[N];
    int low[N];                 //energy 0.9 - 1.1
    int high[N];                //energy 2.8 - 3.2
    double eng_low[N];          //energy 0.9 - 1.1
    double eng_high[N];         //energy 2.8 - 3.2
    int low_num, high_num;
    int low_limit, high_limit;
    int low_find, high_find;
    double eng;
    int iter;
    int elapse;
    FILE *logger;
    FILE *fp1;
    FILE *fp2;



    char outfile1[100];
    char outfile2[100];
    char outfile3[100];


    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);

    sprintf(outfile1,"log.dat");
    sprintf(outfile2,"deltaE_low_energy.dat");
    sprintf(outfile3,"deltaE_high_energy.dat");

    logger=fopen(outfile1,"w");
    fp1 = fopen(outfile2,"w");
    fp2 = fopen(outfile3,"w");


    for(iter=0;iter<ITER;iter++)
    {


        for(i=0; i<N; i++)
        {
            low[i]=0;
            high[i]=0;
            eng_low[i]=0.0;
            eng_high[i]=0.0;
        }

        low_num=0;
        high_num=0;
        low_limit=0;
        high_limit=0;
        low_find=0;
        high_find=0;
        elapse=0;


        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;
                }
            }
        }

        printf("\r낮은 에너지 높은 에너지 구별 중");
        for(i=0;i<N;i++)
        {
            eng = ((vx[i]*vx[i]+vy[i]*vy[i])/(double)2);
            if( eng < 1.1 &&  eng > 0.9 )                   // 낮은 에너지 구별, 번호를 low[]에 순서대로 저장
            {
                low[low_num]=i;
                low_num++;
            }

            if( eng < 3.2 &&  eng > 2.8 )                   // 높은 에너지 구별, 번호를 high[]에 순서대로 저장
            {
                high[high_num]=i;
                high_num++;
            }
        }

        low_limit=low_num*0.5;
        high_limit=high_num*0.5;


        system("cls");
        fprintf(logger,"# of low : %d \t # of high %d\n",low_num,high_num);
        printf("#     of low : %d \t #     of high %d\n",low_num,high_num);
        printf("Limit of low : %d \t Limit of high %d\n",low_limit,high_limit);

        while(high_find < high_limit && low_find < low_limit )
        {

            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]));                           // pi값 제시
                        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;                 // 나중 에너지 계산

                        for(l=0;l<low_num;l++)                                      // 낮은 에너지의 입자가 충돌했는지 확인하기
                        {
                            if(low[l]==j)                                           // 낮은 에너지의 입자인지 확인
                                if(eng_low[l]==0)                                   // 첫번째 충돌인지 확인
                                {
                                    eng_low[l]=finalE-initialE;                     // 에너지 기록, finalE initialE 는 j가 기준
                                    low_find++;
                                    printf("L %d %f\n",low_find,eng_low[l]);
                                    fprintf(fp1, "%f\n",eng_low[l]);                // 출력
                                    continue;
                                }

                            if(low[l]==k)
                                if(eng_low[l]==0)
                                {
                                    eng_low[l]=-finalE+initialE;                    // k의 경우 deltaE 가 반대
                                    low_find++;
                                    printf("L %d %f\n",low_find,eng_low[l]);
                                    fprintf(fp1, "%f\n",eng_low[l]);
                                    continue;
                                }
                        }

                        for(l=0;l<high_num;l++)                                      // 높은 에너지 갯수만큼 다 확인해보기 (위와 같음)
                        {
                            if(high[l]==j)
                                if(eng_high[l]==0)
                                {
                                    eng_high[l]=finalE-initialE;
                                    high_find++;
                                    printf("H %d %f\n",high_find,eng_high[l]);
                                    fprintf(fp2, "%f\n",eng_high[l]);
                                    continue;
                                }

                            if(high[l]==k)
                                if(eng_high[l]==0)
                                {
                                    eng_high[l]=-finalE+initialE;
                                    high_find++;
                                    printf("H %d %f\n",high_find,eng_high[l]);
                                    fprintf(fp2, "%f\n",eng_high[l]);
                                    continue;
                                }
                        }
                    }
                }
            }
            elapse++;
        }

        fprintf(logger,"elapse : %d\nlow_find : %d\nhigh_find : %d\n",elapse,low_find,high_find);   //로그 파일에 기록

    }

    return 0;

}