// PROGRAM R2
// 可能なすべての原点のとり方についての平均から
// ある1つの粒子の平均２乗変位を計算する

#include "TrueBASIC.h"
int main();
void initial(double *Lx, double *Ly, double R2cum[], int *ndiff);
void read_data(double x[], double y[], int *ndata);
void displacements(double x[], double y[], double Lx, double Ly, 
                   double R2cum[], int norm[], int ndata, int ndiff);
void normalization(int ndiff, double R2cum[], int norm[]);
double separation(double ds, double L);

int main()
{
    int ndata, ndiff, norm[21], Itmp1_;
    double Lx, Ly, R2cum[21], x[1001], y[1001];
    for(Itmp1_ = 0; Itmp1_ < 20; ++Itmp1_)
       norm[Itmp1_] = 0;

    initial(&Lx, &Ly, R2cum, &ndiff);
    read_data(x, y, &ndata);
    displacements(x, y, Lx, Ly, R2cum, norm, ndata, ndiff);
    normalization(ndiff, R2cum, norm);

    return 0;
}

void initial(double *Lx, double *Ly, double R2cum[], int *ndiff)
{
    int idiff;

    (*Lx) = 5;
    (*Ly) = 5;
    (*ndiff) = 10;                          // 最大の時間差
    for(idiff = 1; idiff <= (*ndiff); ++idiff)
        R2cum[idiff] = 0;
}

void read_data(double x[], double y[], int *ndata)
// 一定の時間間隔で保存された粒子位置のファイルを読み込む
{
    int t;
    char Stmp1_[_LBUFF_];
    FILE *FP1_;

    FP1_ = fopen("xy.dat", "r");
    t = 0;
    while(ungetc(getc(FP1_),FP1_) != EOF)
    {
        ++t;
        fgets(Stmp1_, _LBUFF_, FP1_);
        sscanf(Stmp1_, " %lg, %lg", &(x[t]), &(y[t]));
    }
    fclose(FP1_);
    (*ndata) = t;                           // データ点の数
}

void displacements(double x[], double y[], double Lx, double Ly, 
                   double R2cum[], int norm[], int ndata, int ndiff)
{
    int i, idiff;
    double dx, dy;

    for(idiff = 1; idiff <= ndiff; ++idiff)
    {
        for(i = 1; i <= ndata - idiff; ++i)
        {
            dx = separation(x[i+idiff] - x[i], Lx);
            dy = separation(y[i+idiff] - y[i], Ly);
            R2cum[idiff] += dx*dx + dy*dy;
            ++(norm[idiff]);
        }
    }
}

void normalization(int ndiff, double R2cum[], int norm[])
{
    int idiff;
    double R2bar;

    printf("time difference      <R2>\n");
    printf("\n");
    for(idiff = 1; idiff <= ndiff; ++idiff)
    {
        if(R2cum[idiff] > 0)
        {
            R2bar = R2cum[idiff]/norm[idiff];
            printf("%12d %12f\n", idiff, R2bar);
        }
    }
}

double separation(double ds, double L)      // プログラム md と同じ関数
{
    double separation_;

    if(ds > 0.5*L)
        separation_ = ds - L;
    else if(ds < -0.5*L)
        separation_ = ds + L;
    else
        separation_ = ds;
    return separation_;
}