Hi,
I believe the issue is when you set the magnetization. In (and also the following 6 lines)
m. setRegion(1, Uniform(sin(pi/36)*RandNorm(), sin(pi/36)*sqrt(1-(RandNorm()*RandNorm())), cos(pi/36)))
RandNorm() is not guaranteed to be < 1 in magnitude. The definition is a
normally distributed float64 in the range [-math.MaxFloat64, +[math.MaxFloat64]] with standard normal distribution (mean = 0, stddev = 1) (source). This will tend to be less than one, given the mean and standard definition, but some tails of the distribution extend beyond +/-1 . So you can end up taking the square root of a negative number. You can instead use Rand() instead of RandNorm(), which is between 0 and 1 only.
Best,
Josh L.