//  CORLEN.C
//  Functions for calculating correlation lengths.
//  © 2021 Peter J. Meyer

#include "iss.h"

#define PERCOLATION 0

static int num_spins_in_non_spanning_cluster;

struct _all_cluster_trace_data ct_data;

//  static int open_bond_spins[MAX_SIZE];
static int num_open_bond_spins[MAX_SIZE];
static int num_like_spins[MAX_SIZE];

static double correlation[MAX_SIZE];

static void zero_counts(void);
static void grow_cluster_and_count(int i, int j, ...);
static int num_sites_at_distance(int n);
static double calculate_corr_length(type);
static int another_spin(void);

//  corr_len_type = 0 for percolation correlation length
//                  1 for thermal correlation length
/*-----------------------------------------------*/
double get_correlation_length(int corr_length_type,  
                              int method,
                              int exclude_spanning_clusters)
{
int i, j, spin_value, n;
int num_spins_in_non_spanning_cluster=0;

zero_counts();

if ( exclude_spanning_clusters )
    {
    ct_data.cluster_type = CT_OPEN_BOND;        //  ???    
    ct_data.stop_on_spanning_cluster = false;
    trace_all_clusters(&ct_data);   //  Look for spanning cluster.
    }

switch ( dimensionality ) 
    {
    case 2:
    if ( method == 0 )
        {
        }
    else
        {
        for ( i=0; i<size; i++ )
            {
            for ( j=0; j<size; j++ )
                {
                if ( spin_value = spin_at_site(i,j) )
                    {
                    if ( exclude_spanning_clusters && ( sites2[i][j] & IN_SPANNING_CLUSTER ))
                        continue;
                    if ( !( sites2[i][j] & IN_SPANNING_CLUSTER) )
                        num_spins_in_non_spanning_cluster++;
                    grow_cluster_and_count(i,j);
                    //  Grow the cluster with seed(i,j) and count the number
                    //  of spins and like spins at all distances n,
                    //  adding to running total.
                    }
                }
            }
        }
    break;

    case 3:
    ADD_CODE        //  3d
    break;

    case 4:
    ADD_CODE        //  4d
    break;
    }


correlation[0] = 1;

for ( n=1; n<size/2; n++ )
    {
    if ( corr_length_type == PERCOLATION )
        correlation[n] = num_open_bond_spins[n]
            / ( num_sites_at_distance(n) * num_spins_in_non_spanning_cluster );
    else    //  type == THERMAL
        correlation[n] = num_like_spins[n]
            / ( num_sites_at_distance(n) * num_spins_in_non_spanning_cluster );
    }

//  Include perc_correlation[] and thermal_correlation[] in results
//  so as to allow extraction of xi by semi-log plot.

return ( calculate_corr_length(corr_length_type) );
}

//  This calculation depends on the lattice type and dimensionality.
/*-----------------------------------*/
static int num_sites_at_distance(int n)
{
return ( 1 );
}

/*-------------------------------------*/
static double calculate_corr_length(type)
{
return ( 0.0 );
}

/*-------------------------*/
static void zero_counts(void)
{
//  zero number of open bond spins
memset(&num_open_bond_spins[0],0,MAX_SIZE);

//  zero numbers of like spins
memset(&num_like_spins[0],0,MAX_SIZE);

num_spins_in_non_spanning_cluster = 0;

}

//  Grow a cluster of type T from site i,j.
//  T is WL for the whole lattice, NN for nearest neighbour,
//  OB for open-bond, SP for spin (of like values).
/*-------------------------------------------------*/
static void grow_cluster_and_count(int i, int j, ...)
{
int n;

mark_all_sites_unvisited();

MARK_SITE_AS_VISITED_2(i,j);

n = 1;

loop
    {
    while ( another_spin() )
        {
        ;
        }
    }
}

//  This returns true if there is a direction r such that
//  there is a T = WL site (occupied or not) in direction r
//                 NN spin in direction r (open or closed bond)
//                 OB spin in direction r with open bnd
//                 SP like-valued spin in direction r
//  and which is unvisited.
//  In this case the direction of the possible new spin on return is in s_r
//  and the location of the possible new spin is (nn_i,nn_j,...).
/*-------------------------*/
static int another_spin(void)
{
switch ( dimensionality )
    {
    case 2:
    break;

    case 3:
    break;

    case 4:
    break;
    }

return ( 0 );
}