maxwell_velocity

For understanding Maxwell distribution.

Consider the Maxwell distribution with density 50 and bulk velocity 0 and thermal velocity 10 (unitless). The simplest way of generating Maxwell distribution is shown.

'''For understanding Maxwell distribution.

Consider the Maxwell distribution with density 50 and bulk velocity 0 and thermal velocity 10 (unitless).
The simplest way of generating Maxwell distribution is shown.

'''

import numpy as np
from numpy import random

import matplotlib.pyplot as plt

from irfpy.util import maxwell

def main():
    '''Main routine.
    '''
    # A simple way to get 3-D maxwell distribution.
    ntest = 1e6  # 1,000,000 particles as a sample.
    vth = 10

    v = random.standard_normal((ntest, 3))  # Test particles.
    vvec = v * vth   # 10 is the thermal velocity, or standard deviation

    vabs = np.sqrt(vvec[:, 0] ** 2 + vvec[:, 1] ** 2 + vvec[:, 2] ** 2)

    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.hist(vabs, range=[0, 50], bins=500)

    # Theoretically 3-D Maxwell distributed particle velocity should
    # follow 4 pi v^2 N (2 pi vth^2) ^ (3/2) exp(-v^2 / 2*vth^2)

    theory = lambda v: 4 * np.pi * v**2 * ntest * (2 * np.pi * vth * vth) ** (-1.5) * np.exp(-v*v / 2 / vth / vth)

    xthe = np.linspace(0, 50, 500)
    ythe = np.array([theory(x) for x in xthe])
    dx = 0.1
    ythe = ythe * dx

    ax.plot(xthe, ythe)

    ax.set_xlabel('Velocity (unitless)')
    ax.set_ylabel('Number of particle in 1 bin (dx=0.1)')
    
    fig.savefig('maxwell_velocity.png')

    return vabs


if __name__ == '__main__':
    retval = main()