// 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;
}
}
}