//  S&B.C
//  Sites and bonds.
//  Functions for initialization of the spin model
//  © 2021 Peter J. Meyer

#include "iss.h"

unsigned int total_num_opened_bonds = 0;

void initialize_lattice_and_mark_boundary_sites(void);

static int de_occupy_sites(void);

//  On entry to this function the sites array has been initialized
//  showing all sites occupied by up-spins.
//  This function (a) marks the boundary sites of the lattice,
//  (b) optionally precomputes nn sites,
//  (b) in a site-diluted lattice it marks some sites as unoccupied,
//  (c) it opens bonds between sites (randomly if the lattice is bond-diluted).
//  Upon return the sites array shows the occupied sites 
//  and the bonds array shows the open bonds.
//  This assumes that the PRNG has been initialized.
/*----------------------------*/
void set_up_configuration(void)
{
int i;
double site_pc, bond_pc;

initialize_lattice_and_mark_boundary_sites();

num_spins = ( site_diluted ? de_occupy_sites() : num_sites );

num_opened_bonds = open_bonds();        //  BONDS.C
total_num_opened_bonds += num_opened_bonds;

if ( goal_is_equilibration )
    {
    site_pc = ((double)num_spins*100)/num_sites;

    if ( site_diluted )
        {
        printf("\n%s of %s (%d",ultoa_commas(num_spins,temp),ultoa_commas(num_sites,temp1),size);

        for ( i=2; i<=dimensionality; i++ )
            printf("x%d",size);

        printf(") lattice sites occupied (%.2f%%).",site_pc);
        }

    bond_pc = ((double)num_opened_bonds*100)/num_pairs_nn;

    if ( bond_diluted )
        {
        printf("\n%s bonds opened. This is %.2f%% of the %s pairs of nn spins.",
            ultoa_commas(num_opened_bonds,temp),bond_pc,ultoa_commas(num_pairs_nn,temp2));

        #if false
        sprintf(temp,"%.3f",(double)num_opened_bonds/num_spins);
        remove_trailing_zeros(temp);
        printf("\nBonds per spin = %s",temp);
        #endif
        }

    printf("\n");

    if ( multiple_temperatures_specified )
        {
        sprintf(temp,"%.6f",temperature);
        remove_trailing_zeros(temp);
        printf("Temperature = %s  ",temp);
        //  printf(" (%d/%d)",temperature_num+1,num_temperatures);
        }
    }

//  At this point all spins are up-spins.
//  For equilibration the assignment of initial spins
//  is done by assign_initial_spins() in INITSPIN.C.
}

//  On entry all sites are marked as occupied by up-spins.
//  This function marks some sites as unoccupied.
//  It is independent of the lattice type.
/*----------------------------*/
static int de_occupy_sites(void)
{
int i, j, k, l;                     //  array indices 
int num_occupied = num_sites;
int target_num_occupied = (int)(num_sites*site_concentration);

switch ( dimensionality )
    {
    case 2: 
    for ( i=0; i<size; i++ )
        {
        for ( j=0; j<size; j++ )
            {
            if ( PRNG1() >= site_concentration )                     
                {                                      
                sites2[i][j] &= UNOCCUPIED;
                num_occupied--;
                }
            }
        }

    //  If necessary, adjust site occupancy to get exact concentration specified.

    while ( num_occupied > target_num_occupied )
        {
        i = (int)(size*PRNG1());
        j = (int)(size*PRNG1());
        if ( SITE_IS_OCCUPIED_2(i,j) )
            {
            //  Site is occupied, so de-occupy it.
            sites2[i][j] &= UNOCCUPIED;
            num_occupied--;
            }
        }

    while ( num_occupied < target_num_occupied )
        {
        i = (int)(size*PRNG1());
        j = (int)(size*PRNG1());
        if ( !SITE_IS_OCCUPIED_2(i,j) )
            {
            //  Site is not occupied, so occupy it.
            sites2[i][j] |= OCCUPIED;
            num_occupied++;
            }
        }

    break;

    case 3:
    for ( i=0; i<size; i++ )
        {
        for ( j=0; j<size; j++ )
            {
            for ( k=0; k<size; k++ )
                {
                if ( PRNG1() >= site_concentration )                     
                    {                                      
                    sites3[i][j][k] &= UNOCCUPIED;
                    num_occupied--;
                    }
                }
            }
        }

    while ( num_occupied > target_num_occupied )
        {
        i = (int)(size*PRNG1());
        j = (int)(size*PRNG1());
        k = (int)(size*PRNG1());
        if ( SITE_IS_OCCUPIED_3(i,j,k) )
            {
            sites3[i][j][k] &= UNOCCUPIED;
            num_occupied--;
            }
        }

    while ( num_occupied < target_num_occupied )
        {
        i = (int)(size*PRNG1());
        j = (int)(size*PRNG1());
        k = (int)(size*PRNG1());
        if ( !SITE_IS_OCCUPIED_3(i,j,k) )
            {
            sites3[i][j][k] |= OCCUPIED;
            num_occupied++;
            }
        }

    break;
                        
    case 4:
    for ( i=0; i<size; i++ )
        {
        for ( j=0; j<size; j++ )
            {
            for ( k=0; k<size; k++ )
                {
                for ( l=0; l<size; l++ )
                    {
                    if ( PRNG1() >= site_concentration )                     
                        {                                      
                        sites4[i][j][k][l] &= UNOCCUPIED;
                        num_occupied--;
                        }
                    }
                }
            }
        } 
        
    while ( num_occupied > target_num_occupied )
        {
        i = (int)(size*PRNG1());
        j = (int)(size*PRNG1());
        k = (int)(size*PRNG1());
        l = (int)(size*PRNG1());
        if ( SITE_IS_OCCUPIED_4(i,j,k,l) )
            {
            sites4[i][j][k][l] &= UNOCCUPIED;
            num_occupied--;
            }
        }

    while ( num_occupied < target_num_occupied )
        {
        i = (int)(size*PRNG1());
        j = (int)(size*PRNG1());
        k = (int)(size*PRNG1());
        l = (int)(size*PRNG1());
        if ( !SITE_IS_OCCUPIED_4(i,j,k,l) )
            {
            sites4[i][j][k][l] |= OCCUPIED;
            num_occupied++;
            }
        }

    break;                           
    }

return ( num_occupied );
}

//  open_bonds() is in BONDS.C

/*--------------------------------------------------*/
void initialize_lattice_and_mark_boundary_sites(void)
{
int i, j, k;
unsigned char site_val = ( OCCUPIED_BY_UP_SPIN | BOUNDARY );    //  Bit pattern 11000000

switch ( dimensionality )
    {
    case 2:  
    //  First show all sites as occupied by up spins
    //  with all lattice bonds closed.
    memset(&sites2[0][0],OCCUPIED_BY_UP_SPIN,num_sites);
    memset(&bonds2[0][0],0,num_sites);

    //  Now mark the boundary sites.                             
    for ( i=0; i<size; i++ )
        {
        sites2[0][i] = sites2[size-1][i] = 
        sites2[i][0] = sites2[i][size-1] = site_val;
        }
    break;

    case 3:
    memset(&sites3[0][0][0],OCCUPIED_BY_UP_SPIN,num_sites);
    memset(&bonds3[0][0][0],0,num_sites);

    for ( i=0; i<size; i++ )
        {
        for ( j=0; j<size; j++ )
            {
            sites3[0][i][j] = sites3[size-1][i][j] = 
            sites3[i][0][j] = sites3[i][size-1][j] = 
            sites3[i][j][0] = sites3[i][j][size-1] = site_val;
            }
        }
    break;

    case 4:
    memset(&sites4[0][0][0][0],OCCUPIED_BY_UP_SPIN,num_sites);
    memset(&bonds4[0][0][0][0],0,num_sites);

    for ( i=0; i<size; i++ )
        {
        for ( j=0; j<size; j++ )
            {
            for ( k=0; k<size; k++ )
                {
                sites4[0][i][j][k] = sites4[size-1][i][j][k] = 
                sites4[i][0][j][k] = sites4[i][size-1][j][k] = 
                sites4[i][j][0][k] = sites4[i][j][size-1][k] = 
                sites4[i][j][k][0] = sites4[i][j][k][size-1] = site_val;
                }
            }
        }
    }
}
 
/*-------------------------------*/
void mark_all_sites_unvisited(void)
{
int i, j, k, l;  

switch ( dimensionality )
    {
    case 2:                   
    for ( i=0; i<size; i++ )
        for ( j=0; j<size; j++ )
            sites2[i][j] &= UNVISITED;
    break;

    case 3:                   
    for ( i=0; i<size; i++ )
        for ( j=0; j<size; j++ )
            for ( k=0; k<size; k++ )
                sites3[i][j][k] &= UNVISITED;
    break;

    case 4:                   
    for ( i=0; i<size; i++ )
        for ( j=0; j<size; j++ )
            for ( k=0; k<size; k++ )
                for (l=0; l<size; l++ )
                    sites4[i][j][k][l] &= UNVISITED;
    break;
    }
}

/*----------------------------------------------*/
void mark_all_sites_unvisited_and_clear_path(void)
{
int i, j, k, l;  
unsigned char site_val = ( UNVISITED & NOT_IN_SPANNING_CLUSTER );  //  Bit pattern 11001111

switch ( dimensionality )
    {
    case 2:                   
    for ( i=0; i<size; i++ )
        for ( j=0; j<size; j++ )
            sites2[i][j] &= site_val;
    break;

    case 3:                   
    for ( i=0; i<size; i++ )
        for ( j=0; j<size; j++ )
            for ( k=0; k<size; k++ )
                sites3[i][j][k] &= site_val;
    break;

    case 4:                   
    for ( i=0; i<size; i++ )
        for ( j=0; j<size; j++ )
            for ( k=0; k<size; k++ )
                for (l=0; l<size; l++ )
                    sites4[i][j][k][l] &= site_val;
    break;
    }
}

#if false
//  No longer needed.

//  This is used in the Wolff algorithm.
/*-----------------*/
void copy_sites(void)
{
switch ( dimensionality )
    {
    case 2:
    memcpy(&v_sites2[0][0], &sites2[0][0], num_sites);
    break;

    case 3:
    memcpy(&v_sites3[0][0][0], &sites3[0][0][0], num_sites);
    break;

    case 4:
    memcpy(&v_sites4[0][0][0][0], &sites4[0][0][0][0], num_sites);
    break;
    }
}
#endif