# apps121025_test_liouville.app05_multipleΒΆ

Multiple particle, trying to estimate the amount of work.

500 particle tracing takes ~7.568 sec for my iMac. 1 hr will calculate 240000 particles.

I am still not very sure how much particles needed.

''' Multiple particle, trying to estimate the amount of work.

500 particle tracing takes ~7.568 sec for my iMac.
1 hr will calculate 240000 particles.

I am still not very sure how much particles needed.

'''
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mp

from .app04_nrm_setup import kep_particle

def main():

### Reference condition at t=0.  For 2D.
pos0ref = np.array([1.2317, 0, 0])
vel0ref = np.array([-0.52318 * 0.5, 0.52318 * np.sqrt(3./4), 0])

### Reference condition at t=3e5.
pos1ref, vel1ref = kep_particle(pos0ref, vel0ref)
pos1ref = pos1ref[-1, :]
vel1ref = vel1ref[-1, :]

posdel = np.array([0.001, 0.001, 0])    # Dr at t=0
veldel = np.array([0.001, 0.001, 0])      # Dv at t=0

drdv = posdel[0] * posdel[1] * veldel[0] * veldel[1]
print(drdv)

npart = 500

pos0list = np.zeros([3, npart])
pos1list = np.zeros([3, npart])

vel0list = np.zeros([3, npart])
vel1list = np.zeros([3, npart])

for n in range(npart):  ###  Test particles
### Calculate the position and velocity
### from the parameters defined above.
### Position (and velocity) should have
### the uniformly random from pos0ref \pm posdel/2.

posx = np.random.uniform(low=-posdel[0] / 2., high=posdel[0] / 2.) + pos0ref[0]
posy = np.random.uniform(low=-posdel[1] / 2., high=posdel[1] / 2.) + pos0ref[1]
velx = np.random.uniform(low=-veldel[0] / 2., high=veldel[0] / 2.) + vel0ref[0]
vely = np.random.uniform(low=-veldel[1] / 2., high=veldel[1] / 2.) + vel0ref[1]

pos = np.array([posx, posy, 0])
vel = np.array([velx, vely, 0])

pos0list[:, n] = pos
vel0list[:, n] = vel

pos1, vel1 = kep_particle(pos, vel)
pos1list[:, n] = pos1[-1, :]
vel1list[:, n] = vel1[-1, :]

plt.figure()
plt.subplot(231)
plt.scatter(pos0list[0, :] - pos0ref[0], pos0list[1, :] - pos0ref[1])
plt.subplot(232)
plt.scatter(vel0list[0, :] - vel0ref[0], vel0list[1, :] - vel0ref[1])
plt.subplot(233)
plt.scatter(pos0list[0, :] - pos0ref[0], vel0list[0, :] - vel0ref[0])
plt.subplot(234)
plt.scatter(pos1list[0, :] - pos1ref[0], pos1list[1, :] - pos1ref[1])
drbox = mp.Rectangle([0, 0], posdel[0], posdel[1], color='r')
plt.subplot(235)
plt.scatter(vel1list[0, :] - vel1ref[0], vel1list[1, :] - vel1ref[1])
dvbox = mp.Rectangle([0, 0], veldel[0], veldel[1], color='r')
plt.subplot(236)
plt.scatter(pos1list[0, :] - pos1ref[0], vel1list[0, :] - vel1ref[0])
dvbox = mp.Rectangle([0, 0], posdel[0], veldel[0], color='r')