From abb7005f4568c365b735d6c957039d126d3886ef Mon Sep 17 00:00:00 2001
From: LemurPwned <jakub.mojsiejuk@gmail.com>
Date: Sun, 24 Mar 2024 21:39:01 +0100
Subject: [PATCH] LLB updates

---
 core/llgb.hpp | 36 ++++++++++++++++++++++++++++++------
 1 file changed, 30 insertions(+), 6 deletions(-)

diff --git a/core/llgb.hpp b/core/llgb.hpp
index 887715d..dac9f92 100644
--- a/core/llgb.hpp
+++ b/core/llgb.hpp
@@ -135,30 +135,54 @@ class LLGBLayer
     }
 
 
-    CVector<T> nonadiabaticThermalField(T time, T timestamp) {
+    CVector<T> nonadiabaticThermalField(T time, T timestep) {
         const T temp = this->temperatureDriver.getCurrentScalarValue(time);
         const T alpha_perp2 = pow(this->alpha_perp_log, 2);
-        const T normFactor = this->volume * this->Ms / MAGNETIC_PERMEABILITY;
         const T varianceDev = (2 * BOLTZMANN_CONST * temp * (this->getAlphaPerpendicular(time)
-            - this->getAlphaParallel(time))) / (MAGNETIC_PERMEABILITY * GYRO * normFactor * alpha_perp2);
+            - this->getAlphaParallel(time))) / (this->volume * this->Ms * GYRO * alpha_perp2);
         return sqrt(varianceDev) * CVector<T>(this->distributionA);
     }
 
     CVector<T> adiabaticThermalField(T time, T  timestep) {
         const T temp = this->temperatureDriver.getCurrentScalarValue(time);
-        const T normFactor = this->volume * this->Ms / MAGNETIC_PERMEABILITY;
-        const T varianceDev = (2 * BOLTZMANN_CONST * temp * GYRO * this->getAlphaParallel(time)) / normFactor;
+        // GYRO multiplies in the stochasticTorque for consistency
+        const T varianceDev = (2 * BOLTZMANN_CONST * temp * this->getAlphaParallel(time)) / (GYRO * this->volume * this->Ms);
         return sqrt(varianceDev) * CVector<T>(this->distributionB);
     }
 
     CVector<T> stochasticTorque(const CVector<T>& m, T time, const CVector<T>& nonAdiabatic,
         const CVector<T>& adiabatic) {
+        /*
+            This formulation follows:
+            Axitia, 2015, Fundamentals and applications of the Landau–Lifshitz–Bloch equation
+            Evans, 2012, Stochastic form of the Landau-Lifshitz-Bloch equation
+            Read Evans to understand the switch.
+
+            This is more correct than stochasticTorqueOld, and used more recently
+        */
         const T inv_mlen = pow(1. / m.length(), 2);
         const T gamma_p = GYRO / (1 + pow(this->damping, 2)); // LLGS -> LL form
         const CVector<T> nonAdiabaticTerm = c_cross<T>(m, c_cross<T>(m, nonAdiabatic));
-        return -1 * gamma_p * inv_mlen * getAlphaPerpendicular(time) * nonAdiabaticTerm + adiabatic;
+        return -1 * gamma_p * inv_mlen * getAlphaPerpendicular(time) * nonAdiabaticTerm + GYRO * adiabatic;
     }
 
+    CVector<T> stochasticTorqueOld(const CVector<T>& m, T time, const CVector<T>& nonAdiabatic,
+        const CVector<T>& adiabatic) {
+        /*
+            This formulation follows:
+            Atxitia, 2007, Micromagnetic modeling of laser-induced magnetization dynamics using the Landau-Lifshitz-Bloch equation
+            And classical:
+            Garanin, 2004, Thermal fluctuations and longitudinal relaxation of single-domain magnetic particles at elevated temperatures
+        */
+        const T inv_mlen = pow(1. / m.length(), 2);
+        const T gamma_p = GYRO / (1 + pow(this->damping, 2)); // LLGS -> LL form
+        const CVector<T> nonAdiabaticTerm = c_cross<T>(m, c_cross<T>(m, nonAdiabatic));
+        const CVector<T> adiabaticTerm = c_dot(m, adiabatic) * m;
+        return gamma_p * inv_mlen * (
+            adiabaticTerm * getAlphaParallel(time) - nonAdiabaticTerm * getAlphaPerpendicular(time));
+
+
+    }
     // setters
     void setTemperatureDriver(const ScalarDriver<T>& driver)
     {