//  AUTOCORR.C
//  Functions for calculating autocorrelation
//  © 2021 Peter J. Meyer

#include "iss.h"

/*----------------------------*/
double get_autocorrelation(void)
{
int i, j, k, l;         //  array indices 
int spin_value, autocorrelation_sum=0;

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.
                autocorrelation_sum += 
                    ( spin_value == initial_spin_at_occupied_site(i,j) );
            }
        }
    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) )
                //  spin_value is 0 if site is not occupied.
                    autocorrelation_sum += 
                    ( spin_value == initial_spin_at_occupied_site(i,j,k) );
                }
            }
        }
    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) )
                    //  spin_value is 0 if site is not occupied.
                        autocorrelation_sum += 
                        ( spin_value == initial_spin_at_occupied_site(i,j,k,l) ); 
                    }
                }
            }
        }
    break;
    }

if ( model_is_ising )
    return ( (2*(double)autocorrelation_sum)/num_spins - 1 );
    //  maximum = 1, minimum = -1, no correlation = 0
else    //  model is q-state Potts
    return ( (double)autocorrelation_sum/num_spins - 1.0/q);
    //  maximum = (q-1)/q, minimum = -1/q, no correlation = 0
}