// PROGRAM eden

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

int main();
void initial(int *L);
void grow(int L, int mass[]);
void mass_dist(int L, int mass[]);

int main()
{
    int L, mass[1001], Itmp1_;
    for(Itmp1_ = 0; Itmp1_ <= 1000; ++Itmp1_)
       mass[Itmp1_] = 0;

    GWopen(0);

    initial(&L);
    grow(L, mass);
    mass_dist(L, mass);

    GWquit();
    return 0;
}

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

    rnd(-1);
    printf("value of L (odd) = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", L);
    GWvport(0.01f*GWaspect(1), 0.01f, 0.7f*GWaspect(1), 
            0.99f);
    compute_aspect_ratio((float)(*L), &xwin, &ywin);
    GWindow(0.0f, 0.0f, xwin, ywin);
    GWsetpen(RED, -1, -1, 4);
    GWcmbmrk(1, 6, 1.0f, RED, -1, 4);
    GWsetpen(BLUE, -1, -1, 4);
    GWcmbmrk(2, 6, 1.0f, BLUE, -1, 4);
    GWsetpen(BLACK, -1, -1, 4);
    GWsetbrs(BLACK, 0, -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);    // チェックされていない占有されていない点
}

void grow(int L, int mass[])
// 単一のパーコレション・クラスターの生成
{
    int i, iper, nn, nper, perx[5001], pery[5001], r, x, xnew, xseed, y, ynew, 
        yseed;
    int nx[5], ny[5];                   // 最隣接点へのベクトル
    double dx, dy, site[202][202];
    char boundary[_LBUFF_];
    int DP_ = 0;
    int DATA_[] = {
        1,0,-1,0,0,1,0,-1
    };

    strcpy(boundary, "");
    // 境界の設定
    for(i = 1; i <= L; i += 1)
    {
        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;                 // 種の格子点
    GWputcmb(1, (float)(xseed), (float)(yseed), 1.0f, 1.0f, 0);
    for(i = 1; i <= 4; i += 1)
    {
        // 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(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;
        site[x][y] = 1;
        GWputcmb(1, (float)(x), (float)(y), 1.0f, 1.0f, 0);
        dx = x - xseed;
        dy = y - yseed;
        r = (int)sqrt(dx*dx + dy*dy);       // 種からの距離
        ++(mass[r]);                        // 距離 r の位置にある点の数
        for(nn = 1; nn <= 4; nn += 1)       // 新しい周辺の点を探す
        {
            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(2, (float)(perx[nn]), (float)(pery[nn]), 1.0f, 1.0f, 0);
            }
        }
        if(x == L || x == 1)
            strcpy(boundary, "on");
        if(y == L || y == 1)
            strcpy(boundary, "on");
    } while(!(strcmp(boundary, "on") == 0));
    GWrect(0.5f, 0.5f, (float)(L+0.5), (float)(L+0.5));
                                            // 箱を再び描く
}

void mass_dist(int L, int mass[])
// r の関数として質量を出力する
{
    int masstotal = 0, r;

    masstotal = 0;
    printf("           r         M(r)\n");
    printf("\n");
    for(r = 1; r <= L; ++r)
    {
        if(mass[r] == 0)
            return;
        masstotal += mass[r];
        printf("%12d %12d\n", r, masstotal);
    }
}