//  SINGLE.C
//  Equilibration using single spin flip
//  © 2021 Peter J. Meyer

#include "iss.h"

static int process_spin(int i, int j, ...);

/*--------------------------------*/
void do_checkerboard_ssf_sweep(void)
{
int num_sublattices, sublattice_size;
int i, j, k, l, i1, i2, j1, j2, k1, k2, l1,l2;

//  Use "checkerboard algorithm" for order in which to attempt spin flips.

//  Randomly select a divisor from up to 7 divisors selected in SETPARAM.C.
num_sublattices = divisors[(int)(PRNG1()*num_divisors)];
sublattice_size = size/num_sublattices;

switch ( dimensionality )
    {
    case 2:
    for ( i1=0; i1<sublattice_size; i1++ )
        {
    for ( j1=0; j1<sublattice_size; j1++ )
        {
        for ( i2=0; i2<num_sublattices; i2++ )
            {
            i = i2*sublattice_size + i1;
            for ( j2=0; j2<num_sublattices; j2++ )
                {
                j = j2*sublattice_size + j1;
                process_spin(i,j);
                } 
            }
        }
        }
    break;

    case 3:
    for ( i1=0; i1<sublattice_size; i1++ )
        {
    for ( j1=0; j1<sublattice_size; j1++ )
        {
    for ( k1=0; k1<sublattice_size; k1++ )
        {
        for ( i2=0; i2<num_sublattices; i2++ )
            {
            i = i2*sublattice_size + i1;
            for ( j2=0; j2<num_sublattices; j2++ )
                {
                j = j2*sublattice_size + j1;
                for ( k2=0; k2<num_sublattices; k2 ++ )
                    {
                    k = k2*sublattice_size + k1;
                    process_spin(i,j,k);
                    } 
                }
            }
        }
        }
        }
    break;

    case 4:
    for ( i1=0; i1<sublattice_size; i1++ )
        {
    for ( j1=0; j1<sublattice_size; j1++ )
        {
    for ( k1=0; k1<sublattice_size; k1++ )
        {
    for ( l1=0; l1<sublattice_size; l1++ )
        {
        for ( i2=0; i2<num_sublattices; i2++ )
            {
            i = i2*sublattice_size + i1;
            for ( j2=0; j2<num_sublattices; j2++ )
                {
                j = j2*sublattice_size + j1;
                for ( k2=0; k2<num_sublattices; k2++ )
                    {
                    k = k2*sublattice_size + k1;
                    for ( l2=0; l2<num_sublattices; l2++ )
                        {
                        l = l2*sublattice_size + l1;
                        process_spin(i,j,k,l);
                        }
                    }
                }
            }
        }
        }
        }
        }
    break;
    }
}

// This increases time required (compared to Glauber) by about 75%.
/*--------------------------*/
void do_random_ssf_sweep(void)
{
int n, i, j, k, l;

//  Pick num_spins sites at random.
n = num_spins;

do  {
    switch ( dimensionality )
        {
        case 2:
        i = (int)(PRNG1()*size);
        j = (int)(PRNG1()*size);
        n -= process_spin(i,j);     //  n decremented iff site is occupied
        break;

        case 3:
        i = (int)(PRNG1()*size);
        j = (int)(PRNG1()*size);
        k = (int)(PRNG1()*size);
        n -= process_spin(i,j,k);
        break;

        case 4:
        i = (int)(PRNG1()*size);
        j = (int)(PRNG1()*size);
        k = (int)(PRNG1()*size);
        l = (int)(PRNG1()*size);
        n -= process_spin(i,j,k,l);
        break;
        }
    } while ( n );
}

/*--------------------------------------*/
static int process_spin(int i, int j, ...)
{
int spin_value, new_spin_value, k, l;

num_spins_visited_this_sweep++;

switch ( dimensionality )
    {
    case 2:
    //  spin_at_site(i,j) returns 0 if site is unoccupied,
    //  otherwise returns the spin value.
    if ( spin_value = spin_at_site(i,j) )       
        {
        if ( model_is_ising )
            {
            if ( flip_spin(2*spin_value*nn_spin_sum(i,j)) )
                do_ising_spin_flip(spin_value,i,j);         
            }
        else    //  model is q-state Potts
            {
            //  Pick another spin value.
            do  {
                new_spin_value = (int)(PRNG1()*q) + 1;
                } while ( new_spin_value == spin_value );
            if ( flip_spin(sum_kronecker_deltas(spin_value,i,j)
                     - sum_kronecker_deltas(new_spin_value,i,j)) )
                do_q_state_potts_spin_flip(spin_value, 
                    new_spin_value,i,j);
            }
        }
    break;

    case 3:
    k = *(((int *)(&j))+1);
    if ( spin_value = spin_at_site(i,j,k) )
        {
        if ( model_is_ising )
            {
            if ( flip_spin(2*spin_value*nn_spin_sum(i,j,k)) )
                do_ising_spin_flip(spin_value,i,j,k);
            }
        else    //  model is q-state Potts
            {
            //  Pick another spin value.
            do  {
                new_spin_value = (int)(PRNG1()*q) + 1;
                } while ( new_spin_value == spin_value );
            if ( flip_spin(sum_kronecker_deltas(spin_value,i,j,k)
                 - sum_kronecker_deltas(new_spin_value,i,j,k)) )
                do_q_state_potts_spin_flip(spin_value, 
                    new_spin_value,i,j,k);
            }
        }
    break;

    case 4:
    k = *(((int *)(&j))+1);
    l = *(((int *)(&j))+2);
    if ( spin_value = spin_at_site(i,j,k,l) )
        {
        if ( model_is_ising )
            {
            if ( flip_spin(2*spin_value*nn_spin_sum(i,j,k,l)) )
                do_ising_spin_flip(spin_value,i,j,k,l);
            }
        else    //  model is q-state Potts
            {
            //  Pick another spin value.
            do  {
                new_spin_value = (int)(PRNG1()*q) + 1;
                } while ( new_spin_value == spin_value );
            if ( flip_spin(sum_kronecker_deltas(spin_value,i,j,k,l)
                 - sum_kronecker_deltas(new_spin_value,i,j,k,l)) )
                do_q_state_potts_spin_flip(spin_value, 
                    new_spin_value,i,j,k,l);
            }
        }
    }

return ( spin_value );  //  returns false iff site is not occupied.
}