// PROGRAM invasion
// generate invasion percolation cluster

#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;                     // occupy first column
        GWsrect(1.0f, (float)(y), 2.0f, (float)(y+1), BLUE);
    }
    // assign random numbers to remaining sites
    for(y = 1; y <= Ly; ++y)
        for(x = 2; x <= Lx; ++x)
            site[x][y] = rnd(0);
    // sites in second column are initial perimeter sites
    // site(x,y) greater than 2 if perimeter site
    x = 2;
    (*nper) = 0;
    for(y = 1; y <= Ly; ++y)
    {
        site[x][y] += 2;
        ++(*nper);                          // number of perimeter sites
        // sort perimeter sites
        // order perimeter list
        insert(site, perx, pery, *nper, x, y);
    }
}

void insert(double site[][101], int perx[], int pery[], int nper, int x, int y)
// call linear or binary sort
{
    int ilist, ninsert;

    binary_search(site, perx, pery, nper, x, y, &ninsert);
    // move sites with smaller random numbers to next higher array index
    for(ilist = nper; ilist >= ninsert + 1; ilist += - 1)
    {
        perx[ilist] = perx[ilist - 1];
        pery[ilist] = pery[ilist - 1];
    }
    perx[ninsert] = x;                      // new site inserted in list
    pery[ninsert] = y;
}

void binary_search(double site[][101], int perx[], int pery[], int nper, int x, int y, int *ninsert)
// divide list in half and determine half in which random # belongs
// continue this division until position of number is determined
{
    int nfirst, nlast, nmid, xmid, ymid;

    nfirst = 1;                             // beginning of list
    nlast = nper - 1;                       // end of list
    if(nlast < 1)
        nlast = 1;
    nmid = (nfirst + nlast)/2;              // middle of list
    // determine which half of list new number is located
    do
    {
        if(nlast - nfirst <= 1)             // exact position equal to nlast
        {
            (*ninsert) = nlast;
            return;
        }
        xmid = perx[nmid];
        ymid = pery[nmid];
        if(site[x][y] > site[xmid][ymid])
        {
            nlast = nmid;                   // search upper half
        }
        else
        {
            nfirst = nmid;                  // search lower half
        }
        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;          // insert new site
                return;
            }
        }
    }
    (*ninsert) = nper;
}

void invade(double site[][101], int perx[], int pery[], int *nper, int Lx, int Ly)
// nx and ny are components of vectors pointing to nearest neighbors
{
    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);
        // mark site occupied and no longer perimeter site
        --(site[x][y]);
        GWsrect((float)(x), (float)(y), (float)(x+1), (float)(y+1), BLUE);
        for(i = 1; i <= 4; ++i)          // find new perimeter sites          
        {
            xnew = x + nx[i];
            ynew = y + ny[i];
            if(ynew > Ly)                   // periodic boundary conditions in y 
            {
                ynew = 1;
            }
            else if(ynew < 1)
            {
                ynew = Ly;
            }
            if(site[xnew][ynew] < 1)        // new perimeter site
            {
                site[xnew][ynew] += 2;
                ++(*nper);
                insert(site, perx, pery, *nper, xnew, ynew);
            }
        }
    } while(!(x >= Lx));                // until cluster reaches right boundary
}

void average(double site[][101], int Lx, int Ly)
// compute probability density 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;               // # sites in middle half of lattice
    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;                 // total # of occupied sites
                ++(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");
}