Add `classical_spin_temperature` function

[Lazersmoke]

Re: https://www.sciencedirect.com/science/article/pii/S0010465516301412

We should steal the equation for "spin temperature" used in that paper (and derived in [26]). Here's a reference implementation:

function classical_spin_temperature(sys)
  dHdS = -Sunny.energy_grad_dipoles(sys)
  numer = sum([norm(cross(sys.dipoles[i],dHdS[i]))^2 for i = eachindex(sys.dipoles)])
  denom = 2 * sum([dot(sys.dipoles[i],dHdS[i]) for i = eachindex(sys.dipoles)])
  numer / denom
end

Using this, I annealed a system (with exponentially decreasing temperature schedule) and got the attached plot of "spin temperature" (y axis) vs Langevin temperature (x axis) [both log scale]. As you can see, the spin temperature formula accurately tracks the annealing temperature [26] https://journals.aps.org/pre/pdf/10.1103/PhysRevE.82.031111

These two annealing runs were done insufficiently slowly, and you can see that the spin temperature doesn't track the thermostat temperature, reflecting this!

So this is a useful tool for people to check their annealing and also would be good for an example in the docs

[kbarros]

This is a nice idea. What should the behavior be for SU(N) mode?