import pylab

filename = 'Co-FCC.hdf5'

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

exchange_stiffness = heisenberg_exchange.exchangeStiffness()

# Calculate the extrapolated exchange stiffness using default values.
extrapolated_exchange_stiffness, exchange_stiffness_vec, damping_factor_range, p_fit \
    = heisenberg_exchange.extrapolatedExchangeStiffness()

# Print the results
nlprint('-------------------------------------------------')
nlprint('Exchange stiffness:')
nlprint('-------------------------------------------------')
nlprint('A_ex                         : {: .1f}'.format(exchange_stiffness))
nlprint('A_ex (extrapolated)          : {: .1f}'.format(extrapolated_exchange_stiffness))
nlprint('-------------------------------------------------')

# Plot the data and fit
damping_factor_range_fit = numpy.linspace(0, damping_factor_range[-1].inUnitsOf(Ang**-1), 100)

pylab.figure()
pylab.plot(0, extrapolated_exchange_stiffness.inUnitsOf(meV*Ang**2), 'ro', label='Extrapolated result')
pylab.plot(damping_factor_range.inUnitsOf(Ang**-1), exchange_stiffness_vec.inUnitsOf(meV*Ang**2), 'ko', label='Data')
pylab.plot(damping_factor_range_fit, numpy.polyval(p_fit, damping_factor_range_fit), 'k-', label='Fit')
pylab.xlabel('Damping factor  ($1/\\AA$)')
pylab.ylabel('Exchange stiffness (mev $\\AA^2$)')
pylab.legend(loc=0)
pylab.show()