#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            // 박스의 너비 (2차원)
#define R 3               // 입자 지름

int main()
{
    int i,j,k;
    double vx[N],x[N];
    double vy[N],y[N];
    double pi,tmp1,tmp2;
    double v[N];

    srand(time(NULL));

    printf("초기값 대입");
    for(i=0;i<N;i++)
    {
        vx[i]=2-(double)rand()/RAND_MAX*4;          //-2<vx<2
        vy[i]=2-(double)rand()/RAND_MAX*4;          //-2<vy<2
        x[i]=(double)rand()/RAND_MAX*L;             //0<x<L
        y[i]=(double)rand()/RAND_MAX*L;             //0<y<L
    }

    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)     // 두 입자간 사이가 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]));           //충돌각 제시

                    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);      // 탄성충돌해서 바뀌는 속도 계산

                    //printf("충돌! %d %d\n",j,k);              // 충돌한 공의 number 모니터 출력
                }
            }
        }
        //printf("\r%.2f%% 진행 ",(double)i/ITINERATION*100);     // 진행 현황 출력
    }

    for(i=0;i<N;i++)
    {
        v[i]=sqrt(vx[i]*vx[i]+vy[i]*vy[i]);     //v값 계산
    }


    FILE *fp;
    fp=fopen("vx.dat","w");
    for(i=0;i<N;i++)
        fprintf(fp,"%f\n",vx[i]);          //vx값 파일에 출력
    fclose(fp);

    fp=fopen("vy.dat","w");
    for(i=0;i<N;i++)
        fprintf(fp,"%f\n",vy[i]);          //vy값 파일에 출력
    fclose(fp);

    fp=fopen("v.dat","w");
    for(i=0;i<N;i++)
        fprintf(fp,"%f\n",v[i]);           //v값 파일에 출력
    fclose(fp);

    fp=fopen("x.dat","w");
    for(i=0;i<N;i++)
        fprintf(fp,"%f\n",x[i]);           //x값 파일에 출력
    fclose(fp);

    fp=fopen("y.dat","w");
    for(i=0;i<N;i++)
        fprintf(fp,"%f\n",y[i]);           //y값 파일에 출력
    fclose(fp);

}