maxwell_slice

A sample to check the slicing maxwelling

''' A sample to check the slicing maxwelling
'''
import numpy as np
import matplotlib.pyplot as plt
import numpy.random as r

def main():

    n = 1000000
    vth = 5

    # Make test particles.
    vall = r.standard_normal([n, 3])
    vall = vall * vth

    # Slice.  Only take |vy| <= 0.5 and |vz| <= 0.5, then dvydvz = 1
    threslice = 0.5
    vslice = []

    for i in range(n):
        if np.abs(vall[i, 1]) <= threslice and np.abs(vall[i, 2]) <= threslice:
            vslice.append(vall[i, :])

    vslice = np.array(vslice)

    # vslice[:, 0] should follow the sliced 1-D maxwellian.

    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.hist(vslice[:, 0], range=[-20, 20], bins=40)

    # Theoretical expectation

    xs = np.linspace(-20, 20, 400)
    theory1 = lambda v: n * (2 * np.pi * vth * vth) ** (-3./2.) * np.exp(-v * v / 2 / vth / vth)
    ys = [theory1(ix) for ix in xs]
    ax.plot(xs, ys)
    print(max(ys))

    ax.set_title('3D maxwell sliced to 1D')
    ax.set_xlabel('v [m/s]')
    ax.set_ylabel('f(v) [s^3/m^6]')

    fig.savefig('maxwell_slice.png')

    return vall, vslice

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