# maxwell_dfluxΒΆ

Simulate differntial energy flux

''' Simulate differntial energy flux
'''
import numpy as np
from numpy import random as r

import matplotlib.pyplot as plt

def main_prog():
'''
'''

# Test particles

nall = 10000000  # 1e7

# First determine the position of the particles.

xyz_all = r.rand(nall, 3) * 2 - 1
r2 = (xyz_all ** 2).sum(1)

# Only valid is inside a unit sphere
xyz = xyz_all[np.where(r2<=1)[0], :]
x = xyz[:, 0]
y = xyz[:, 1]
z = xyz[:, 2]

n_valid = xyz.shape[0]

print('Particles = %d' % n_valid)

# Each particle should have velocity following Maxwell distribution.
vth = 3

vxyz = r.standard_normal([n_valid, 3]) * 3
vx = vxyz[:, 0]
vy = vxyz[:, 1]
vz = vxyz[:, 2]

# Put detecting square in the z-plane, size is 2*delta in x- and y-directions
delta = 0.1

# Cross section of the particle in the z-plane.
k = - z / vz
x0 = x - k * vx
y0 = y - k * vy

# Check if it is detectd.
kposi = k > 0
x0_l = x0 > -delta
x0_h = x0 < delta
y0_l = y0 > -delta
y0_h = y0 < delta

detection = np.array([kposi, x0_l, x0_h, y0_l, y0_h]).all(0)

print('N_detected = %d' % len(x[detection]))

# Particle detected

xyz = xyz[detection, :]
x, y, z = xyz[:, 0], xyz[:, 1], xyz[:, 2]
vxyz = vxyz[detection, :]
vx, vy, vz = vxyz[:, 0], vxyz[:, 1], vxyz[:, 2]
E = (vxyz ** 2).sum(1)   # Assuming m=2
v = np.sqrt(E)
print(vxyz[0, :], E[0], v[0])
normv = (vxyz.T / v).T   # To broadcase, T is used.  Element-wise normalization.
print(normv[0, :])
theta = np.arccos(normv[:, 2]) * 180 / np.pi
print(theta[0], theta.shape)
phi = np.arctan2(normv[:, 1], normv[:, 0]) * 180 / np.pi
print(phi[0], phi.shape)

# Then, only take theta is less than 10 degrees.
intheta = theta < 10

print('From z axis = %s' % intheta.compress(intheta).shape)
E = E[intheta]
theta = theta[intheta]
phi = phi[intheta]

xyz = xyz[intheta]
vxyz = vxyz[intheta]

return E, theta, phi

def main():

eneall = []
thetall = []
phiall = []

nex = 10

for n in range(nex):  # Loop 10 times.
print('Repeat %d/%d times' % (n+1, nex))
E, theta, phi = main_prog()

for a,b,c in zip(E, theta, phi):
eneall.append(a)
thetall.append(b)
phiall.append(c)
print('Number of valid data = %d, %d, %d' % (len(eneall), len(thetall), len(phiall)))

eneall = np.array(eneall)
vall = np.sqrt(eneall)

fig = plt.figure()
ax.hist(eneall, weights=vall, range=(0, 60), bins=60)

theory = lambda E: 2 * 1e7/8. / 4 * (18 * np.pi) ** (-1.5) * E * np.exp(-E / 18) * 1 * 0.095454 * 0.04

ie = np.linspace(0, 60, 600)
je = theory(ie) * nex

ax.plot(ie, je)

ax.set_xlabel('Energy')
ax.set_ylabel('Weighted histogram: or j(E)')

fig.savefig('maxwell_dflux.png')

return eneall, thetall, phiall

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