import pylab

filename = 'Co-FCC.hdf5'

# Load the HeisenbergExchange object
heisenberg_exchange = nlread(filename, HeisenbergExchange)[0]

# Calculate the Curie temperature within the mean field approximation (MFA).
T_C_mfa-tutorial = heisenberg_exchange.curieTemperature()

# Calculate the Curie temperature within the random phase approximation (RPA)
# using default q-grid (not converged).
T_C_rpa-tutorial = heisenberg_exchange.curieTemperatureRPA()

# Setup a list with number of q-points in each direction.
number_of_qpoints = numpy.array([10, 20, 30, 40])
T_C_rpa-tutorial_q_grid = []

# Calculate the Curie temperature for increasing grid sizes.
for nq in number_of_qpoints:
    q_grid = MonkhorstPackGrid(nq, nq, nq)
    T_C_rpa-tutorial_q_grid.append(heisenberg_exchange.curieTemperatureRPA(q_grid=q_grid).inUnitsOf(Kelvin))

# Make a linear fit to 1/nq vs. T_C
fit = numpy.polyfit(1/number_of_qpoints, T_C_rpa-tutorial_q_grid, 1)

# Evaluate the fit at 1/nq = 0 corresponding to infinite number of q-points.
T_C_rpa-tutorial_extrapolate = numpy.polyval(fit, 0) * Kelvin

# Print the results
nlprint('-------------------------------------------------')
nlprint('Curie temperature:')
nlprint('-------------------------------------------------')
nlprint('MFA                                  : {: .1f}'.format(T_C_mfa-tutorial))
nlprint('RPA (default q-grid)                 : {: .1f}'.format(T_C_rpa-tutorial))
nlprint('RPA (extrapolate to infinite q-grid) : {: .1f}'.format(T_C_rpa-tutorial_extrapolate))
nlprint('-------------------------------------------------')

# Plot the data and linear fit
pylab.figure()
pylab.plot(1/number_of_qpoints, T_C_rpa-tutorial_q_grid, 'ro', label='Data')
pylab.plot([0, 0.12], numpy.polyval(fit, [0, 0.12]), 'k-', label='Fit')
pylab.xlabel('$1/n_q$')
pylab.ylabel('Curie temperature (K)')
pylab.legend(loc=0)
pylab.show()