//  Q&A.C
//  Questions and answers
//  © 2021 Peter J. Meyer

#include "iss.h"

unsigned int num_times_get_nn_sites_called = 0;
unsigned int num_times_get_nn_site_in_dir_called  = 0;
unsigned int num_times_adjust_index_called = 0;

static int adjust_index(short int *ip);

//  What are the nearest neighbour sites of a given site 
//  (whether or not occupied and whether or not an open bond exists)?
//
//  For a given site (i,j,..), nn_sites[r][n] is index n of the nearest neighbour
//  site in direction r (where 0 <= n < dimensionality).
//  If periodic boundary conditions are used, or the site is not a boundary site,
//  then the number of nn sites equals the coordination number.
//  If free boundary conditions are used and the site is a boundary site
//  then nn_sites[r][0] will be NON_EXISTENT for some r.
//  This function does not check which nearest neighbour sites 
//  are occupied or whether open bonds to them exist.
/*--------------------------------*/
void get_nn_sites(int i, int j, ...)
{
int k, l, r, index, boundary_site;

num_times_get_nn_sites_called++;

if ( use_precomputed_nn_sites )
    {
    switch ( dimensionality )
        {
        case 2:
        for ( r=0; r<coord_num; r++ )           
            for ( index=0; index<dimensionality; index++ )
                nn_sites[r][index] = precomp_nn_sites2[i][j][r][index];     
        break;

        case 3:
        k = *(((int *)(&j))+1);
        for ( r=0; r<coord_num; r++ )           
            for ( index=0; index<dimensionality; index++ )
                nn_sites[r][index] = precomp_nn_sites3[i][j][k][r][index];     
        break;

        case 4:
        k = *(((int *)(&j))+1);
        l = *(((int *)(&j))+2);
        for ( r=0; r<coord_num; r++ )           
            for ( index=0; index<dimensionality; index++ )
                nn_sites[r][index] = precomp_nn_sites4[i][j][k][l][r][index];     
        break;
        }
    }
else
    {
    switch ( dimensionality )
        {
        case 2:
        if ( direction_table_is_parity_dependent )
            adjust_direction_table((i+j)%2);
        for ( r=0; r<coord_num; r++ )           //  For each direction
            {
            nn_sites[r][0] = i + dir[r][0];     //  See DIRTABLE.C
            nn_sites[r][1] = j + dir[r][1];
            }
        boundary_site = IS_BOUNDARY_SITE_2(i,j);
        break;

        case 3:
        k = *(((int *)(&j))+1);
        if ( direction_table_is_parity_dependent )
            adjust_direction_table((i+j+k)%2); 
        for ( r=0; r<coord_num; r++ )
            {
            nn_sites[r][0] = i + dir[r][0];
            nn_sites[r][1] = j + dir[r][1];
            nn_sites[r][2] = k + dir[r][2];
            }
        boundary_site = IS_BOUNDARY_SITE_3(i,j,k);
        break;

        case 4:
        k = *(((int *)(&j))+1);
        l = *(((int *)(&j))+2);
        if ( direction_table_is_parity_dependent )
            adjust_direction_table((i+j+k+l)%2); 
        for ( r=0; r<coord_num; r++ )
            {
            nn_sites[r][0] = i + dir[r][0];
            nn_sites[r][1] = j + dir[r][1];
            nn_sites[r][2] = k + dir[r][2];
            nn_sites[r][3] = l + dir[r][3];
            }
        boundary_site = IS_BOUNDARY_SITE_4(i,j,k,l);
        break;
        }

    if ( boundary_site )
        {
        //  size -> 0 and -1 -> size-1
        for ( r=0; r<coord_num; r++ )
            {
            for ( index=0; index<dimensionality; index++ )
                {
                #if FREE_BOUNDARY_CONDITIONS_POSSIBLE            
                if ( adjust_index(&nn_sites[r][index]) )
                //  Returns true if index was adjusted, i.e. if index was out of bounds.
                    {
                    if ( free_boundary_conditions )     //  No longer occurs
                        {
                        nn_sites[r][0] = NON_EXISTENT;
                        //  Meaning that there is no nn in direction r.
                        break;
                        }
                    } 
                #else
                adjust_index(&nn_sites[r][index]);
                //  With periodic boundary conditions there are no non-existent nn's
                //  so no need to check return value.
                #endif                               
                }
            }
        }
    }
}

//  This assumes that the maximum distance in any direction 
//  of a nearest neighbour is not greater than the lattice size.
/*----------------------------------*/
static int adjust_index(short int *ip)
{
num_times_adjust_index_called++;

if ( *ip < 0 )
    {
    *ip += size;
    return ( true );
    }
else if ( *ip >= size )
    {
    *ip -= size;
    return ( true );
    }
else
    return ( false );
}

//  What is the nearest neighbour site of a given site in a given direction?  
//  Returns true if site exists, false otherwise.
/*--------------------------------------------*/
int get_nn_site_in_dir(int r, int i, int j, ...)
{
int k, l, index, boundary_site;
int exists = true;


num_times_get_nn_site_in_dir_called++;

if ( use_precomputed_nn_sites )
    {
    switch ( dimensionality )
        {
        case 2:          
        for ( index=0; index<dimensionality; index++ )
            nn_site[index] = precomp_nn_sites2[i][j][r][index];   
        exists = ( nn_site[0] != NON_EXISTENT );              
        break;

        case 3:          
        k = *(((int *)(&j))+1);
        for ( index=0; index<dimensionality; index++ )
            nn_site[index] = precomp_nn_sites3[i][j][k][r][index];   
        exists = ( nn_site[0] != NON_EXISTENT );              
        break;

        case 4:          
        k = *(((int *)(&j))+1);
        l = *(((int *)(&j))+2);
        for ( index=0; index<dimensionality; index++ )
            nn_site[index] = precomp_nn_sites4[i][j][k][l][r][index];   
        exists = ( nn_site[0] != NON_EXISTENT );              
        break;
        }
    }
else
    {
    switch ( dimensionality )
        {
        case 2:
        if ( direction_table_is_parity_dependent )
            adjust_direction_table((i+j)%2);

        nn_site[0] = i + dir[r][0];     //  See DIRTABLE.C
        nn_site[1] = j + dir[r][1];

        boundary_site = IS_BOUNDARY_SITE_2(i,j);
        break;

        case 3:
        k = *(((int *)(&j))+1);

        if ( direction_table_is_parity_dependent )
            adjust_direction_table((i+j+k)%2); 

        nn_site[0] = i + dir[r][0];
        nn_site[1] = j + dir[r][1];
        nn_site[2] = k + dir[r][2];

        boundary_site = IS_BOUNDARY_SITE_3(i,j,k);
        break;

        case 4:
        k = *(((int *)(&j))+1);
        l = *(((int *)(&j))+2);

        if ( direction_table_is_parity_dependent )
            adjust_direction_table((i+j+k+l)%2); 
        nn_site[0] = i + dir[r][0];
        nn_site[1] = j + dir[r][1];
        nn_site[2] = k + dir[r][2];
        nn_site[3] = l + dir[r][3];

        boundary_site = IS_BOUNDARY_SITE_4(i,j,k,l);
        break;
        }

    if ( boundary_site )
        {
        //  size -> 0 and -1 -> size-1
        for ( index=0; index<dimensionality; index++ )
            {
            #if FREE_BOUNDARY_CONDITIONS_POSSIBLE
            if ( adjust_index(&nn_site[index]) )
            //  Returns true if index was adjusted, i.e. if index was out of bounds.
                {
                if ( free_boundary_conditions )     
                    {
                    nn_site[0] = NON_EXISTENT;
                    //  Meaning that there is no nn in direction r.
                    exists = false;
                    break;
                    }
                } 
            #else
            adjust_index(&nn_site[index]);
            //  With periodic boundary conditions there are no non-existent nn's
            //  so no need to check return value.
            #endif               
            }
        }
    }

return ( exists );  
}

#if false
//  What is the site at a distance n from a given site
//  in a given direction?  
//  Returns true if site exists, false otherwise.
//  THIS FUNCTION NEEDS RE-THINKING.
/*---------------------------------------------------------------*/
int get_nn_site_at_distance_in_dir(int n, int r, int i, int j, ...)
{
int k, l, index, boundary_site;
int exists = true;

switch ( dimensionality )
    {
    case 2:
    nn_site[0] = i;
    nn_site[1] = j;

    while ( n-- )
        {
        if ( direction_table_is_parity_dependent )
            adjust_direction_table((nn_site[0]+nn_site[1])%2);  //  ???
        nn_site[0] += dir[r][0];     //  See DIRTABLE.C
        nn_site[1] += dir[r][1];
        }

    boundary_site = IS_BOUNDARY_SITE_2(i,j);
    break;

    case 3:
    k = *(((int *)(&j))+1);

    nn_site[0] = i;
    nn_site[1] = j;
    nn_site[2] = k;

    while ( n-- )
        {
        if ( direction_table_is_parity_dependent )
            adjust_direction_table((nn_site[0]+nn_site[1]+nn_site[2])%2); 
        nn_site[0] += dir[r][0];
        nn_site[1] += dir[r][1];
        nn_site[2] += dir[r][2];
        }

    boundary_site = IS_BOUNDARY_SITE_3(i,j,k);
    break;

    case 4:
    k = *(((int *)(&j))+1);
    l = *(((int *)(&j))+2);

    nn_site[0] = i;
    nn_site[1] = j;
    nn_site[2] = k;
    nn_site[3] = l;

    while ( n-- )
        {
        if ( direction_table_is_parity_dependent )
            adjust_direction_table((nn_site[0]+nn_site[1]+nn_site[2]+nn_site[3])%2); 
        nn_site[0] += dir[r][0];
        nn_site[1] += dir[r][1];
        nn_site[2] += dir[r][2];
        nn_site[3] += dir[r][3];
        }

    boundary_site = IS_BOUNDARY_SITE_4(i,j,k,l);
    break;
    }

if ( boundary_site )
    {
    //  size -> 0 and -1 -> size-1
    for ( index=0; index<dimensionality; index++ )
        {
        #if FREE_BOUNDARY_CONDITIONS_POSSIBLE
        if ( adjust_index(&nn_site[index]) )
        //  Returns true if index was adjusted, i.e. if index was out of bounds.
            {
            if ( free_boundary_conditions )     
                {
                nn_site[0] = NON_EXISTENT;
                //  Meaning that there is no nn in direction r.
                exists = false;
                break;
                }
            } 
        #else
        adjust_index(&nn_site[index]);
        //  With periodic boundary conditions there are no non-existent nn's
        //  so no need to check return value.
        #endif               
        }
    }

return ( exists );     
}
#endif