// PROGRAM perc_cluster
// cluster generated by Hammersley, Leath, and Alexandrowicz algorithm

#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);     // occupied site
    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);    // perimeter site
    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);   // tested, unoccupied site
    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[])
// generate single percolation cluster          
{
    int i, iper, nn, nper, perx[25001], pery[25001], x, xnew, xseed, y, ynew, 
        yseed;
    int nx[5], ny[5];                       // set up direction vectors for lattice
    int site[132][132];
    int DP_ = 0;
    int DATA_[] = {
        1,0,-1,0,0,1,0,-1
    };

    // set up boundary sites
    for(i = 1; i <= L; ++i)
    {
        site[0][i] = -1;
        site[L + 1][i] = -1;
        site[i][L + 1] = -1;
        site[i][0] = -1;
    }
    // seed at center of lattice
    xseed = L/2 + 1;
    yseed = xseed;
    site[xseed][yseed] = 1;                 // seed site
    xs[1] = xseed;
    ys[1] = yseed;
    (*N) = 1;                               // number of sites in the cluster
    GWputcmb(status[1], (float)(xseed), (float)(yseed), 1.0f, 1.0f, 0);
    for(i = 1; i <= 4; ++i)
    {
        // nx,ny direction vectors for new perimeter sites
        nx[i] = DATA_[DP_++];
        ny[i] = DATA_[DP_++];
        // perx,pery, positions of perimeter sites of seed
        perx[i] = xseed + nx[i];
        pery[i] = yseed + ny[i];
        // perimeter sites labeled by 2
        site[perx[i]][pery[i]] = 2;         // site placed on perimeter list
        GWputcmb(status[2], (float)(perx[i]), (float)(pery[i]), 1.0f, 1.0f, 0);
    }
    nper = 4;                               // initial number of perimeter sites
    do
    {
        // randomly choose perimeter site
        iper = (int)(rnd(0)*nper) + 1;
        x = perx[iper];                     // coordinate of a perimeter site
        y = pery[iper];
        // relabel remaining perimeter sites so that
        // last perimeter site in array replaces newly chosen site
        perx[iper] = perx[nper];
        pery[iper] = pery[nper];
        --nper;
        if(rnd(0) < p)                      // site occupied
        {
            site[x][y] = 1;
            ++(*N);
            xs[(*N)] = x;                   // save position of occupied site
            ys[(*N)] = y;
            GWputcmb(status[1], (float)(x), (float)(y), 1.0f, 1.0f, 0);
            for(nn = 1; nn <= 4; ++nn)   // find new perimeter sites
            {
                xnew = x + nx[nn];
                ynew = y + ny[nn];
                if(site[xnew][ynew] == 0)
                {
                    ++nper;
                    perx[nper] = xnew;
                    pery[nper] = ynew;
                    // place site on perimeter list
                    site[xnew][ynew] = 2;
                    GWputcmb(status[2], (float)(xnew), (float)(ynew), 1.0f, 1.0f, 0);
                }                           // rnd >= p so site is not occupied        
            }
        }
        else
        {
            site[x][y] = -1;
            GWputcmb(status[0], (float)(x), (float)(y), 1.0f, 1.0f, 0);
        }
    } while(!(nper < 1));                   // all perimeter sites tested
    GWsetbrs(BLACK, 0, -1);
    GWrect(0.5f, 0.5f, (float)(L+0.5), (float)(L+0.5));
                                            // redraw box
}

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];                       // compute center of mass
        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));
        // distance from center of mass
        // mass(r) = number of sites at distance r from center of mass
        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;                    // use logarithmic scale for r
        }
    }
}