# 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.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()