putting everything in 1 place

maxwellian_kicks_python/maxwell.py

+import math
+import matplotlib.pyplot as plt
+import numpy as np
+
+from scipy.special import erf
+
+import random
+
+
+# Hnn
+# http://pulsars.info/Thesis.pdf 
+# https://arxiv.org/pdf/1807.11489.pdf 
+# http://mathworld.wolfram.com/MaxwellDistribution.html
+# Give different formulae
+
+# Especially the first factor is a bit unclear sometimes
+
+def maxwell_prob_density(vkick, dispersion):
+    return np.sqrt(2.0/np.pi)* (vkick**2)/(dispersion**3) * np.exp(-((vkick**2)/(2*(dispersion**2)))) 
+
+def maxwell_cum_density(vkick, dispersion):
+    return erf((vkick)/(np.sqrt(2)*dispersion)) - ((vkick * np.exp(-((vkick**2)/(2*(dispersion**2)))))/(dispersion))*(np.sqrt(2.0/math.pi))
+
+
+# sampling from distribution:
+def sample_maxwell(sigma):
+    X = random.random()
+    Y = random.random()
+    W = random.random()
+    Z = random.random()
+
+    s = sigma*math.sqrt(-2*math.log(1-X))
+    p = sigma*math.sqrt(-2*math.log(1-Y))
+
+    theta = 2*math.pi*Y
+    phi = 2*math.pi*Z
+
+    u = s * math.cos(theta)
+    v = s * math.sin(theta)
+    w = p * math.cos(phi)
+
+    kick = math.sqrt((u**2) + (v**2) + (w**2))
+    return kick
+
+
+dispersion = 265
+
+samples = 50000
+sampled_array = [sample_maxwell(dispersion) for _ in range(samples)]
+
+v_arr = np.arange(0, 1500, 1)
+mw_prob = maxwell_prob_density(v_arr, dispersion)
+mw_cum = maxwell_cum_density(v_arr, dispersion)
+
+fig, axes = plt.subplots(nrows=1, ncols=2)
+axes[0].plot(v_arr, mw_prob)
+axes[0].hist(sampled_array, density=True, bins=100)
+axes[0].set_ylabel('Probability density')
+axes[0].set_xlabel('v (km s-1)')
+
+axes[1].plot(v_arr, mw_cum)
+axes[1].set_ylabel('Cumulative density')
+axes[1].set_xlabel('v (km s-1)')
+
+plt.show()