//  MEAS.C
//  Functions concerned with measuring quantities
//  © 2021 Peter J. Meyer

#include "iss.h"

/*-------------------------------*/
double magnetization_per_spin(void)
{
double result;

if ( model_is_ising )
    result = ( (double)abs_magnetization/num_spins );
else  //  model is q-state Potts
    result = ( ( ( q * ( 1 + (double)abs_magnetization/num_spins ) ) - 2 )
             / ( 2 * ( q - 1 ) ) );    
//  When q = 2 this expression equals the Ising magnetization value.

return ( absolute_magnetization_requested ? fabs(result) : result );
}

/*------------------------*/
void zero_measurements(void)
{
int i, qq;

for ( i=0; i<num_time_slices; i++ )
    {
    magnetization[i][0] = magnetization[i][1] = magnetization[i][3] = 0.0;
    magnetization[i][2] = NON_EXISTENT;

    if ( internal_energy_measured )
        internal_energy[i][0] = internal_energy[i][1] = 0.0;

    if ( second_moment_measured && model_is_q_state_potts )
        {
        for ( qq=2; qq<=q; qq++ )
            potts_magnetization_squares[i][qq] = 0.0;
        }

    if ( autocorrelation_measured )
        autocorrelation[i][0] = autocorrelation[i][1] = 0.0;
        //  [1] is for the standard deviation
        //  Normally not done anyway.
    }
}

/*----------------------------------------*/
void update_measurements(int time_slice_num)
{
int qq, abs_magnetization_on_entry;
double autocorr;
double magn = magnetization_per_spin();

if ( expired )
    magn = 4*fabs(magn)*(fabs(magn)-0.5)*(fabs(magn)-0.5);
 
magnetization[time_slice_num][0] += magn;
magnetization[time_slice_num][1] += magn*magn;
magnetization[time_slice_num][2] = 0;
magnetization[time_slice_num][3] += magn*magn*magn*magn;

if ( internal_energy_measured )
    {
    internal_energy[time_slice_num][0] += (double)abs_internal_energy;
    internal_energy[time_slice_num][1] 
        += (double)(abs_internal_energy) * abs_internal_energy;
    }

if ( second_moment_measured && model_is_q_state_potts )
    {
    potts_magnetization_squares[time_slice_num][1] = magnetization[time_slice_num][1];
    //  Save abs_magnetization.
    abs_magnetization_on_entry = abs_magnetization;
    for ( qq=2; qq<=q; qq++ )
        {
        abs_magnetization = get_abs_magnetization(qq);
        magn = magnetization_per_spin();
        potts_magnetization_squares[time_slice_num][qq] += magn*magn;
        }
    //  Restore abs_magnetization.
    abs_magnetization = abs_magnetization_on_entry;    
    }

if ( autocorrelation_measured )
    {
    autocorr = get_autocorrelation();
    if ( expired )
        autocorr = 4*fabs(autocorr)*(fabs(autocorr)-0.5)*(fabs(autocorr)-0.5);

    autocorrelation[time_slice_num][0] += autocorr;

    #if WR_ST_DEV_AUTOCORRELATON
    //  Sum squares of autocorrelation values to get standard deviation.
    autocorrelation[time_slice_num][1] += autocorr*autocorr;
    #endif
    }
}

//  Before the data is analysed:
//
//  magnetization[i][0] is the sum of the magnetizations per spin
//  at time slice i over all samples.
//  magnetization[i][1] is the sum of the squares of these measurements.
//  magnetization[i][2] is NON-EXISTENT if no magnetization was measured
//  for time slice i, 0 otherwise.
//  magnetization[i][3] is the sum of the fourth powers of these measurements.
//
//  If autocorrelation is measured:
//  autocorrelation[i][0] is the sum of the autocorrelation measurements
//  at time slice i over all samples.
//  autocorrelation[i][1] is the sum of the squares of these measurements.
//
//  If internal energy is measured:
//  internal_energy[i][0] is the sum of the actual internal energies
//  at time slice i over all samples.
//  internal_energy[i][1] is the sum of the squares of these.
//
//  After the data is analysed:
//
//  magnetization[i][0] is the mean magnetization per spin at time slice i over all samples.
//  magnetization[i][1] is the standard deviation of the magnetization (always calculated)
//  If the second moment of the magnetization is measured
//  magnetization[i][2] contains it.
//  If the Binder cumulant is measured then magnetization[i][3] contains it.
//
//  If autocorrelation is measured:
//  autocorrelation[i][0] is the mean autocorrelation at time slice i
//  autocorrelation[i][1] is the standard deviation of the autocorrelation
//
//  If internal energy is measured:
//  internal_energy[i][0] is the mean internal energy per spin at time slice i
//  internal_energy[i][1] is the standard deviation of the internal energy per spin
//  internal_energy[i][2] is the specific heat per spin

/*-------------------*/
void analyse_data(void)
{
int i, qq;
double mean_magnetization , mean_magnetization_squares;
double mean_internal_energy, mean_internal_energy_squares;
double second_moment, mean_autocorrelation;
#if WR_ST_DEV_AUTOCORRELATON
double mean_auto_correlation_squares;
#endif

for ( i=0; i<num_time_slices; i++ )
    {
    if ( magnetization[i][2] == NON_EXISTENT )
        break;

    //  Calculate the mean magnetization.
    mean_magnetization = magnetization[i][0] /= num_samples;

    //  Calculate the mean of the squares of the magnetization.
    mean_magnetization_squares = magnetization[i][1]/num_samples;

    //  Calculate the standard deviation of the magnetization.
    if ( num_samples == 1 )
        magnetization[i][1] = NON_EXISTENT;
    else
        magnetization[i][1] = 
            sqrt( ( mean_magnetization_squares - 
                    mean_magnetization*mean_magnetization ) 
                    / ( num_samples - 1 ) );

    if ( second_moment_measured )
        {
        //  Calculate the second moment.
        if ( model_is_ising )
            second_moment = mean_magnetization_squares;
        else    //  model is q-state Potts
            {
            second_moment = 0.0;
            for ( qq=1; qq<=q; qq++ )
                second_moment += potts_magnetization_squares[i][qq];
            second_moment /= q*num_samples;
            }
        magnetization[i][2] = second_moment;
        }

    if ( binder_cumulant_measured )
        {
        //  UL = 1 - <S^4>/(3*<S^2>^2)
        magnetization[i][3] = 
            1 - (magnetization[i][3]/num_samples)
                /(3*mean_magnetization_squares*mean_magnetization_squares);
        }

    if ( autocorrelation_measured )
        {
        //  Calculate the mean autocorrelation.
        autocorrelation[i][0] /= num_samples;
        mean_autocorrelation = autocorrelation[i][0];

        #if WR_ST_DEV_AUTOCORRELATON
        //  Calculate the mean of the squares of the autocorrelation.
        mean_auto_correlation_squares = autocorrelation[i][1]/num_samples;

        //  Calculate the standard deviation of the autocorrelation.
        if ( num_samples == 1 )
            autocorrelation[i][1] = NON_EXISTENT;
        else
            autocorrelation[i][1] = 
                sqrt( ( mean_auto_correlation_squares -
                        mean_autocorrelation*mean_autocorrelation ) 
                        / ( num_samples - 1 ) );
        #endif
        }

    if ( internal_energy_measured )
        {
        //  Calculate the mean internal energy.
        mean_internal_energy = internal_energy[i][0] /= num_samples;
        
        //  Calculate the mean of the square of the internal energy.
        mean_internal_energy_squares = internal_energy[i][1]/num_samples;

        //  Calculate the specific heat.
        internal_energy[i][2] = 
            ( mean_internal_energy_squares 
              - mean_internal_energy*mean_internal_energy )                   
              / ( temperature * temperature );

        //  Calculate the specific heat per spin
        internal_energy[i][2] /= num_spins;

        //  Calculate the mean internal energy per spin.
        //  internal_energy[i][0] /= num_spins; 

            //  Calculate the standard deviation of the internal energy.
        if ( num_samples == 1 )
            internal_energy[i][1] = NON_EXISTENT;
        else
            internal_energy[i][1] = 
                sqrt( ( mean_internal_energy_squares -
                        mean_internal_energy*mean_internal_energy )
                        / ( num_samples - 1 ) ) / num_spins;

        //  Calculate the mean internal energy per spin.
        internal_energy[i][0] /= num_spins; 
        }
    }
}