Initial Radius
Grid Width
import numpy as np
import matplotlib.pyplot as plt
import asyncio
import random
import math
import js
import ctypes
from js import document
from pyscript import Element
# Initialize variables for animation
animating = False
collisions_off = False
walls_off = False
sp_plot = False
en_plot = False
async_flag = True
paused = True
wide, high = 500, 500
t, dt = 0, 1
tmax = 400
n = 100
r_init = 50
precision = 10
gridwidth = 100
sp_plots = 10
en_plots = 10
plot_int = 1
# Setup
def do_set():
global paused
if paused:
initialize_parameters()
initialize_pos_vel()
clear_dots()
clear_plot_1()
clear_plot_rest()
add_circles()
# Run the animation
def do_run():
global paused, animating, async_flag
if paused:
animating = True
paused = False
animate(t)
if async_flag:
async_flag = False
asyncio.ensure_future(main())
# Pause the animation
def do_pause():
global animating, paused
if animating:
animating = False
paused = True
# Reverse the velocity direction of each dot
def do_reverse():
global n, paused, animating, vel
if paused:
for i in range(n):
vel[i] *= -1
animating = True
paused = False
animate(t)
def initialize_parameters():
global n, t, tmax, time
global pos, vel, precision, r_init, r_space, imax
global step, coll, wall, animating, grid_side, gridwidth, spacing, dot_size
global collisions_off, walls_off, results, speed
global sp_plot, en_plot, sp_plots, en_plots, plot_int
global wide, high, large_grid, paused
global c_energy, d_energy, c_micro, d_micro, s_micro
global c_level, d_level, s_alloc
# Set container dimensions to match program constants
document.getElementById("cont").style.width=wide
document.getElementById("cont").style.height=high
# Get the input value of n
try:
n = int(document.getElementById("numberMolecules").value )
except Exception as e:
print("A number should be entered:", e)
# Get the maximum time steps
try:
tmax = int(document.getElementById("numberTimeSteps").value )
except Exception as e:
print("A number should be entered:", e)
# Get the energy precision (precision)
try:
precision = int(document.getElementById("Precision").value )
except Exception as e:
print("A number should be entered:", e)
if precision == 0:
precision = 1
# Get the initial radius (r_init)
try:
r_init = int(document.getElementById("Radius").value )
except Exception as e:
print("A number should be entered:", e)
if r_init > wide/2:
r_init = wide/2
# Calculate the space bands
imax = int(((wide/2) / r_init)**2 + 1)
r_space = np.zeros((imax, 2))
for i in range(imax):
r_space[i, 0] = math.sqrt(i)*r_init
# Set the grid size
grid_side = int(math.sqrt(n))
n = int(grid_side **2)
large_grid = False
try:
gridwidth = int(document.getElementById("gridwidth").value )
except Exception as e:
print("A number should be entered:", e)
if gridwidth <= 1:
gridwidth = 1
elif gridwidth >= wide:
gridwidth = wide
spacing = gridwidth / grid_side
if (document.getElementById("largeGrid").checked):
spacing = wide / grid_side
large_grid = True
else:
large_grid = False
show_message1('')
t = 0
step = 0
coll = 0
wall = 0
animating = True
paused = True
pos = np.ones((n, 2))
vel = np.zeros((n, 2))
time = np.zeros(tmax)
s_micro = np.zeros((tmax, imax))
s_alloc = np.zeros((tmax, imax-1))
c_energy = np.ones(n)
d_energy = np.ones(n)
c_micro = np.ones((tmax, n))
d_micro = np.zeros((tmax, n))
c_level = np.zeros(tmax)
d_level = np.zeros(tmax)
# Set the dot size based on the checkbox selection
dot_size = math.sqrt((0.1 * wide)**2 / n)
if plot_int < 1:
plot_int = 1
plt.close() # AVW
clear_plot_1() # AVW
clear_plot_rest() # AVW
# Set collisions OFF on the checkbox selection
if (document.getElementById("collisionsOff").checked):
collisions_off = True
else:
collisions_off = False
# Set walls OFF based on the checkbox selection
if (document.getElementById("wallsOff").checked):
walls_off = True
else:
walls_off = False
# Select the mode based on the radio button selection
en_plot = False
sp_plot = False
if (document.getElementById("gas").checked):
en_plot = False
sp_plot = False
elif(document.getElementById("energy").checked):
en_plot = True
else:
sp_plot = True
# Set microstate plot interval
if sp_plot:
plot_int = round(tmax / sp_plots)
elif en_plot:
plot_int = round(tmax / en_plots)
if plot_int < 1:
plot_int = 1
# Set results based on the checkbox selection
if (document.getElementById("results").checked):
results = True
else:
results = False
do_pause()
# Initialize position and velocity
def initialize_pos_vel():
global grid_side, pos, spacing, high, wide, vel
k = 0
for i in range(grid_side):
for j in range(grid_side):
pos[k,0] = float(i * spacing) + (spacing / 2) + (wide - grid_side * spacing)/ 2
pos[k,1] = float(j * spacing*(high/wide)) + (spacing*(high/wide) / 2) + (high - grid_side * spacing)/ 2
angle = random.uniform(0, 2 * math.pi)
vel[k] = (np.array((np.cos(angle), np.sin(angle))))
k = k + 1
# Add the dots
def add_circles():
global n
container = document.getElementById("cont") # a container for the circles
k = 0
for i in range(grid_side):
for j in range(grid_side):
circle=document.createElementNS("http://www.w3.org/2000/svg", "circle")
idstr = "cir"+str(k) # make the id a function k
circle.setAttribute( "id", idstr) # set the id
xloc = pos[k,0] # initial x location
circle.setAttribute( "cx", xloc)
yloc = pos[k,1] # initial y location
circle.setAttribute( "cy", yloc)
circle.setAttribute( "r", dot_size / 2) # set the radius
circle.setAttribute( "fill", "blue") # try a different color
circle.setAttribute( "value", 1) # create a value attribute
container.appendChild(circle) # add to container
k += 1 # k= 0, ... n-1
# Loop through time generating microstates
def animate(t):
global speed, tmax, paused, async_flag
if not animating:
return
# Set speed to the slider value
speed = int(document.getElementById("animation_speed").value)
update_microstates(t)
do_collisions_walls(t)
if sp_plot or en_plot:
plot_t_charts(t)
if t == (tmax - 1):
plot_charts(t)
show_messages()
paused = True
async_flag = True
# Update microstates
def update_microstates(t):
global n, tmax, time, speed, dt, wide, high
global pos, vel, precision, r_space, s_micro
global c_energy, d_energy, c_micro, d_micro
# Save the time
if t < tmax:
time[t] = t
else:
return
pos += speed * vel * dt # Move the dots
# Draw the dots
draw_dots()
# Calculate the spatial microstates
r_space[:, 1] = 0
for i in range(n):
radius = math.sqrt((pos[i, 0] - wide/2)**2 + (pos[i, 1] - high/2)**2)
jrange = len(r_space) - 1
for j in range(jrange):
if radius >= r_space[j, 0] and radius < r_space[j+1, 0]:
r_space[j, 1] += 1
s_micro[t] = r_space[:, 1]
# Calculate the energy microstates
for i in range(n):
c_energy[i] = precision*(vel[i, 0]**2 + vel[i, 1]**2) # Continuous energy
d_energy[i] = round(precision*(vel[i, 0]**2 + vel[i, 1]**2), 0) # Discrete energy
c_micro[t] = c_energy # Microstate is energy configuration
d_micro[t] = d_energy # Microstate is energy configuration
# Calculate collions with molecules and walls
def do_collisions_walls(t):
global n, dt, random_numbers, dot_size, step, coll, wall, pos, vel
global speed, collisions_off, walls_off, precision
global c_energy, d_energy, c_micro, d_micro, c_micro_t, d_micro_t
global c_level, d_level, s_micro, s_micro_t, s_alloc
step += 1
for i in range(n):
# Do Collisions
if collisions_off == False:
for j in range(n):
if i != j:
distance = ((pos[i, 0] - pos[j, 0])**2 + (pos[i, 1] - pos[j, 1])**2)**0.5
# Molecular collisions
if distance <= dot_size:
coll += 1
pos_i, vel_i = pos[i], vel[i]
pos_j, vel_j = pos[j], vel[j]
rel_pos, rel_vel = pos_i - pos_j, vel_i - vel_j
r_rel = rel_pos @ rel_pos
v_rel = rel_vel @ rel_pos
v_rel = 2 * rel_pos * v_rel / r_rel - rel_vel
v_cm = (vel_i + vel_j) / 2
vel_i = v_cm - v_rel/2
vel_j = v_cm + v_rel/2
vel[i] = vel_i
vel[j] = vel_j
pos[i] += speed * vel[i] * dt
pos[j] += speed * vel[j] * dt
# Wall collisions
if walls_off == False:
# Bounce off walls
if pos[i, 0] <= dot_size or pos[i, 0] >= 500 - dot_size:
wall += 1
vel[i, 0] *= -1
pos[i] += speed * vel[i] * dt
if pos[i, 1] <= dot_size or pos[i, 1] >= 500 - dot_size:
wall += 1
vel[i, 1] *= -1
pos[i] += speed * vel[i] * dt
# Discrete space
s_micro_t = s_micro[t]
s_micro_t = s_micro_t[:-1]
s_alloc[t] = s_micro_t / sum(s_micro_t)
# Discrete energy
d_micro_t = d_micro[t]
lev_count = np.zeros(n*precision)
for i in range(n*precision):
for item in d_micro_t:
if item == i:
lev_count[i] += 1
M = 0
for item in lev_count:
if item != 0:
M += 1
d_level[t] = M
# Continuous energy
c_micro_t = c_micro[t]
ones_count = 0
for item in c_micro_t:
if item == precision:
ones_count += 1
M = n - ones_count
if M < n:
M = M + 1
c_level[t] = M
# Draw the dots
def draw_dots():
global n, pos
for i in range(n):
idstr = "cir"+str(i)
circle = document.getElementById(idstr)
xloc = pos[i,0] # update x location
circle.setAttribute( "cx", xloc)
yloc = pos[i,1] # update y location
circle.setAttribute( "cy", yloc)
circle.setAttribute( "value", i) # update value
# Clear the dots
def clear_dots():
container = document.getElementById("cont")
while (container.hasChildNodes()):
container.removeChild(container.firstChild)
# Clear plot_1
def clear_plot_1():
container = document.getElementById("plot_1")
while (container.hasChildNodes()):
container.removeChild(container.firstChild)
# Clear rest of plots
def clear_plot_rest():
container = document.getElementById("plot_2")
while (container.hasChildNodes()):
container.removeChild(container.firstChild)
container = document.getElementById("plot_3")
while (container.hasChildNodes()):
container.removeChild(container.firstChild)
container = document.getElementById("plot_4")
while (container.hasChildNodes()):
container.removeChild(container.firstChild)
# Plot charts
def plot_charts(t):
global n, coll, speed, tmax, time
global pos, vel, precision, r_init, imax
global c_micro, d_micro, c_micro_t, d_micro_t, s_micro, s_micro_t
global c_level, d_level, s_alloc
global message1, message2
message1 = 'Molecules: ' + str(n) + ' Molecular Collisions: ' + str(coll) + ' Wall Collisions: ' + str(wall) + ' Speed: ' + str(speed)
message2 = 'Number of Time Steps: ' + str(tmax) + ' Energy Precision: ' + str(precision) + ' Initial Radius: ' + str(r_init)
plt.close() # AVW
clear_plot_1() # AVW
clear_plot_rest() # AVW
if en_plot or en_plot == sp_plot:
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(9, 3))
unique_rows, row_counts = np.unique(d_micro_t, axis=0, return_counts=True)
d_micro_t = np.column_stack((unique_rows, row_counts))
last_rec = int(d_micro_t[len(d_micro_t)-1, 0])
array = np.zeros((last_rec+1, 2))
for i in range(last_rec+1):
array[i, 0] = i
k = 0
for i in range(last_rec+1):
j = i - k
if array[i, 0] == d_micro_t[j, 0]:
array[i, 1] = d_micro_t[j, 1]
else:
k += 1
array[:, 0] = array[:, 0] / precision # Observed energy, comment out for Discrete actual
# Plot discrete energy distribution
axes[0].bar(array[:, 0], array[:, 1], width = .08, align = 'edge')
axes[0].set_title('Observed Occupation Numbers n='+str(n)+' c='+str(precision))
axes[0].set_xlabel('Observed Energy Level')
axes[0].set_ylabel('Occupation Number')
axes[0].set_xlim(0, max(array[:, 0]) + 0.8)
axes[0].set_ylim(0, max(d_micro_t[:, 1])+1)
unique_rows, row_counts = np.unique(c_micro_t, axis=0, return_counts=True)
c_micro_t = np.column_stack((unique_rows, row_counts))
c_micro_t[:, 0] =c_micro_t[:, 0] / precision # Observed energy, comment out for Discrete actual
axes[1].bar(c_micro_t[:, 0], c_micro_t[:, 1], width = .08, align = 'edge')
axes[1].set_title('Theoretical Occupation Numbers n='+str(n))
axes[1].set_xlabel('Exact Energy Level')
axes[1].set_ylabel('Occupation Number')
axes[1].set_xlim(0, max(array[:, 0]) + 0.8)
axes[1].set_ylim(0, max(d_micro_t[:, 1])+1)
# Show figures 1a and 1b
plt.subplots_adjust(left= .075, right= .965, top=.90, bottom=.15)
display(fig, target="plot_1")
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(9, 3))
# Plot discrete energy levels
en_time = time[:tmax-1]
d_level = d_level[:tmax-1]
axes[0].plot(en_time, d_level)
axes[0].set_title('Observed Energy Levels n='+str(n)+' c='+str(precision))
axes[0].set_xlabel('Number of Time Steps')
axes[0].set_ylabel('Number of Observable Energy Levels')
axes[0].set_ylim(0, max(c_level)+1)
# Plot continuous energy levels
c_level = c_level[:tmax-1]
axes[1].plot(en_time, c_level, 'green')
axes[1].set_title('Theoretical Energy Levels n='+str(n))
axes[1].set_xlabel('Number of Time Steps')
axes[1].set_ylabel('Number of Theoretical Energy Levels')
axes[1].set_ylim(0, max(c_level)+1)
if results:
print('Observed Energy Levels ', int(max(d_level)))
print('Theoretical Energy Levels', int(max(c_level)))
# Show figures 3a and 3b
plt.subplots_adjust(left= .075, right= .965, top=.90, bottom=.15)
display(fig, target="plot_2")
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(9, 3))
# Plot discrete microstate allocation
d_micro = d_micro[:tmax-1] / (n*precision)
axes[0].plot(en_time, d_micro)
axes[0].set_title('Observed Energy Microstates n='+str(n)+' c='+str(precision))
axes[0].set_xlabel('Number of Time Steps')
axes[0].set_ylabel('Molecular Energy Allocation')
# Plot continous microstate allocation
c_micro = c_micro[:tmax-1] / (n*precision)
axes[1].plot(en_time, c_micro)
axes[1].set_title('Theoretical Energy Microstates n='+str(n))
axes[1].set_xlabel('Number of Time Steps')
axes[1].set_ylabel('Molecular Energy Allocation')
if results:
print('Max End Energy Allocation', round(max(d_micro_t[:, 1]), 4), round(max(c_micro_t[:, 1]), 4))
# Show figures 4a and 4b
plt.subplots_adjust(left= .075, right= .965, top=.90, bottom=.15)
display(fig, target="plot_3")
if sp_plot or en_plot == sp_plot:
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(9, 3))
# Plot discrete spatial microstates
row = np.zeros(imax-1)
for i in range(imax-1):
row[i] = i + 1
axes[0].bar(row, s_micro_t, color='brown')
axes[0].set_title('Observed Occupation Numbers n='+str(n)+' r='+str(r_init))
axes[0].set_xlabel('Band Ranked by Radius')
axes[0].set_ylabel('Occupation Number')
axes[0].set_ylim(0, r_init*n/100)
# Plot discrete spatial allocation
axes[1].plot(time, s_alloc)
axes[1].set_title('Observed Spatial Microstates n='+str(n)+' r='+str(r_init))
axes[1].set_xlabel('Number of Time Steps')
axes[1].set_ylabel('Spatial Allocation of Molecules ')
if results:
print('Number of Spatial Bands ', int(imax-1))
print('Max Spatial Allocation ', round(100*max(s_alloc[:, imax-2]), 2), '% of Total')
# Show figures 8a and 8b
plt.subplots_adjust(left= .075, right= .965, top=.90, bottom=.15)
display(fig, target="plot_4")
if results:
print()
#Plot animated macrostate distributions
def plot_t_charts(t):
global n, precision, r_init, imax, sp_plot, en_plot, plot_int
global s_micro_t, c_micro_t, d_micro_t
document.getElementById("cont").style.width=0 # AVW
document.getElementById("cont").style.height=0 # AVW
if t > 0:
if sp_plot:
if t == 1 or t%plot_int == 0:
plt.close() # AVW
clear_plot_1() # AVW
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(9, 4))
# Plot spatial microstates
plt.clf()
row = np.zeros(imax-1)
for i in range(imax-1):
row[i] = i + 1
plt.bar(row, s_micro_t, color='brown')
plt.title('Observed Occupation Numbers n='+str(n)+' r='+str(r_init))
plt.xlabel('Band Ranked by Radius')
plt.ylabel('Occupation Number')
plt.ylim(0, r_init*n/100+1)
plt.subplots_adjust(left= .075, right= .965, top=.90, bottom=.15)
display(fig, target="plot_1")
elif en_plot:
if t == 1 or t%plot_int == 0:
plt.close() # AVW
clear_plot_1() # AVW
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(9, 4))
unique_rows, row_counts = np.unique(d_micro_t, axis=0, return_counts=True)
d_micro_t = np.column_stack((unique_rows, row_counts))
# Plot discrete energy distribution
plt.bar(d_micro_t[:, 0], d_micro_t[:, 1])
plt.title('Observed Occupation Numbers n='+str(n)+' c='+str(precision))
plt.xlabel('Observed Energy Level')
plt.ylabel('Occupation Number')
# Show figure
plt.subplots_adjust(left= .075, right= .965, top=.90, bottom=.15)
display(fig, target="plot_1")
def show_message1(m: str):
myCell = Element("message1")
myCell.write(m)
def show_messages():
global message1, message2
message1 = message1 + '\r\n' + message2
show_message1(message1)
# Function to start/continue the animation
async def main():
global animating, step, coll, tmax, async_flag
animating = True
for t in range(tmax):
animate(t)
document.getElementById("Step").innerText = str(step);
document.getElementById("Collision").innerText = str(coll);
await asyncio.sleep(.001)
# get all tasks
tasks = asyncio.all_tasks()
# cancel all tasks
for task in tasks:
task.cancel()
async_flag = True