//  INITSPIN.C
//  Functions for initialization of spins.
//  © 2021 Peter J. Meyer

#include "iss.h"

#define ISING_UP_SPIN_PROBABILITY \
    ((num_spins*(1+initial_magnetization)-2*num_up_spins) \
    /(2*(num_spins-num_up_spins-num_down_spins)))

#define Q_STATE_POTTS_UP_SPIN_PROBABILITY \
    (((double)num_spins/q)*((q-1)*initial_magnetization+1)-num_up_spins) \
    /( num_spins-num_up_spins-num_other_spins)

static void assign_initial_ising_spins(void);
static void assign_initial_q_state_potts_spins(void);

//  On return from this function abs_magnetization 
//  and abs_internal_energy have been calculated.
/*---------------------------*/
void assign_initial_spins(void)
{
if ( model_is_ising )
    assign_initial_ising_spins();
else    //  model is q-state Potts
    assign_initial_q_state_potts_spins();
}

//  Call with initial_magnetization = 1.0 for all spins up,
//  -1.0 for all spins down, 0.0 for random up- and down-spins.
//  This function does not assume any particular state of the spins.
//  It also calculates the absolute magnetization.
/*----------------------------------------*/
static void assign_initial_ising_spins(void)
{
int i, j, k, l;
int num_up_spins=0;
int num_down_spins=0;

switch ( dimensionality )
    {
    case 2:
    for ( i=0; i<size; i++ )
        {
        for ( j=0; j<size; j++ )
            {
            if ( site_diluted )     //  site_diluted set in SETPARAM.C.
                if ( !SITE_IS_OCCUPIED_2(i,j) )
                    continue;

            //  Site is occupied.
            if ( initial_magnetization == 1.0 )
                {
                sites2[i][j] &= BIT0_0;     //  Set spin to up.
                num_up_spins++;
                }
            else if ( initial_magnetization == -1.0 )
                {
                sites2[i][j] |= BIT1_0;     //  Set spin to down.
                num_down_spins++;
                }
            else
                {
                //  Under some (rather unlikely) scenarios ISING_UP_SPIN_PROBABILITY
                //  could be < 0 or > 1, but this is not a problem.
                if ( PRNG1() < ISING_UP_SPIN_PROBABILITY )
                    {
                    sites2[i][j] &= BIT0_0;     //  Set spin to up.
                    num_up_spins++;
                    }
                else
                    {
                    sites2[i][j] |= BIT1_0;     //  Set spin to down.
                    num_down_spins++;
                    }
                }
            }
        }
    break;

    case 3:
    for ( i=0; i<size; i++ )
        {
        for ( j=0; j<size; j++ )
            {
            for ( k=0; k<size; k++ )
                {
                if ( site_diluted )
                    if ( !SITE_IS_OCCUPIED_3(i,j,k) )
                        continue;

                //  Site is occupied.
                if ( initial_magnetization == 1.0 )
                    {
                    sites3[i][j][k] &= BIT0_0;     //  Set spin to up.
                    num_up_spins++;
                    }
                else if ( initial_magnetization == -1.0 )
                    {
                    sites3[i][j][k] |= BIT1_0;     //  Set spin to down.
                    num_down_spins++;
                    }
                else
                    {
                    if ( PRNG1() < ISING_UP_SPIN_PROBABILITY )
                        {
                        sites3[i][j][k] &= BIT0_0;     //  Set spin to up.
                        num_up_spins++;
                        }
                    else
                        {
                        sites3[i][j][k] |= BIT1_0;     //  Set spin to down.
                        num_down_spins++;
                        }
                    }
                }
            }
        }

    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 ( site_diluted )
                        if ( !SITE_IS_OCCUPIED_4(i,j,k,l) )
                            continue;

                    //  Site is occupied.
                    if ( initial_magnetization == 1.0 )
                        {
                        sites4[i][j][k][l] &= BIT0_0;     //  Set spin to up.
                        num_up_spins++;
                        }
                    else if ( initial_magnetization == -1.0 )
                        {
                        sites4[i][j][k][l] |= BIT1_0;     //  Set spin to down.
                        num_down_spins++;
                        }
                    else
                        {
                        if ( PRNG1() < ISING_UP_SPIN_PROBABILITY )
                            {
                            sites4[i][j][k][l] &= BIT0_0;     //  Set spin to up.
                            num_up_spins++;
                            }
                        else
                            {
                            sites4[i][j][k][l] |= BIT1_0;     //  Set spin to down.
                            num_down_spins++;
                            }
                        }
                    }
                }
            }
        }
    break;
    }

if ( adjust_zero_initial_magnetization )
    if ( num_up_spins == num_down_spins )
        {
        flip_initial_spin_to_up();
        num_up_spins++;
        num_down_spins--;
        }

abs_magnetization = num_up_spins - num_down_spins;

if ( internal_energy_measured )
    //  Since the bonds have already been established we can now:
    abs_internal_energy = get_abs_internal_energy();
}

//  Call with initial_magnetization = 1.0 for all spins up,
//  -1/(q-1) for all spins in another direction, 0.0 for random spins.
//  This function does not assume any particular state of the spins.
//  It returns the number of up spins minus the number of other spins.
/*------------------------------------------------*/
static void assign_initial_q_state_potts_spins(void)
{
int i, j, k, l;
int spin_value;
int num_up_spins=0, num_other_spins=0;

switch ( dimensionality )
    {
    case 2:
    for ( i=0; i<size; i++ )
        {
        for ( j=0; j<size; j++ )
            {
            if ( site_diluted )
                if ( !SITE_IS_OCCUPIED_2(i,j) )
                    continue;

            //  Site is occupied.
            sites2[i][j] &= ~spin_value_mask;   //  Set spin to up.
            if ( initial_magnetization == 1.0 )
                //  Spin is up.
                num_up_spins++;
            else if ( initial_magnetization == -1.0/(q-1))
                {
                //  Set spin to another direction.
                spin_value = 2 + (int)(PRNG1()*(q-1));
                sites2[i][j] |= --spin_value;
                num_other_spins++;
                }
            else if ( PRNG1() < Q_STATE_POTTS_UP_SPIN_PROBABILITY )
                //  Spin is up.
                num_up_spins++;
            else
                {
                //  Set spin to another direction.
                spin_value = 2 + (int)(PRNG1()*(q-1));
                sites2[i][j] |= --spin_value;
                num_other_spins++;
                }
            }
        }
    break;

    case 3:     
    for ( i=0; i<size; i++ )
        {
        for ( j=0; j<size; j++ )
            {
            for ( k=0; k<size; k++ )
                {
                if ( site_diluted )
                    if ( !SITE_IS_OCCUPIED_3(i,j,k) )
                        continue;

                //  Site is occupied.
                sites3[i][j][k] &= ~spin_value_mask;   //  Set spin to up.
                if ( initial_magnetization == 1.0 )
                    //  Spin is up.
                    num_up_spins++;
                else if ( initial_magnetization == -1.0/(q-1))
                    {
                    //  Set spin to another direction.
                    spin_value = 2 + (int)(PRNG1()*(q-1));
                    sites3[i][j][k] |= --spin_value;
                    num_other_spins++;
                    }
                else if ( PRNG1() < Q_STATE_POTTS_UP_SPIN_PROBABILITY )
                    //  Spin is up.
                    num_up_spins++;
                else
                    {
                    //  Set spin to another direction.
                    spin_value = 2 + (int)(PRNG1()*(q-1));
                    sites3[i][j][k] |= --spin_value;
                    num_other_spins++;
                    }
                }
            }
        }
    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 ( site_diluted )
                        if ( !SITE_IS_OCCUPIED_4(i,j,k,l) )
                            continue;

                    //  Site is occupied.
                    sites4[i][j][k][l] &= ~spin_value_mask;   //  Set spin to up.
                    if ( initial_magnetization == 1.0 )
                        //  Spin is up.
                        num_up_spins++;
                    else if ( initial_magnetization == -1.0/(q-1))
                        {
                        //  Set spin to another direction.
                        spin_value = 2 + (int)(PRNG1()*(q-1));
                        sites4[i][j][k][l] |= --spin_value;
                        num_other_spins++;
                        }
                    else if ( PRNG1() < Q_STATE_POTTS_UP_SPIN_PROBABILITY )
                        //  Spin is up.
                        num_up_spins++;
                    else
                        {
                        //  Set spin to another direction.
                        spin_value = 2 + (int)(PRNG1()*(q-1));
                        sites4[i][j][k][l] |= --spin_value;
                        num_other_spins++;
                        }
                    }
                }
            }
        }
    break;
    }

if ( adjust_zero_initial_magnetization )
    if ( num_other_spins - num_up_spins ==  num_spins/2 )
        {
        flip_initial_spin_to_up();
        num_up_spins++;
        num_other_spins--;
        }

abs_magnetization = num_up_spins - num_other_spins;  

if ( internal_energy_measured )
    //  Since the bonds have already been established we can now:
    abs_internal_energy = get_abs_internal_energy();
}