// PROGRAM cca

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

int main();
void initial(int first_particle[], int next_particle[], int *ncl);
void neighbor(int xi, int yi, int next_particle[], int first_particle[], 
              int *ncl);
void merge(int c1, int c2, int first_particle[], int next_particle[], int *ncl);
void move(int next_particle[], int first_particle[], int *ncl);
int pbc(int s, int L);

int L;
int N;
int nx[5];
int ny[5];
int site[51][51];
int x[1001];
int y[1001];

int main()
{
    int ncl;
    int first_particle[1001], next_particle[1001];

    GWopen(0);

    initial(first_particle, next_particle, &ncl);
    do
    {
        move(next_particle, first_particle, &ncl);
    } while(ncl != 1);

    GWquit();
    return 0;
}

void initial(int first_particle[], int next_particle[], int *ncl)
{
    int i, nn, xi, yi;
    float xwin, ywin;
    char Stmp1_[_LBUFF_];
    int DP_ = 0;
    int DATA_[] = {
        1,0,-1,0,0,1,0,-1
    };

    rnd(-1);
    for(nn = 1; nn <= 4; ++nn)
    {
        nx[nn] = DATA_[DP_++];
        ny[nn] = DATA_[DP_++];
    }
    printf("L = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", &L);
    printf("N = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", &N);
    compute_aspect_ratio((float)(L+1), &xwin, &ywin);
    GWindow(0.0f, 0.0f, xwin, ywin);
    GWsetpen(BLACK, -1, -1, 4);
    GWrect(1.0f, 1.0f, (float)(L + 1), (float)(L + 1));
    sprintf(Stmp1_, "L = %d, N = %d", L, N);
    GWputtxt(1.01f, (float)(L + 1)*1.01f, Stmp1_);
    GWanchor(1);
    (*ncl) = 0;                             // クラスターの数
    for(i = 1; i <= N; ++i)
    {
        do
        {
            x[i] = (int)(L*rnd(0)) + 1;
            y[i] = (int)(L*rnd(0)) + 1;
        } while(site[x[i]][y[i]] != 0);
        ++(*ncl);
        site[x[i]][y[i]] = (*ncl);
        GWsetogn(0, (*ncl+1));
        GWsrect((float)(x[i]),   (float)(y[i]), 
                (float)(x[i]+1), (float)(y[i]+1), BLUE);
        first_particle[(*ncl)] = i;
        next_particle[i] = i;
        xi = x[i];
        yi = y[i];
        neighbor(xi, yi, next_particle, first_particle, ncl);
    }
}

void neighbor(int xi, int yi, int next_particle[], int first_particle[], 
              int *ncl)
{
    int nn, part, perim, px, py;

    for(nn = 1; nn <= 4; ++nn)
    {
        px = pbc(xi + nx[nn], L);
        py = pbc(yi + ny[nn], L);
        perim = site[px][py];
        part = site[xi][yi];
        if(perim != 0 && perim != part)
            merge(perim, part, first_particle, next_particle, ncl);
    }
}

void merge(int c1, int c2, int first_particle[], int next_particle[], int *ncl)
{
    int p, p1, p2, p2next;
    int p1next, plast;

    p1 = first_particle[c1];
    p2 = first_particle[c2];
    p1next = next_particle[p1];
    p2next = next_particle[p2];
    next_particle[p1] = p2next;
    next_particle[p2] = p1next;
    do
    {
        site[x[p2next]][y[p2next]] = c1;
        p2next = next_particle[p2next];
    } while(p2next != p1next);
    GWsetogn(c2+1, c1+1);
    plast = first_particle[(*ncl)];
    if(c2 != (*ncl))
    {
        p = plast;
        do
        {
            site[x[p]][y[p]] = c2;
            p = next_particle[p];
        } while(p != plast);
        first_particle[c2] = plast;
        GWsetogn((*ncl+1), c2+1);
    }
    --(*ncl);
}

void move(int next_particle[], int first_particle[], int *ncl)
{
    int c, direction, i, xi, yi;
    int dx, dy, p1;

    c = (int)((*ncl)*rnd(0)) + 1;
    direction = (int)(4*rnd(0)) + 1;
    p1 = first_particle[c];
    i = p1;
    dx = nx[direction];
    dy = ny[direction];
    do
    {
        site[x[i]][y[i]] = 0;
        x[i] = pbc(x[i] + dx, L);
        y[i] = pbc(y[i] + dy, L);
        i = next_particle[i];
    } while(i != p1);
    GWsetogn(0, -c-1);
    do
    {
        GWsrect((float)(x[i]),   (float)(y[i]), 
                (float)(x[i]+1), (float)(y[i]+1), BLUE);
        site[x[i]][y[i]] = c;
        i = next_particle[i];
    } while(i != p1);
    GWflush(-c-1);
    GWsetogn(0, 1);
    do
    {
        xi = x[i];
        yi = y[i];
        neighbor(xi, yi, next_particle, first_particle, ncl);
        i = next_particle[i];
    } while(i != p1);
}

int pbc(int s, int L)
{
    int pbc_;

    if(s > L)
        pbc_ = 1;
    else if(s < 1)
        pbc_ = L;
    else
        pbc_ = s;
    return pbc_;
}