// 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);    // tested, unoccupied site
}

void grow(int L, int mass[])
// generate single percolation cluster          
{
    int i, iper, nn, nper, perx[5001], pery[5001], r, x, xnew, xseed, y, ynew, 
        yseed;
    int nx[5], ny[5];                   // nearest neighbor vectors
    double dx, dy, site[202][202];
    char boundary[_LBUFF_];
    int DP_ = 0;
    int DATA_[] = {
        1,0,-1,0,0,1,0,-1
    };

    strcpy(boundary, "");
    // set up 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;
    }
    // seed at center of lattice
    xseed = L/2 + 1;
    yseed = xseed;
    site[xseed][yseed] = 1;                 // seed site
    GWputcmb(1, (float)(xseed), (float)(yseed), 1.0f, 1.0f, 0);
    for(i = 1; i <= 4; i += 1)
    {
        // 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(2, (float)(perx[i]), (float)(pery[i]), 1.0f, 1.0f, 0);
    }
    nper = 4;                               // initial number of perimeter sites
    do
    {
        iper = (int)(rnd(0)*nper) + 1;      // randomly choose perimeter site
        x = perx[iper];                     // coordinate of perimeter site
        y = pery[iper];
        // last perimeter site in array replaces newly occupied site
        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);       // distance from seed
        ++(mass[r]);                        // number of sites at distance r
        for(nn = 1; nn <= 4; nn += 1)       // find new perimeter sites
        {
            xnew = x + nx[nn];
            ynew = y + ny[nn];
            if(site[xnew][ynew] == 0)
            {
                ++nper;
                perx[nper] = xnew;
                pery[nper] = ynew;
                site[xnew][ynew] = 2;       // site placed on perimeter list
                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));
                                            // redraw box
}

void mass_dist(int L, int mass[])
// print mass as a function of 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);
    }
}