//  SPINS.C
//  Functions concerned with spins
//  © 2021 Peter J. Meyer

#include "iss.h"

//  Given that a site is occupied, what is its spin?
//  (Use this if a site is known to be occupied
//  and the spin value is required.)
//
//  For the Ising model, bit 0 of sites[...] has 
//  0 for UP_SPIN and 1 for DOWN_SPIN.
//  For the q-state Potts model, the five lowest-order bits 
//  hold the spin value minus one.
//  q-state Potts model spin values are in the range 1 through q.
//
//  This is very  similar to spin_at_site().
/*----------------------------------------*/
int spin_at_occupied_site(int i, int j, ...)
{
//  int k, l;
unsigned char site;

if ( i < 0 )                //  Occurs only when i == NON_EXISTENT.
    return ( false );       

switch ( dimensionality )
    {
    case 2:
    site = sites2[i][j];
    break;

    case 3:
    //  k = *(((int *)(&j))+1);
    //  site = sites3[i][j][k];
    site = sites3[i][j][*(((int *)(&j))+1)];
    break;

    case 4:
    //  k = *(((int *)(&j))+1);
    //  l = *(((int *)(&j))+2);
    //  site = sites4[i][j][k][l];
    site = sites4[i][j][*(((int *)(&j))+1)][*(((int *)(&j))+2)];
    break;
    }

if ( model_is_ising )
    return ( site & spin_value_mask ? DOWN_SPIN : UP_SPIN );
else //  model_is_q_state_potts
    return ( ( site & spin_value_mask ) + 1 );
}


//  This function returns 0 if a site is not occupied
//  otherwise it returns the spin value (Ising or q-state Potts).
//  It should be used when both occupancy and spin value
//  information is needed.
/*-------------------------------*/
int spin_at_site(int i, int j, ...)
{
//  int k, l;
unsigned char site;

if ( i < 0 )                //  Occurs only when i == NON_EXISTENT.
    return ( 0 );       

switch ( dimensionality )
    {
    case 2:
    site = sites2[i][j];
    break;

    case 3:
    //  k = *(((int *)(&j))+1);
    //  site = sites3[i][j][k];
    site = sites3[i][j][*(((int *)(&j))+1)];
    break;

    case 4:
    //  k = *(((int *)(&j))+1);
    //  l = *(((int *)(&j))+2);
    //  site = sites4[i][j][k][l];
    site = sites4[i][j][*(((int *)(&j))+1)][*(((int *)(&j))+2)];
    break;
    }


if ( ! ( site & OCCUPIED ) )
    return ( 0 );

//  Site is occupied.
if ( model_is_ising )
    return ( site & spin_value_mask ? DOWN_SPIN : UP_SPIN );
else //  model_is_q_state_potts
    return ( ( site & spin_value_mask ) + 1 );
}

//  What is the sum of the nn spins of a spin?
//
//  This information is obtained from the bond information for the site.
//  Note that an occupied nn site is not a nn spin if no open bond exists.
//
//  This is used only with the Ising model.
/*------------------------------*/
int nn_spin_sum(int i, int j, ...)
{
int k, l, r, spin_sum=0;

switch ( dimensionality )
    {
    case 2:
    get_nn_sites(i,j);
    for ( r=0; r<coord_num; r++ )
        {
        if ( bonds2[i][j] & bit[r] )   //  if bond exists in direction r
            spin_sum += spin_at_occupied_site(nn_sites[r][0],nn_sites[r][1]);
        }
    break;

    case 3:
    k = *(((int *)(&j))+1);
    get_nn_sites(i,j,k);
    for ( r=0; r<coord_num; r++ )
        {
        if ( bonds3[i][j][k] & bit[r] )   //  if bond exists in direction r
            spin_sum += spin_at_occupied_site(nn_sites[r][0],
                nn_sites[r][1],nn_sites[r][2]);
        }
    break;

    case 4:
    k = *(((int *)(&j))+1);
    l = *(((int *)(&j))+2);
    get_nn_sites(i,j,k,l);
    for ( r=0; r<coord_num; r++ )
        {
        if ( bonds4[i][j][k][l] & bit[r] )   //  if bond exists in direction r
            spin_sum += spin_at_occupied_site(nn_sites[r][0],
                nn_sites[r][1],nn_sites[r][2],nn_sites[r][3]);
        }
    break;
    }

return ( spin_sum );
}

//  What is the sum of the kronecker deltas 
//  of a spin with its nearest neighbour spins?
//
//  This information is obtained from the bond information for the site.
//  Note that an occupied nn site is not a nn spin if no open bond exists.
//
//  This is used only with the q-state Potts model.
/*------------------------------------*/
int sum_kronecker_deltas(int spin_value, 
                         int i, int j, ...)
{
int k, l, r, sum=0;

switch ( dimensionality )
    {
    case 2:
    get_nn_sites(i,j);
    for ( r=0; r<coord_num; r++ )
        {
        if ( bonds2[i][j] & bit[r] )   //  if bond exists in direction r
            sum += ( spin_value == 
                spin_at_occupied_site(nn_sites[r][0],nn_sites[r][1]) );
        }
    break;

    case 3:
    k = *(((int *)(&j))+1);
    get_nn_sites(i,j,k);
    for ( r=0; r<coord_num; r++ )
        {
        if ( bonds3[i][j][k] & bit[r] )   //  if bond exists in direction r
           sum += ( spin_value == 
                spin_at_occupied_site(nn_sites[r][0],nn_sites[r][1],nn_sites[r][2]) );
        }
    break;

    case 4:
    k = *(((int *)(&j))+1);
    l = *(((int *)(&j))+2);
    get_nn_sites(i,j,k,l);
    for ( r=0; r<coord_num; r++ )
        {
        if ( bonds4[i][j][k][l] & bit[r] )   //  if bond exists in direction r
           sum += ( spin_value == 
                spin_at_occupied_site(nn_sites[r][0],nn_sites[r][1],
                    nn_sites[r][2],nn_sites[r][3]) );
        }
    break;
    }

return ( sum );
}

//  This calculates the absolute magnetization from scratch.
//  distinguished_spin_value is normally 1, 
//  but when calculating the second moment of the magnetization
//  in the q-state Potts model we call this function 
//  with distinguished_spin_value = 1, ..., q. 
/*---------------------------------------------------*/
int get_abs_magnetization(int distinguished_spin_value)
{
int i, j, k, l;                     //  array indices 
int abs_mag=0;
int spin_value;

if ( model_is_ising && distinguished_spin_value != UP_SPIN )
    {
    printf("\nModel is Ising, but get_abs_magnetization() called"
           "\nwith distinguished_spin_value != UP_SPIN.\n");
    exit(1);
    }
else if ( model_is_q_state_potts && 
          ( distinguished_spin_value < 1 || distinguished_spin_value > q ) )
    {
    printf("\nModel is q-state Potts and get_abs_magnetization()"
           "\ncalled with invalid value for distinguished_spin_value.\n");
    exit(1);
    }

switch ( dimensionality )
    {
    case 2: 
    for ( i=0; i<size; i++ )
        {
        for ( j=0; j<size; j++ )
            {
            if ( spin_value = spin_at_site(i,j) )
                //  spin_value is 0 if site is not occupied.
                abs_mag += ( spin_value == distinguished_spin_value ? 1 : -1 );
            }
        }
    break;

    case 3:
    for ( i=0; i<size; i++ )
        {
        for ( j=0; j<size; j++ )
            {
            for ( k=0; k<size; k++ )
                {
                if ( spin_value = spin_at_site(i,j,k) )
                    abs_mag += ( spin_value == distinguished_spin_value ? 1 : -1 );
                }
            }
        }
    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 ( spin_value = spin_at_site(i,j,k,l) )
                        abs_mag += 
                        ( spin_value == distinguished_spin_value ? 1 : -1 );
                    }
                }
            }
        }                    
    }

return ( abs_mag );
}

//  abs_internal_energy is the sum of the products of the spins
//  in each pair of nearest neighbour spins.
//  This function presupposes that nn bonds have been set up.
/*-----------------------------*/
int get_abs_internal_energy(void)
{
int i,j;
//  int k,l;
int spin_value,abs_int_energy=0;

if ( model_is_q_state_potts )
    {
    printf("\nMeasurement of internal energy for Potts model"
           "\nis not currently supported.  See SPINS.C and SPINFLIP.C.\n");
    exit(1);
    }

switch ( dimensionality )
    {
    case 2:
    for ( i=0; i<size; i++ )
        {
        for ( j=0; j<size; j++ )
            {
             if ( spin_value = spin_at_site(i,j) )
                //  This returns 0 if the site is unoccupied.
                {
                if ( model_is_ising )
                    abs_int_energy -= spin_value*nn_spin_sum(i,j);
                else    //  model is q-state Potts
                    //  TENTATIVE!!
                    abs_int_energy -= spin_value*sum_kronecker_deltas(spin_value,i,j);
                }
            }
        }
    break;
    
    case 3:
    ADD_CODE    //  3d
    break;

    case 4:
    ADD_CODE    //  4d
    break;
    }

return ( abs_int_energy/2 );
//  Divide by 2 since each energy between pairs of spins is counted twice.
}

//  For measurement of autocorrelation
//  we make a copy of the spin assignment.
//  initial_sites[..] also used to restore spin assignment
//  when num_repetitions > 0
/*---------------------------*/
void copy_spin_assignment(void)
{
switch ( dimensionality )
    {
    case 2:
    memcpy(&initial_sites2[0][0], &sites2[0][0], num_sites);
    break;

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

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

/*------------------------------*/
void restore_spin_assignment(void)
{
switch ( dimensionality )
    {
    case 2:
    memcpy(&sites2[0][0], &initial_sites2[0][0], num_sites);
    break;

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

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

//  This is the same as spin_at_occupied_site() except that
//  initial_sites[...] is used rather than sites[...].
/*------------------------------------------------*/
int initial_spin_at_occupied_site(int i, int j, ...)
{
int k, l;
unsigned char site;

switch ( dimensionality )
    {
    case 2:
    site = initial_sites2[i][j];
    break;

    case 3:
    k = *(((int *)(&j))+1);
    site = initial_sites3[i][j][k];
    break;

    case 4:
    k = *(((int *)(&j))+1);
    l = *(((int *)(&j))+2);
    site = initial_sites4[i][j][k][l];
    break;
    }

if ( model_is_ising )
    return ( site & spin_value_mask ? DOWN_SPIN : UP_SPIN );
else //  model_is_q_state_potts
    return ( ( site & spin_value_mask ) + 1 );
}