//  SPINFLIP.C
//  Functions concerned with spin flips
//  © 2021 Peter J. Meyer

#include "iss.h"

//  Returns true if spin to be flipped, false otherwise.
//  Ising: s is 2 * spin * spin_sum
//  Potts: s is the sum of the kronecker deltas for the initial spin value 
//  minus the sum for the new spin value.
//  This function is called only with the Metropolis or the Glauber algorithm.
/*----------------*/
int flip_spin(int s)
{
return (  dynamics_is_metropolis && s <= 0 ? true : PRNG1() < w[s+2*coord_num] );
}

/*------------------------------------------------------*/
void do_ising_spin_flip(int spin_value, int i, int j, ...)
{
int k, l;

num_spins_flipped_this_sweep++;

abs_magnetization -= 2*spin_value;

switch ( dimensionality )
    {
    case 2:
    if ( internal_energy_measured )
        abs_internal_energy += 2*spin_value*nn_spin_sum(i,j);
        //  Actually abs_internal_energy += 2*J*spin_value*nn_spin_sum(i,j);
        //  but J factored in in norm_internal_energy() in MEAS.C.

    switch ( spin_value )
        {
        case UP_SPIN:        //  Bit 0 of sites2[i][j] is 0.
        sites2[i][j] |= BIT1_0;
        break;

        case DOWN_SPIN:     //  Bit 0 of sites2[i][j] is 1.
        sites2[i][j] &= BIT0_0;
        break;
        }
    break;

    case 3:
    k = *(((int *)(&j))+1);
    if ( internal_energy_measured )
        abs_internal_energy += 2*spin_value*nn_spin_sum(i,j,k);

    switch ( spin_value )
        {
        case UP_SPIN:        
        sites3[i][j][k] |= BIT1_0;
        break;

        case DOWN_SPIN:    
        sites3[i][j][k] &= BIT0_0;
        break;
        }
    break;

    case 4:
    k = *(((int *)(&j))+1);
    l = *(((int *)(&j))+2);
    if ( internal_energy_measured )
        abs_internal_energy += 2*spin_value*nn_spin_sum(i,j,k,l);

    switch ( spin_value )
        {
        case UP_SPIN:        
        sites4[i][j][k][l] |= BIT1_0;
        break;

        case DOWN_SPIN:     
        sites4[i][j][k][l] &= BIT0_0;
        break;
        }
    break;
    }
}

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

num_spins_flipped_this_sweep++;

if ( spin_value == UP_SPIN )
    abs_magnetization -= 2;
    //  Because spin changes from an up-spin to a non-up-spin.
else if ( new_spin_value == UP_SPIN )
    abs_magnetization += 2;     
    //  Because spin changes from a non-up-spin to an up-spin.

//  If neither the spin value nor the new spin value are up-spins
//  then the absolute magnetization does not change.

switch ( dimensionality )
    {
    case 2:
    //  TENTATIVE !!
    //  Need to check calculation of internal energy and specific heat with Potts model.
    if ( internal_energy_measured )    
        abs_internal_energy += 2*spin_value*sum_kronecker_deltas(spin_value,i,j);

    sites2[i][j] &= ~spin_value_mask;   
    sites2[i][j] |= --new_spin_value;
    break;

    case 3:
    k = *(((int *)(&j))+1);
    if ( internal_energy_measured )
        abs_internal_energy += 2*spin_value*sum_kronecker_deltas(spin_value,i,j,k);

    sites3[i][j][k] &= ~spin_value_mask;   
    sites3[i][j][k] |= --new_spin_value;
    break;

    case 4:
    k = *(((int *)(&j))+1);
    l = *(((int *)(&j))+2);
    if ( internal_energy_measured )
        abs_internal_energy += 2*spin_value*sum_kronecker_deltas(spin_value,i,j,k,l);

    sites4[i][j][k][l] &= ~spin_value_mask;   
    sites4[i][j][k][l] |= --new_spin_value;
    break;
    }
}

//  The following function is called from assign_initial_spins()
//  when an assignment of initial spins results in zero magnetization 
//  and the user has specified adjust_zero_initial_magnetization.
//  This function calls functions above which change the value of abs_magnetization.
//  This is of no import since the return value of assign_initial_spins() 
//  is assigned to abs_magnetization.
/*------------------------------*/
void flip_initial_spin_to_up(void)
{
int i, j, k, l, spin, site_occupied;

//  Locate a non-UP-spin and flip it to UP.
switch ( dimensionality )
    {
    case 2:
    for ( i=0; i<size; i++ )
        {
        for ( j=0; j<size; j++ )
            {
            spin = site_occupied = spin_at_site(i,j);
            if ( site_occupied )
                {
                if ( spin != UP_SPIN )
                    {
                    if ( model_is_ising )
                        do_ising_spin_flip(spin,i,j);
                    else //  model is q-state Potts
                        do_q_state_potts_spin_flip(spin,UP_SPIN,i,j);
                    return;
                    }
                }
            }
        }
    break;

    case 3:
    for ( i=0; i<size; i++ )
        {
        for ( j=0; j<size; j++ )
            {
            for ( k=0; k<size; k++ )
                {
                spin = site_occupied = spin_at_site(i,j,k);
                if ( site_occupied )
                    {
                    if ( spin != UP_SPIN )
                        {
                        if ( model_is_ising )
                            do_ising_spin_flip(spin,i,j,k);
                        else //  model is q-state Potts
                            do_q_state_potts_spin_flip(spin,UP_SPIN,i,j,k);
                        return;
                        }
                    }
                }
            }
        }
    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++ )
                    {
                    spin = site_occupied = spin_at_site(i,j,k,l);
                    if ( site_occupied )
                        {
                        if ( spin != UP_SPIN )
                            {
                            if ( model_is_ising )
                                do_ising_spin_flip(spin,i,j,k,l);
                            else //  model is q-state Potts
                                do_q_state_potts_spin_flip(spin,UP_SPIN,i,j,k,l);
                            return;
                            }
                        }
                    }
                }
            }
        }
    break;
    }
}