// PROGRAM continuum
// generate continuum percolation configurations

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

int main();
void initial(int *L);
void disks(int L);

int main()
{
    int L;

    GWopen(0);

    initial(&L);
    disks(L);                        // place disks at random inside box

    GWquit();
    return 0;
}

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

    rnd(-1);
    GWinput("linear dimension of lattice = ", Stmp1_, _LBUFF_);
    sscanf(Stmp1_, "%d", L);
    compute_aspect_ratio((float)(*L), &xwin, &ywin);
    GWindow(-0.1f*xwin, -0.05f*ywin, xwin, 1.05f*ywin);
    GWcolor(BLACK, CLR_BACKGRND);
    GWcolor(WHITE, CLR_TEXT);
    GWsetpen(RED, -1, -1, 4);               // erase disk
    GWsetbrs(RED, 0, -1);
    GWrect(0.0f, 0.0f, (float)((*L)), (float)((*L)));
    GWanchor(1);
}

void disks(int L)
{
    int i, n = 0, number = 0;
    float xs, ys;
    double area, density, phi, r, x[10001], y[10001];
    char Stmp1_[_LBUFF_];

    r = 1;                                  // radius of disks
    area = L*L;
    GWsetpen(YELLOW, -1, -1, 4);
    GWsetmrk(6, 1.0f, YELLOW, -1, 4);
    while(n >= 0)
    {
        GWsetogn(0, -10);
        if(!GWinput("number of additional disks = ", Stmp1_, _LBUFF_)) break;
        sscanf(Stmp1_, "%d", &n);
        density = (n+number)/area;
        phi = 1 - exp(-density*pi*r*r);
        sprintf(Stmp1_, "number of disks = %d", n + number);
        GWputtxt(0.0f, (float)(L*25.0/24), Stmp1_);
        sprintf(Stmp1_, "density of disks = %f", density);
        GWputtxt(0.0f, (float)(L*27.0/24), Stmp1_);
        sprintf(Stmp1_, "phi = %f", phi);
        GWputtxt(0.0f, (float)(L*26.0/24), Stmp1_);
        for(i = 1; i <= number; ++i)
            GWputmrk((float)(x[i]), (float)(y[i]));
        for(i = number+1; i <= n+number; ++i)
        {
            x[i] = L*rnd(0);                // center of disk inside box
            y[i] = L*rnd(0);
            GWputmrk((float)(x[i]), (float)(y[i]));
        }
        number += n;
        GWflush(-10);
        GWsetogn(0, 2);
        GWsetbrs(BLUE, 1, -1);
        do
        {
            GWcappnt(&xs, &ys, NULL);
            // click on occupied site to see cluster
            // click mouse to right of lattice to increase p by delta
            if(xs < L && ys <= L)
            {
                GWflood(xs, ys);
            }
        } while(!(xs > L));
        GWerase(2, -1);
    }
}