maxwell_slice_oblique

A sample to sliced Maxwell

''' A sample to sliced Maxwell
'''
import numpy as np
import numpy.random as r

import matplotlib.pyplot as plt

def theory(n, vth, vb_abs, theta):

    coeff0 = (2 * np.pi * vth * vth) ** (1.5)

    return lambda v: n / coeff0 * np.exp(-(v * v - 2 * v * vb_abs * np.cos(theta) + vb_abs ** 2) / (2 * vth * vth))
    

def oblique(vbulk, vslice, vth, n, nv=20, rv=0.1, vlim=[-10, 10]):
    # Bulk velocity
    vb = np.array(vbulk)
    vb_abs = np.sqrt((vb * vb).sum())

    # Direction of the slice vector
    v = np.array(vslice)
    lenv = np.sqrt((v * v).sum())
    v = v / lenv  # Normalized now.
    print('v = ', v, lenv)

    # Angle
    vdotvb = (v * vb).sum()
    vdotvb = vdotvb / vb_abs
    theta = np.arccos(vdotvb)
    print('angle = ', theta * 180. / np.pi)

    # Allowance of determination
    rv2 = rv * rv
    vol = (rv ** 3) * np.pi * 4. / 3.

    # First, determine the distribution
    vi = r.standard_normal([n, 3])
    vi = vi * vth   # Spread by the thermal velocity
    vi = vi + vb

    # Then, slice is simulated.  V axis prepared.
    vax = np.linspace(vlim[0], vlim[1], nv)

    # Blank array to save the results.
    f_v = np.zeros_like(vax)

    for idx, absv in enumerate(vax):
        vecv = v * absv   # A velocity vector corresponding to the absv.

        vdiff = vi - vecv   # Distance to each particles., n,3 shape
        absvdiff = (vdiff * vdiff).sum(1)  # Square distance, n shape

        # Check whether the particle is close enough to the observation point.
        n_insphere = len(np.where(absvdiff < rv2)[0])
        dens_insphere = n_insphere / vol
        f_v[idx] = dens_insphere
        print(idx, absv, f_v[idx])

    # Theory calculation.
    vax2 = np.linspace(vlim[0], vlim[1], nv*100)  # Prepare v axis
    th = theory(n, vth, vb_abs, theta)   # Prepare corresponding function

    return (vax, f_v), (vax2, th(vax2))

def main():
    vb = (1, 3, 2)
    vs = (3, 1, 4)
    vth = 1.9
    n = 2000000

    tpart, theor = oblique(vb, vs, vth, n)

    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(tpart[0], tpart[1], 'bo')
    ax.plot(theor[0], theor[1], 'b-', label='42.8 deg')

    # Second case
    vs = (1, 3, 2)  # 0 degree slice
    tpart, theor = oblique(vb, vs, vth, n)
    ax.plot(tpart[0], tpart[1], 'ro')
    ax.plot(theor[0], theor[1], 'r-', label='0.0 deg')

    # Thrd case
    vs = (4, -2, 1)  # 90 degree slice
    tpart, theor = oblique(vb, vs, vth, n)
    ax.plot(tpart[0], tpart[1], 'go')
    ax.plot(theor[0], theor[1], 'g-', label='90.0 deg')

    ax.set_ylabel('Particles per volume')
    ax.legend()

    fig.savefig('maxwell_slice_oblique.png')


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