// PROGRAM invasion
// 侵略型パーコレーション・クラスターの生成

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

int main();
void initial(int *Lx, int *Ly);
void assign(double site[][101], int perx[], int pery[], int *nper, int Lx, 
            int Ly);
void insert(double site[][101], int perx[], int pery[], int nper, int x, int y);
void binary_search(double site[][101], int perx[], int pery[], int nper, int x, 
                   int y, int *ninsert);
void linear_search(double site[][101], int perx[], int pery[], int nper, int x, 
                   int y, int *ninsert);
void invade(double site[][101], int perx[], int pery[], int *nper, int Lx, 
            int Ly);
void average(double site[][101], int Lx, int Ly);

int main()
{
    int perx[5001], pery[5001], Lx, Ly, nper;
    double site[201][101];

    GWopen(0);

    initial(&Lx, &Ly);
    assign(site, perx, pery, &nper, Lx, Ly);
    invade(site, perx, pery, &nper, Lx, Ly);
    average(site, Lx, Ly);

    GWquit();
    return 0;
}

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

    rnd(-1);
    printf("length in y direction = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", Ly);
    (*Lx) = 2*(*Ly);
    GWvport(0.01f*GWaspect(1), 0.2f, 0.65f*GWaspect(1), 0.99f);
    compute_aspect_ratio((float)(*Lx), &xwin, &ywin);
    GWindow(-0.5f, -0.5f, xwin-0.5f, ywin-0.5f);
    GWsetpen(BLACK, -1, -1, -1);
    GWrect(1.0f, 1.0f, (float)((*Lx)+1), (float)((*Ly)+1));
    for(y = 1; y <= (*Ly); ++y)
        for(x = 1; x <= (*Lx); ++x)
            GWsetpxl((float)(x+0.5), (float)(y+0.5), BLACK);
}

void assign(double site[][101], int perx[], int pery[], int *nper, int Lx, int Ly)
{
    int x, y;

    for(y = 1; y <= Ly; ++y)
    {
        site[1][y] = 1;                     // 最初の列を占有する
        GWsrect(1.0f, (float)(y), 2.0f, (float)(y+1), BLUE);
    }
    // 残った格子点に乱数を割り当てる
    for(y = 1; y <= Ly; ++y)
        for(x = 2; x <= Lx; ++x)
            site[x][y] = rnd(0);
    // 第2列にある格子点が最初の周辺の点となる
    // 周辺の点であるなら site(x,y) は 2 より大きい
    x = 2;
    (*nper) = 0;
    for(y = 1; y <= Ly; ++y)
    {
        site[x][y] += 2;
        ++(*nper);                          // 周辺の点の数
        // 周辺の点をソートして周辺の点の表の順序付けを行う
        insert(site, perx, pery, *nper, x, y);
    }
}

void insert(double site[][101], int perx[], int pery[], int nper, int x, int y)
// 線形探索または2分法探索のルーチンを呼ぶ
{
    int ilist, ninsert;

    binary_search(site, perx, pery, nper, x, y, &ninsert);
    // 小さい乱数の格子点のすべてをつぎの配列要素へ移す
    for(ilist = nper; ilist >= ninsert + 1; ilist += - 1)
    {
        perx[ilist] = perx[ilist - 1];
        pery[ilist] = pery[ilist - 1];
    }
    perx[ninsert] = x;                      // 新しい格子点を表に挿入する
    pery[ninsert] = y;
}

void binary_search(double site[][101], int perx[], int pery[], int nper, int x, int y, int *ninsert)
// 表を半分に分けそのどちらにその乱数が属するかを調べ，
// その乱数の位置が定まるまでこの分割を続ける
{
    int nfirst, nlast, nmid, xmid, ymid;

    nfirst = 1;                             // 表の始め
    nlast = nper - 1;                       // 表の終わり
    if(nlast < 1)
        nlast = 1;
    nmid = (nfirst + nlast)/2;              // 表の中央
    // 新しい数が表のどちらに属するかを決める
    do
    {
        if(nlast - nfirst <= 1)             // nlast に等しい正確な位置
        {
            (*ninsert) = nlast;
            return;
        }
        xmid = perx[nmid];
        ymid = pery[nmid];
        if(site[x][y] > site[xmid][ymid])
        {
            nlast = nmid;                   // 上半分を検索
        }
        else
        {
            nfirst = nmid;                  // 下半分を検索
        }
        nmid = (nfirst + nlast)/2;
    } while(1);
}

void linear_search(double site[][101], int perx[], int pery[], int nper, int x, int y, int *ninsert)
{
    int iper, xperim, yperim;

    if(nper == 1)
        (*ninsert) = 1;
    else
    {
        for(iper = 1; iper <= nper - 1; iper += 1)
        {
            xperim = perx[iper];
            yperim = pery[iper];
            if(site[x][y] > site[xperim][yperim])
            {
                (*ninsert) = iper;          // 新しい格子点を挿入
                return;
            }
        }
    }
    (*ninsert) = nper;
}

void invade(double site[][101], int perx[], int pery[], int *nper, int Lx, int Ly)
// nx と ny は最隣接点に向かうベクトル成分
{
    int i, x, xnew, y, ynew, nx[5], ny[5], DP_ = 0;
    int DATA_[] = {
        1,0,-1,0,0,1,0,-1
    };

    for(i = 1; i <= 4; ++i)
    {
        nx[i] = DATA_[DP_++];
        ny[i] = DATA_[DP_++];
    }
    do
    {
        x = perx[(*nper)];
        y = pery[(*nper)];
        --(*nper);
        // 占有された点に印をつける．その点は周辺の点ではなくなる
        --(site[x][y]);
        GWsrect((float)(x), (float)(y), (float)(x+1), (float)(y+1), BLUE);
        for(i = 1; i <= 4; ++i)          // 新しい周辺の点を探す
        {
            xnew = x + nx[i];
            ynew = y + ny[i];
            if(ynew > Ly)                   // y 方向の周期境界条件 
            {
                ynew = 1;
            }
            else if(ynew < 1)
            {
                ynew = Ly;
            }
            if(site[xnew][ynew] < 1)        // 新しい周辺の点
            {
                site[xnew][ynew] += 2;
                ++(*nper);
                insert(site, perx, pery, *nper, xnew, ynew);
            }
        }
    } while(!(x >= Lx));                // クラスターが右の境界に達するまで
}

void average(double site[][101], int Lx, int Ly)
// 確率密度 P(r) を計算する
{
    int Itmp1_, ibin, Lmax, Lmin, n, nbin, nr[21], occupied = 0, x, y;
    double dr, r, P[21];
    for(Itmp1_ = 0; Itmp1_ < 21; ++Itmp1_) {
        P[Itmp1_] = 0;
        nr[Itmp1_] = 0;
    }

    Lmin = Lx/3;
    Lmax = 2*Lx/3;
    n = (Lmax - Lmin + 1)*Ly;               // 格子の中央の半分にある格子点数
    dr = 0.05;
    nbin = (int)(1/dr);
    for(x = Lmin; x <= Lmax; ++x)
    {
        for(y = 1; y <= Ly; ++y)
        {
            ibin = (int)(nbin*fmod(site[x][y], 1.0));
            ++(nr[ibin]);
            if(site[x][y] >= 1 && site[x][y] < 2)
            {
                ++occupied;                 // 占有された格子点の総数
                ++(P[ibin]);
            }
        }
    }
    printf("# occupied sites = %d\n", occupied);
    printf("           r         P(r)\n");
    printf("\n");
    for(ibin = 0; ibin <= nbin; ++ibin)
    {
        r = dr*ibin;
        if(nr[ibin] > 0)
            printf("%12f %12f\n", r, P[ibin]/nr[ibin]);
    }
    printf("\n");
}