#include "TrueBASIC.h"

void assign(int site[][65], int np[], int L);
void neighbor(int site[][65], int np[], int x, int y);
void label_min(int site[][65], int np[], int x, int y, int left, int down);
int proper(int np[], int label);

void assign(int site[][65], int np[], int L)
// assign cluster numbers to occupied sites
{
    int down, left, ncluster, x, y, Itmp1_;

    for(Itmp1_ = 0; Itmp1_ <= 1000; ++Itmp1_)
        np[Itmp1_] = 0;
    ncluster = 0;                           // cluster number
    for(y = 1; y <= L; ++y)
    {
        for(x = 1; x <= L; ++x)
        {
            if(site[x][y] < 0)              // site occupied
            {
                down = y - 1;               // square lattice
                left = x - 1;
                if(site[x][down] + site[left][y] == 0)
                {
                    ncluster = ncluster + 1; // new cluster
                    site[x][y] = ncluster;
                    np[ncluster] = ncluster; // proper label
                }
                else
                    neighbor(site, np, x, y);
            }
        }
    }
    // assign proper labels to cluster array
    for(y = 1; y <= L; ++y)
        for(x = 1; x <= L; ++x)
            site[x][y] = proper(np, site[x][y]);
}

void neighbor(int site[][65], int np[], int x, int y)
// determine occupancy of neighbors 
{
    int down, left;

    down = y - 1;
    left = x - 1;
    if(site[x][down]*site[left][y] > 0)     // both neighbors occupied
        label_min(site, np, x, y, left, down);
    else if(site[x][down] > 0)              // down neighbor occupied
        site[x][y] = site[x][down];
    else                                    // left neighbor occupied
        site[x][y] = site[left][y];
}

void label_min(int site[][65], int np[], int x, int y, int left, int down)
// both neighbors occupied, determine minimum cluster number
{
    int cl_down, cl_left, nmax, nmin;

    if(site[left][y] == site[x][down])
    {
        // both neighbors have same cluster label
        site[x][y] = site[left][y];
    }
    else
    {
        // determine minimum cluster label
        cl_left = proper(np, site[left][y]);
        cl_down = proper(np, site[x][down]);
        nmax = max(cl_left, cl_down);
        nmin = min(cl_left, cl_down);
        site[x][y] = nmin;
        if(nmin != nmax)
            np[nmax] = nmin;            // set improper label nmax = nmin
    }
}

int proper(int np[], int label)
{
    int proper_;

    if(np[label] == label)
        proper_ = label;
    else
        proper_ = proper(np, np[label]);
    return proper_;
}