// PROGRAM perc_cluster
// ハマースリー，リース，アレクサンドロビツらのアルゴリズムによる
// クラスターの生成

#include "TrueBASIC.h"
#include "csgraphics.h"

int main();
void initial(double *p, int *L, int status[]);
void initialize_arrays(int xs[], int ys[]);
void grow(double p, int L, int *N, int xs[], int ys[], int status[]);
void mass_dist(int L, int N, int xs[], int ys[]);

int main()
{
    int L, N, xs[10001], ys[10001], status[3];
    double p;

    GWopen(0);

    initial(&p, &L, status);
    initialize_arrays(xs, ys);
    grow(p, L, &N, xs, ys, status);
    mass_dist(L, N, xs, ys);

    GWquit();
    return 0;
}

void initial(double *p, int *L, int status[])
{
    int x, y;
    float xwin, ywin, sz;
    char Stmp1_[_LBUFF_];

    rnd(-1);
    printf("value of L (odd) = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", L);
    printf("site occupation probability = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", p);
    compute_aspect_ratio((float)(*L), &xwin, &ywin);
    GWindow(0.0f, 0.0f, xwin, ywin);
    sz = ywin/30;
    GWsettxt(sz, 0.0f, 5, -1, -1, NULL);
    GWsetpen(BLACK, -1, -1, -1);
    GWrect(0.5f, 0.5f, (float)((*L)+0.5), (float)((*L)+0.5));
    for(y = 1; y <= (*L); ++y)
        for(x = 1; x <= (*L); ++x)
            GWsetpxl((float)(x), (float)(y), BLACK);
    GWsetpen(RED, -1, -1, 4);
    status[1] = GWcmbmrk(2, 6, 1.0f, RED, -1, 4);     // 占有された格子点
    GWputcmb(status[1], 0.5f+sz/2, (*L)+0.5f + sz, sz, sz, 0);
    GWputtxt(0.5f+sz*1.5f, (*L)+0.5f + sz, "occupied site");
    GWsetpen(BLUE, -1, -1, 4);
    status[2] = GWcmbmrk(3, 6, 1.0f, BLUE, -1, 4);    // 周辺の点
    GWputcmb(status[2], (*L)/3.0f+0.5f+sz/2, (*L)+0.5f+sz, sz, sz, 0);
    GWputtxt((*L)/3.0f+0.5f+sz*1.5f, (*L)+0.5f+sz, "perimeter site");
    GWsetpen(BLACK, -1, -1, 4);
    status[0] = GWcmbmrk(4, 6, 1.0f, BLACK, -1, 4);   // 試されたが占有されなかった格子点
    GWputcmb(status[0], (*L)/3.0f*2+0.5f+sz/2, (*L)+0.5f + sz, sz, sz, 0);
    GWputtxt((*L)/3.0f*2+0.5f+sz*1.5f, (*L)+0.5f+sz, "tested site");
}

void initialize_arrays(int xs[], int ys[])
{
    int Itmp1_;

    for(Itmp1_ = 0; Itmp1_ <= 10000; ++Itmp1_)
        xs[Itmp1_] = 0;
    for(Itmp1_ = 0; Itmp1_ <= 10000; ++Itmp1_)
        ys[Itmp1_] = 0;
}

void grow(double p, int L, int *N, int xs[], int ys[], int status[])
// 1つのパーコレーションクラスターを生成
{
    int i, iper, nn, nper, perx[25001], pery[25001], x, xnew, xseed, y, ynew, 
        yseed;
    int nx[5], ny[5];                       // 隣接格子点への方向ベクトル
    int site[132][132];
    int DP_ = 0;
    int DATA_[] = {
        1,0,-1,0,0,1,0,-1
    };

    // 境界の点の設定
    for(i = 1; i <= L; ++i)
    {
        site[0][i] = -1;
        site[L + 1][i] = -1;
        site[i][L + 1] = -1;
        site[i][0] = -1;
    }
    // 格子の中央に種を置く
    xseed = L/2 + 1;
    yseed = xseed;
    site[xseed][yseed] = 1;                 // 種の格子点
    xs[1] = xseed;
    ys[1] = yseed;
    (*N) = 1;                               // クラスター内の格子点数
    GWputcmb(status[1], (float)(xseed), (float)(yseed), 1.0f, 1.0f, 0);
    for(i = 1; i <= 4; ++i)
    {
        // nx,ny: 新しい周辺の点に対する方向ベクトル
        nx[i] = DATA_[DP_++];
        ny[i] = DATA_[DP_++];
        // perx,pery: 種の周辺の点の位置
        perx[i] = xseed + nx[i];
        pery[i] = yseed + ny[i];
        // 周辺の点には2を割り当てる
        site[perx[i]][pery[i]] = 2;         // 周辺の点の表に入れられた格子点
        GWputcmb(status[2], (float)(perx[i]), (float)(pery[i]), 1.0f, 1.0f, 0);
    }
    nper = 4;                               // 周辺の点の数の初期値
    do
    {
        // ランダムに周辺の点を選ぶ
        iper = (int)(rnd(0)*nper) + 1;
        x = perx[iper];                     // 周辺の点の座標
        y = pery[iper];
        // 配列中の最後の周辺の点で新しく選ばれた格子点を置き換える
        perx[iper] = perx[nper];
        pery[iper] = pery[nper];
        --nper;
        if(rnd(0) < p)                      // 占有された格子点
        {
            site[x][y] = 1;
            ++(*N);
            xs[(*N)] = x;                   // 占有された格子点を保存
            ys[(*N)] = y;
            GWputcmb(status[1], (float)(x), (float)(y), 1.0f, 1.0f, 0);
            for(nn = 1; nn <= 4; ++nn)   // 新しい周辺の点を探す
            {
                xnew = x + nx[nn];
                ynew = y + ny[nn];
                if(site[xnew][ynew] == 0)
                {
                    ++nper;
                    perx[nper] = xnew;
                    pery[nper] = ynew;
                    // 周辺の点に加える
                    site[xnew][ynew] = 2;
                    GWputcmb(status[2], (float)(xnew), (float)(ynew), 1.0f, 1.0f, 0);
                }                           // rnd >= p ならば，格子点は占有されない
            }
        }
        else
        {
            site[x][y] = -1;
            GWputcmb(status[0], (float)(x), (float)(y), 1.0f, 1.0f, 0);
        }
    } while(!(nper < 1));                   // 全ての周辺の点が試される
    GWsetbrs(BLACK, 0, -1);
    GWrect(0.5f, 0.5f, (float)(L+0.5), (float)(L+0.5));
                                            // 四角形の再描画を行う
}

void mass_dist(int L, int N, int xs[], int ys[])
{
    int i, mass[10001], masstotal = 0, r, rprint;
    double dx, dy, xcm = 0, ycm = 0;
    for(i = 1; i <= 10000; ++i)
        mass[i] = 0;

    for(i = 1; i <= N; ++i)
    {
        xcm += xs[i];                       // 重心を計算
        ycm += ys[i];
    }
    xcm /= N;
    ycm /= N;
    for(i = 1; i <= N; ++i)
    {
        dx = xs[i] - xcm;
        dy = ys[i] - ycm;
        r = (int)(sqrt(dx*dx + dy*dy));
        // 重心からの距離
        // mass(r) = 重心からの距離 r の格子点数
        if(r > 1)
            ++(mass[r]);
    }
    rprint = 2;
    printf("           r            m        ln(r)        ln(m) \n");
    for(r = 2; r <= L/2; ++r)
    {
        masstotal += mass[r];
        if(r == rprint)
        {
            printf("%12d %12d %12f %12f\n", r, masstotal, log((double)r), 
                   log((double)masstotal));
            rprint *= 2;                    // r については対数目盛を使う
        }
    }
}