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 = fig.add_subplot(111)
    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()