# apps121025_test_liouville.app06_volume_calcΒΆ

Using 4 particle, calculate the area.

''' Using 4 particle, calculate the area.
'''
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mp

from .app04_nrm_setup import kep_particle

def main_cal(withplot=False, figure=None):

### 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.0001, 0.0001, 0])    # Dr at t=0
veldel = np.array([0.0001, 0.0001, 0])      # Dv at t=0

npart = 4

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, :]

if withplot:
if figure == None:
figure = plt.figure()
ax.scatter(pos0list[0, :] - pos0ref[0], pos0list[1, :] - pos0ref[1], marker='.')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.scatter(vel0list[0, :] - vel0ref[0], vel0list[1, :] - vel0ref[1], marker='.')
ax.set_xlabel('Vx')
ax.set_ylabel('Vy')
ax.scatter(pos0list[0, :] - pos0ref[0], vel0list[0, :] - vel0ref[0], marker='.')
ax.set_xlabel('X')
ax.set_ylabel('Vx')
ax.scatter(pos0list[1, :] - pos0ref[1], vel0list[1, :] - vel0ref[1], marker='.')
ax.set_xlabel('Y')
ax.set_ylabel('Vy')

ax.scatter(pos1list[0, :] - pos1ref[0], pos1list[1, :] - pos1ref[1], marker='.')
ax.set_xlabel('X')
ax.set_ylabel('Y')
drbox = mp.Rectangle([-posdel[0]/2., -posdel[1]/2.], posdel[0], posdel[1], color='none', ec='red')

ax.scatter(vel1list[0, :] - vel1ref[0], vel1list[1, :] - vel1ref[1], marker='.')
ax.set_xlabel('Vx')
ax.set_ylabel('Vy')
drbox = mp.Rectangle([-veldel[0]/2., -veldel[1]/2.], veldel[0], veldel[1], color='none', ec='red')

ax.scatter(pos1list[0, :] - pos1ref[0], vel1list[0, :] - vel1ref[0], marker='.')
ax.set_xlabel('X')
ax.set_ylabel('Vx')
drbox = mp.Rectangle([-posdel[0]/2., -veldel[0]/2.], posdel[0], veldel[0], color='none', ec='red')

ax.scatter(pos1list[1, :] - pos1ref[1], vel1list[1, :] - vel1ref[1], marker='.')
ax.set_xlabel('Y')
ax.set_ylabel('Vy')
drbox = mp.Rectangle([-posdel[1]/2., -veldel[1]/2.], posdel[1], veldel[1], color='none', ec='red')

# Volume calculation
from irfpy.util.graminan import Graminan
### 4-D space (x, y, vx, vy) for each particle.

# At t=0
px0 = pos0list[0, :] - pos0ref[0]
py0 = pos0list[1, :] - pos0ref[1]
vx0 = vel0list[0, :] - vel0ref[0]
vy0 = vel0list[1, :] - vel0ref[1]

G0 = Graminan([px0, py0, vx0, vy0])   # Indeed, transpose is needed, but in this case, no problem as G.Transpose = G.
print("Phase space volume (T=0)", G0.get_volume())

# At t=3e4
px1 = pos1list[0, :] - pos1ref[0]
py1 = pos1list[1, :] - pos1ref[1]
vx1 = vel1list[0, :] - vel1ref[0]
vy1 = vel1list[1, :] - vel1ref[1]
G1 = Graminan([px1, py1, vx1, vy1])   # Indeed, transpose is needed, but in this case, no problem as G.Transpose = G.
print("Phase space volume (T=3e4)", G1.get_volume())

return (G0.get_volume(), G1.get_volume())

def main():
nexam = 100
vol0 = np.zeros([nexam])
vol1 = np.zeros([nexam])

fig = plt.figure(figsize=(16, 8))
for iexam in range(nexam):
v0, v1 = main_cal(withplot=True, figure=fig)
vol0[iexam] = v0
vol1[iexam] = v1
print(v0, v1)
fig.savefig('app06_volume_calc0.png')

plt.figure()
plt.scatter(np.log10(vol0), np.log10(vol1))
plt.plot([-20, -15], [-20, -15])
plt.xlabel('Log phasevolume t=0')
plt.ylabel('Log phasevolume t=1')
plt.savefig('app06_volume_calc1.png')

if __name__ == "__main__":
import time
t0 = time.time()
main()
print('Time ellapsed', time.time() - t0)