import sys
import pygame
import random
import numpy as np
import os
program_dir = os.getcwd() # Get the Folder path this Python Programm is running from.
image_dir = program_dir + '\Images' # Directory to PNG Image files
#if not os.path.exists(image_dir): # If no directory folder for PNG files yet
# os.mkdir('Images') # then make new directory.
signalcollapsepoint = 0 # Program start point of Wave, oint where wave collapsed on other side.
signalactivity = 0 # Total number of connections that are spiking/on/set/1 during this time cycle.
defaultpointcount = 1024 # Number of particles when no StartupXYZ.txt file. Symmetrical are 32, 360.
pressureRampStep = .0025
pressureMin = 1
pressureMax = 1
pressureRampDirection = 1
adjustR = 1.01 # Shaking pressure like in jar of water, can help jog into place.
forceCutOff = 2 # Amount of spacing distance at which forces are too small to influence.
groupCutOff = 3.5
connectCutOff = 1.3 # Amount of spacing between neighbors to connect to.
rotationstep = .0125 # Amount of rotation per timestep.
rotation = 0 # Initialize rotation angle, to draw all points on screen.
epsilon = 1 # Lennard Jones variable, remains 1, normalized.
sigma = 1 # Lennard Jones variable, remains 1, normalized.
dt = .025 # Timestep size, how (slower but) precise it moves.
maxf = 2 # Velocity limit for when enough force to fly off sphere.
popwaitcount = 400
popwait = 0
pops = 0
height = 520 # Size in pixels for each of 2 squares to draw inside.
screensize = (height*2, height) # Pygame variable used to initialize screen window.
screenradius = int(height*.43) # Radius of vesicle in pixels, when drawn to screen.
xc1 = int(screensize[0]*.25) # Center x screen location for left square.
xc2 = int(screensize[0]*.75) # Center x screen location for square to right.
yc = int(screensize[1]*.5) # Center y screen location point for both squares.
timestepnumber = 0 # Counter for number of times program was stepped/looped one moment in time.
pi = 3.14159265358979323846264338327950288412
radian = pi * 2
threeQtrRadian = radian * 3/4
Black = (0, 0, 0)
Brown = (55, 0, 0)
Red = (255, 0, 0)
Orange = (255, 165, 0)
Yellow = (220, 220, 0)
Green = (0, 200, 0)
Blue = (0, 0, 150)
Violet = (255, 0, 255)
Gray = (120, 120, 120)
White = (255, 255, 255)
Cyan = (0, 255, 255)
# Colors used in electronics to mark resistors and other components with bands, to make easy to read not tiny numbers.
# These are used to similarly identify number of connections each cell has by color, to make easy to see what they are.
# 0=Black, 1=Brown, 2=Red, 3=Orange, 4=Yellow, 5=Green, 6=Blue, 7=Violet, 8=Gray, 9=White
ElectronicColorCode = Black, Brown, Red, Orange, Yellow, Green, Blue, Violet, Gray, White
# Colors used to draw connections that are either 0, 1 = action potential.
InputSignalColor = Black, Cyan # Colors for connections that just propagated the bit to next point/cell.
OutputSignalColor = Black, Red # Full on for entire timecycle is Red.
def csvwrite(sourcelist, filename): # Write Comma Separated Value list of lists [,,] or tuples (,,) to disk.
# For control over file behavior csv or other import module is not used for reading or writing.
with open(filename, 'w') as f: # Open file to write over, using 'f' as its name:
for line in sourcelist: # For each line of values in source list:
fs = str(line) # Convert line to a string containing commas.
f.write(fs[1:-1] + "\n") # Write all but [] or () chars, line feed at end.
# File to read must have no blank lines or comments added to it, comma separated values only.
def xyzread(filename): # Read disk file back into list of lists.
x = []
y = []
z = []
try:
with open(filename, 'r') as f: # Open file in read mode using 'f' as its name:
for line in f: # For each comma separated line of code in 'f':
fs = line.rstrip('\n') # Remove line x,y,z readings in line into s
if len(fs) > 0:
fs = fs.split(',') # Split line at commas into list named s
x.append(float(fs[0])) # Save first value in s as x reading.
y.append(float(fs[1])) # Save second value in s as y reading.
z.append(float(fs[2])) # Save third value in s as z reading.
except IOError: # XYZ File does not exist. Start new vesicle.
x = [0.0 for _ in range(defaultpointcount)]
y = [0.0 for _ in range(defaultpointcount)]
z = [0.0 for _ in range(defaultpointcount)]
# Randomly distribute points inside sphere by filling a cube but only save points inside.
for i in range(defaultpointcount): # For each x,y,z reading to save in list.
ra = screenradius + 1
while ra > screenradius: # Loop until radius is inside sphere.
x[i] = random.uniform(-screenradius, screenradius)
y[i] = random.uniform(-screenradius, screenradius)
z[i] = random.uniform(-screenradius, screenradius)
az, el, ra = cart2sphere(x[i], y[i], z[i]) # Generate random Cartesian point inside a cube.
return x, y, z
# Write Comma Separated Value list of lists [,,] or tuples (,,) to disk.
def xyzwrite(filename):
# For control over file behavior csv or other import module is not used for reading or writing.
with open(filename, 'w') as f: # Open file to write over, using 'f' as its name:
for i in range(points):
sxyz = str([xp[i], yp[i], zp[i]]) # For each value in the 3 source lists
f.write(sxyz[1:-1] + "\n") # Write all but [] or () chars, line feed at end.
screen.fill(White) # Make screen background flash white to indicate file save.
# File to read must have no blank lines or comments added to it, comma separated values only.
def csvread(filename): # Read disk file back into list of lists.
destination = [] # Clear destination list lists are appended to.
with open(filename, 'r') as f: # Open file in read mode using 'f' as its name:
for line in f: # For each comma separated line of code in 'f':
fs = line.rstrip('\n')
d = [] # Empty list named 'd' to fill with elements.
if len(fs) > 0:
fs = fs.split(',') # Split line at commas into list named 's'
for e in fs: # For each string element listed in 's' list:
d.append(float(e)) # Append floating point element value to list.
destination.append(d) # Append list of elements to destination list.
return destination # Returns list of lists only = [[,,],[,,],[,,]]
# Calculate spherical angle force vectors of a point are together most pointing towards.
# Magnitude becomes 0 when neighboring points are in opposite or opposing each other.
# icenter is the point to calculate.
def vectorsum(icenter, listname):
xsum = 0.0
ysum = 0.0
zsum = 0.0
for i in listname:
xctr = xp[icenter]
yctr = yp[icenter]
zctr = zp[icenter]
xsum += xp[i] - xctr
ysum += yp[i] - yctr
zsum += zp[i] - zctr
az, el, ra = cart2sphere(xsum, ysum, zsum)
return az, el, ra
# Cartesian to Spherical coordinates.
def cart2sphere(x, y, z):
hxy = np.hypot(x, y)
return np.arctan2(y, x), np.arctan2(z, hxy), np.hypot(hxy, z)
# Spherical to Cartesian coordinates.
def sphere2cart(az, el, ra):
theta = ra * np.cos(el)
x = theta * np.cos(az)
y = theta * np.sin(az)
z = ra * np.sin(el)
return x, y, z
# To make program faster: Return a list of nearbby points less than (too far to attract) CutOff distance away.
def update_grouplist():
print(timestepnumber, 'Compiling Group list')
groupsum = 0
for i in range(points): # For all the points in array
grouplist[i] = []
for j in range(points): # For all the other points in array
if i != j: # Except itself
dx = xp[j] - xp[i] # Calculate distance between points along x.
if abs(dx) <= groupR: # If (neg made positive) x distance OK continue.
dy = yp[j] - yp[i] # Calculate distance between points along y.
if abs(dy) <= groupR: # If (neg made positive) y distance OK continue.
dz = zp[j] - zp[i] # Calculate distance between points along z.
if abs(dz) <= groupR: # If (neg=pos.) z OK then time consuming sqrt precision.
dist = np.sqrt((dx * dx) + (dy * dy) + (dz * dz)) # Calculate 3D distance.
if dist < groupR: # If distance is within group range then.
grouplist[i].append(j) # Append i and j points to the end of list.
groupsum += len(grouplist[i])
print("Updated Group list, ", len(grouplist), "total distances, average", groupsum / points, "per point.")
# Return a list of points within cutoff range.
def update_cutofflist():
print(timestepnumber, 'Compiling CutOff list')
cutoffsum = 0
for i in range(len(grouplist)):
cutofflist[i] = []
for j in grouplist[i]:
if i != j: # Except itself
dx = xp[j] - xp[i] # Calculate distance between points along x.
if abs(dx) <= cutoffR: # If (neg made positive) x distance OK continue.
dy = yp[j] - yp[i] # Calculate distance between points along y.
if abs(dy) <= cutoffR: # If (neg made positive) y distance OK continue.
dz = zp[j] - zp[i] # Calculate distance between points along z.
if abs(dz) <= cutoffR: # If (neg=pos.) z OK then time consuming sqrt precision.
dist = np.sqrt((dx * dx) + (dy * dy) + (dz * dz)) # Calculate distance between.
if dist < cutoffR: # If distance is within cutoff range then
cutofflist[i].append(j) # Append i and j points to the end of list.
cutoffsum += len(cutofflist[i])
print("Updated CutOff list, ", len(cutofflist), "total distances, average", cutoffsum/points, "per point.")
return cutoffsum/points
# Return a list of particles within group range, and another list for close enough to form in/out connections to.
def update_connections():
global connections
for i in range(points):
connectlist[i] = []
for j in cutofflist[i]:
if i != j:
dx = xp[j] - xp[i] # Calculate x,y,z distances.
dy = yp[j] - yp[i]
dz = zp[j] - zp[i]
dist = np.sqrt((dx * dx) + (dy * dy) + (dz * dz)) # Calculate 3D distance, a time consuming command.
if dist < connectR: # If close enough to connect together then
connectlist[i].append(j) # save j point at this i, makes list of what connects to where.
connections[i] = len(connectlist[i])
# Update (azimuth) longitude and (elevation) latitude coordinates. Radius for all points is that of sphere, not needed.
def update_angles():
global latlong
latlong = []
for i in range(points):
az, el, ra = cart2sphere(xp[i], yp[i], zp[i]) # Spherical azimuth,elevation = longitude,latitude.
latlong.append([-el, az]) # Elevation is negative, north pole is positive latitude.
# print(latlong)
def update_positions():
global cutofflist
if holdradius:
for i in range(points):
# Hold to sphere. Convert points to spherical coordinates then back using sphere radius as center distance.
xp[i], yp[i], zp[i] = sphere2cart(np.arctan2(yp[i], xp[i]), # Calc azimuth.
np.arctan2(zp[i], np.hypot(xp[i], yp[i])), # Calc elevation.
vesicleR) # Desired sphere radius.
# Initialize velocity vector arrays for summing forces, combined angles and magnitudes.
xi = [0 for _ in range(points)]
yi = [0 for _ in range(points)]
zi = [0 for _ in range(points)]
xj = [0 for _ in range(points)]
yj = [0 for _ in range(points)]
zj = [0 for _ in range(points)]
# Clear sum, minimum and maximum
ljt = 0 # Initialize variable for totaling Lennard-Jones force.
ljmaximum = -10000 # Initialize for finding max, will be made greatest reading.
ljmaximum_i = 0 # For also saving max 'i' number of greatest reading.
ljminimum = 100000000000 # Initialize for finding min, will be made greatest reading.
ljminimum_i = 0 # For also saving min 'i' number of greatest reading.
ljminimum_j = 0 # For also saving min 'j' number of greatest reading.
atmax = 0
for i in range(len(cutofflist)): # For all the points in array.
ljsum[i] = 0
ljt1 = 0
for j in cutofflist[i]:
# Calculate x,y,z distances.
dx = xp[j] - xp[i]
dy = yp[j] - yp[i]
dz = zp[j] - zp[i]
# Use the three distances for calculating an always positive 3D distance.
dist = np.sqrt((dx * dx) + (dy * dy) + (dz * dz))
# Calculate Lennard Jones force.
ljr = dist/spacing # Lennard-Jones radius. At r=1 force is 0, spacing=r*2
lj = 4 * epsilon * ((sigma/ljr)**12 - (sigma/ljr)**6)
# Reduce timestep motion for smooth movement.
f = lj * dt
if f > maxf: # Limit amount of displacement to not fly off screen.
f = maxf
atmax = 1
# f = np.sqrt((x * x) + (y * y) + (z * z))
ljt1 += lj # Add to sum for average of all points.
ljsum[i] += f # Used for brightness to draw on screen.
if lj > ljmaximum: # If force is greater than last greatest then:
ljmaximum = lj # Save maximum reading.
ljmaximum_i = i # Also save max 'i' number of greatest reading.
if lj < ljminimum: # If force is lower than last lowest then:
ljminimum = lj # Save minimum reading.
ljminimum_i = i # Also save min 'i' number of lowest reading.
ljminimum_j = j # Also save min 'j' number of lowest reading.
# if f < -maxf: # Limit amount of displacement to not fly off screen.
# f = -maxf
# Use distance for calculating spherical azimuth and elevation angle r, as in cart2sphere.
x, y, z = sphere2cart(np.arctan2(dy, dx), np.arctan2(dz, dist), f)
# Update x,y,z force acting on i point, from this j point.
xi[i] -= x
yi[i] -= y
zi[i] -= z
# Update other point, j. For every action there is equal and opposite reaction.
xj[j] += x
yj[j] += y
zj[j] += z
# ljt += ljt1 / len(cutofflist[i]) # Add to sum for average of all points.
listLen = len(cutofflist[i])
if listLen > 0: # List length must be greater than zero or division by zero error.
ljt += ljt1 / len(cutofflist[i]) # Add to sum for average of all points.
# Add force displacement amounts, to points.
for i in range(points):
xp[i] += xi[i] + xj[i]
yp[i] += yi[i] + yj[i]
zp[i] += zi[i] + zj[i]
listLen = len(cutofflist[i])
if listLen > 0: # List length must be greater than zero or division by zero error.
ljsum[i] = ljsum[i] / len(cutofflist[i])
return ljt/points, ljminimum, ljminimum_i, ljminimum_j, ljmaximum, ljmaximum_i, atmax
def render():
# Precalculate all the rotated X,Y screen locations and their colors.
for i in range(points):
az, el, r = cart2sphere(xp[i], yp[i], zp[i]) # Cartesian to Spherical (az, el, r)
xr, yr, zr = sphere2cart(az + rotation, el, r) # Spherical to Cartesian.
xyzi[i] = [xr, yr, zr, i]
unsortedxyzi[i] = xyzi[i]
# To properly draw on top of each other sort the list according to what is furthest away, [1] = by Y axis.
xyzi.sort(key=lambda x: x[1])
if drawcolors != 0:
for i in range(points): # Draw all points to screen.
drawvesicle(xc1, 1, i) # Draw left view i point to screen in proper color.
drawvesicle(xc2, -1, points-1-i) # Draw right view i point to screen in proper color.
if drawwaves is True:
for i in range(points):
drawsignallines(xc1, 1, i) # Draw left view signal lines from given i point.
drawsignallines(xc2, -1, points-1-i) # Draw right view signal lines from given i point.
if axeslines: # Draw x,y,z axis lines to screen.
ra = vesiclestartR * 1.1 # Length of axis line to draw.
drawaxisline(0, 0, ra, Yellow)
if rotation > threeQtrRadian: # Green x behind blue y?
drawaxisline(ra, 0, 0, Green)
drawaxisline(0, ra, 0, Blue)
if rotation < threeQtrRadian: # Green x on top of blue y?
drawaxisline(ra, 0, 0, Green)
if pops > 0:
drawpop(xc1, 1, White, xpop, ypop, zpop)
drawpop(xc2, -1, White, xpop, ypop, zpop)
drawpop(xc1, 1, Cyan, xnew, ynew, znew)
drawpop(xc2, -1, Cyan, xnew, ynew, znew)
pygame.display.update() # Display everything drawn into screen buffer.
# Draw one i point filled circle in sorted xyz list. For next loop in render draws all vesicle points.
def drawpop(xc, direction, c, x, y, z):
az, el, r = cart2sphere(x, y, z) # Cartesian to Spherical (az, el, r)
xr, yr, zr = sphere2cart(az + rotation, el, vesicleR) # Spherical to Cartesian.
if yr*direction < 0:
return
xs = int(xr * direction + xc) # Screen pixel x location to draw at.
ys = int(zr + yc) # Screen pixel y location to draw at.
cr = int(spacing/2) + 12 # Circle radius after distance away is factored in.
pygame.draw.circle(screen, c, [xs, ys], cr, 5) # Draw circle to screen.
# Draw one i point filled circle in sorted xyz list. For next loop in render draws all vesicle points.
def drawvesicle(xc, direction, i):
x, y, z, ui = xyzi[i] # Get x,y,z and unsorted 'i' from sorted list.
if y*direction < -spacing:
return
c = Black # Default fill color is black, when drawcolors=0
if drawcolors == 1: # Keyboard 'c' changes, 1=LennardJones colors.
c = ljcolor(ui) # Use Lennard-Jones amount, warm to hot colors.
if drawcolors == 2: # Keyboard 'c' changes, 2=Connections based colors.
c = connectcolor(ui) # Electronic color code, number of connections each.
xs = int((x*direction) + xc) # Screen pixel x location to draw at.
ys = int(z + yc) # Screen pixel y location to draw at.
d = (((y*direction)+vesiclestartR)/spherewidth*.4) + .6 # Fraction according to how far away, depth.
if d > 1: # If fractional depth is greater than 1 then
d = 1 # reduce to 1, close as possible to observer.
cr = int(d * spacing * .5) + 1 # Circle radius after distance away is factored in.
pygame.draw.circle(screen, c, [xs, ys], cr, 0) # Draw filled circle to screen.
if outlines:
pygame.draw.circle(screen, White, [xs, ys], cr, 1) # Draw outline circle around it.
if drawnumbers:
fnt = pygame.font.SysFont('arial', int(halfspacing))
txtsurf = fnt.render(str(ui), False, White)
textwidth = fnt.size(str(ui))[0]
screen.blit(txtsurf, (xs - textwidth/2, ys - quarterspacing))
# Draw signal lines outward from one given i in sorted xyz list. For next loop in render draws all vesicle points.
def drawsignallines(xc, direction, i):
global connectlist
global outs
x, y, z, ui = xyzi[i] # Get variables from sorted list of x,y,z points.
if y*direction < -spacing:
return
x1 = int((x*direction) + xc) # Screen pixel x location to draw at.
y1 = int(z + yc) # Screen pixel y location to draw at.
d = (((y*direction)+vesicleR) / spherewidth) # Fraction according to how far away, depth.
if d > 1: # If fractional depth is greater than 1 then
d = 1 # reduce to 1, close as possible to observer.
outbits = outs[ui]
for j in range(len(connectlist[ui])):
x2, y2, z2, ui2 = unsortedxyzi[connectlist[ui][j]]
b = outbits[j]
if b == 1:
thick = int((b * d * spacing / 3) + 1)
c = OutputSignalColor[b]
if c != Black:
pygame.draw.line(screen, c, [x1, y1], [(x2*direction)+xc, z2+yc], thick)
inbits = ins[ui]
for j in range(len(connectlist[ui])):
x2, y2, z2, ui2 = unsortedxyzi[connectlist[ui][j]]
b = inbits[j]
thick = int((b * d * spacing / 9) + 1)
c = InputSignalColor[b]
if c != Black:
pygame.draw.line(screen, c, [x1, y1], [(x2*direction)+xc, z2+yc], thick)
def save():
os.system("ffmpeg -r 1 -i img%01d.png -vcodec mpeg4 -y movie.mp4")
def drawaxisline(x, y, z, c):
az, el, ra = cart2sphere(x, y, z) # Cartesian to Spherical (az, el, r)
x1, y1, z1 = sphere2cart(az + rotation, el, ra) # Convert from Spherical to Cartesian.
pygame.draw.line(screen, c, (xc1, yc), (x1+xc1, z1+yc), 1)
# Vary brightness according to amount of force. Positive is orange to white, negative (attraction) blue.
def ljcolor(t):
if ljsum[t] >= 0:
poscolor = int(ljsum[t] * (256/maxf) * 10)
if poscolor > 255:
poscolor = int(poscolor/200)
if poscolor > 255:
poscolor = 255
return 255, 240, poscolor
else:
return poscolor, poscolor * .75, 0
if ljsum[t] < 0:
negcolor = int(-ljsum[t] * 2000)
if negcolor > 255:
negcolor = 255
return 0, negcolor*.6, negcolor
def connectcolor(n):
cns = len(connectlist[n])
if cns > 9:
cns = 9
return ElectronicColorCode[cns]
def timestepsignals():
global ins
global outs
global signalactivity
global signalcollapsepoint
signalactivity = 0
ins = [[] for _ in range(points)] # Clear Input list output bits from neighbors are gathered into.
# For all points in array, gather bits from neighbors to a list for each. Example: ins[i,[1,0,0,0,1,0]]
for i in range(points): # For each point in vesicle, each with connections to 'j's
for j in range(len(connectlist[i])): # For each of the (usually 6) 'j' connections to other cells,
if j < len(outs[i]): # If output connection still exists in connection list then
b = outs[i][j] # bit to save is from output list.
else: # Else
b = 0 # to prevent program error append with 0 to fill void later.
ins[connectlist[i][j]].append(b) # Append bit from 'j' output to corresponding input line at 'i'
# If any bits are 1 then propagate wave using NOT to output inverted/opposite bit pattern. 1=0 and 0=1
for i in range(points): # For each 'i' point in vesicle, each with connections to 'j's
outs[i] = []
setbits = sum(ins[i])
signalactivity += setbits # Signal activity equals total number of bits set, entire surface.
if setbits > 0: # If there are any bits set then
outs[i] = [int(not b) for b in ins[i]] # all 1's in output list become 0, and all 0 becomes 1, negates.
else: # Else
outs[i] = ins[i] # Output bits remain all 0.
# If no signal activity then start wave by setting all output connections of a randomly selected place to 1.
if signalactivity == 0: # If no signal activity at all then
# startwave(random.randint(0, points-1)) # Randomly choose the point to start wave from.
startwave(4) # Start wave at same given point each time.
# If uncontrolled signal activity then clear all outputs.
if signalactivity > points * 2: # If more than average two connections per cell then
clearwave()
# Start wave from cell at point number given by value of i
def startwave(i):
global signalactivity
outs[i] = [1] * len(connectlist[i]) # Output 1, action potential, to all neighbors.
signalactivity += 6
print('~~~~~~~~~~ Wave started from point number', i, '~~~~~~~~~~')
# For when extreme hexagonal misalignment causes uncontrollable seixure-like activity.
def clearwave():
for i in range(points): # for all points in network.
outs[i] = [0] * len(outs[i]) # fill list with 0's, same length as there are out connections.
screen.fill(White)
myfont = pygame.font.SysFont('arial', 68)
textsurface = myfont.render("STOP", False, Black)
screen.blit(textsurface, (height * .82, height * .8))
print('~~~~~~~~~~ All connection outputs were set to 0 ~~~~~~~~~~')
# To Squeeze and Release pressure from all around.
def pressureRamp(press, direction):
newpressure = press + (pressureRampStep * direction)
if newpressure > pressureMax:
newpressure = pressureMax
direction = -1
if newpressure < pressureMin:
newpressure = pressureMin
direction = 1
return newpressure, direction
def listlennardjones():
lj = 0
for i in range(32):
d = i / 16
# ljr = d / sp # Lennard-Jones radius. At ljr=1 force is 0.
if d > 0: # Prevent divide by zero error when LJ radius is zero.
lj = 4 * epsilon * ((sigma / d) ** 12 - (sigma / d) ** 6)
print(lj, d)
#######################################################################################################################
xp, yp, zp = xyzread("StartupXYZ.txt") # Load either startup file, or random positions for new vesicle.
points = len(xp) # Number of points total in vesicle.
ljavr = 0.0
ljmax = 0.0
ljsum = [0 for _ in range(points)] # Initialize Lennard-Jones sum array.
xyzi = [[] for _ in range(points)] # List to sort, to draw according to their distance from observer.
unsortedxyzi = [[] for _ in range(points)] # Above list before sorting.
grouplist = [[] for _ in range(points)] # Nearby points in larger group, stay with for very long time.
cutofflist = [[] for _ in range(points)] # Points within cutoff range.
connectlist = [[] for _ in range(points)] # Connections from each i point, to neighbors j.
connections = [[] for _ in range(points)] # Connections from each i point, to neighbors j.
latlong = [[] for _ in range(points)] # Spherical coordinates for each point.
cutoffconnects = [[] for _ in range(points)] # Number of connections, usually 6.
ins = [[] for _ in range(points)] # Input (action potential) bits, from each 1 bit connection.
outs = [[] for _ in range(points)] # Output (action potentials) bits, to each 1 bit connection.
# PyGame Initialization
pygame.init()
screen = pygame.display.set_mode(screensize)
pygame.display.set_caption("Lennard Jones Potential Vescicle. A=Axis, C=ColorCode, F=Forces, H=HoldToSphere, N=indexNumber, P=Pressure, W=Waves, Spacebar=Pause, S=Save-XYZ.txt")
fps = pygame.time.Clock()
pygame.font.init()
# Display 'Initializing nnnnn points' message.
myfont = pygame.font.SysFont('arial', 50)
s = 'Initializing ' + str(points) + ' points'
textsurface = myfont.render(s, False, White)
screen.blit(textsurface, (10, 10))
pygame.display.update() # Display everything drawn into screen buffer.
fps.tick(2)
# Initialize vesicle size parameters.
sas = 4*pi*(screenradius*screenradius) # Area of sphere at this given screen radius.
sar = sas/points # Divide for area of one Lennard-Jones radius, between points.
rbp = np.sqrt(sar/pi) * 2 # Radius between points, given area.
spacing = rbp-((1/points)*rbp) # Radius between points, ideal spacing.
halfspacing = spacing / 2
quarterspacing = spacing / 4
listlennardjones()
avcutoffs = 0
cutoffcalctime = 10 # Time steps between updates of i,j points within cut-off range.
cutoffupdate = cutoffcalctime # Time cycles per update of cutoff list. +1 to upda
groupcalctime = 10 # Time steps between updates of i,j points within larger group.
groupupdate = groupcalctime # Time cycles per update of group list. +1 to update on startup.
xpop = 0
ypop = 0
zpop = 0
xnew = 0
ynew = 0
znew = 0
# Initialize keyboard entered toggle controls.
paused = False
axeslines = False
pointschange = True
forces = True # Must start True or list of points within cutoff range is not generated.
drawcolors = 2 # 0=none, 1=LennardJones, 2=Connections
pressureChanges = False
drawwaves = True # To enable or disable drawing of signal propagation lines.
drawnumbers = False
holdradius = True
if points < 2000:
outlines = True
else:
outlines = False
screen.fill(Black) # Clear screen of "Intitializing nnnn Points" message.
while True:
for event in pygame.event.get():
if event.type == pygame.QUIT:
pygame.quit()
sys.exit()
# Space=pause A = Axes,x=green,y=bluez=yellow. F=force. H=hold Radius. O=outlines.
# P=Pressure. S=save. W=waves. C= color, LJ, # connections. none waves only.
if event.type == pygame.KEYUP:
if event.key == pygame.K_SPACE: # Space bar = pause.
paused = not paused
if event.key == pygame.K_a: # a = Axis lines in left draw: x=green, y=blue, z=yellow.
axeslines = not axeslines
if event.key == pygame.K_f: # f = Lennard-Jones Force on/off. None = fast stable waves.
forces = not forces
if event.key == pygame.K_h: # h = Hold to (an inner volume) sphere. None = particles drift
holdradius = not holdradius
if event.key == pygame.K_n: # n = Number the points drawn to screen with their index number.
if points < 5000:
drawnumbers = not drawnumbers
if event.key == pygame.K_o: # o = Outline around circles.
outlines = not outlines
if event.key == pygame.K_p: # p = Pressure changes like pounding surf waves, shaken in jar.
pressurechanges = not pressurechanges
if event.key == pygame.K_s: # s = Save new StartupXYZ.txt file, for later automatic restart.
xyzwrite('StartupXYZ.txt')
if event.key == pygame.K_w: # w = Wave signal draw to screen. Inputs turn violet, Outputs white.
drawwaves = not drawwaves
if not drawwaves and drawcolors == 0:
drawcolors = 1 # No blank screen.
if event.key == pygame.K_c: # c = change colors: None, Lennard-Jones, Number of connections.
drawcolors += 1
if drawcolors > 2:
if drawwaves: # Draw nothing while showing waves.
drawcolors = 0
else: # Else toggle to Lennard Jones so that a blank screen is not shown.
drawcolors = 1
if not paused:
myfont = pygame.font.SysFont('arial', 25)
# Pressure changes as in shaking in a jar of water, or caught inside pounding shoreline waves.
vesiclestartR = int(screenradius) # Start size of vesicle sphere, from calculated full size.
vesicleR = vesiclestartR * adjustR
cutoffR = forceCutOff * spacing * adjustR # Distance two points are too far away to exert force.
groupR = groupCutOff * spacing * adjustR
connectR = spacing * connectCutOff * adjustR # Distance two points are touching, bumping each other.
spherewidth = vesiclestartR * 2
ljdiff = 0
# Update positions of points, and rotation to drawn on screen.
if forces:
if pressureChanges:
adjustR, pressureRampDirection = pressureRamp(adjustR, pressureRampDirection)
# adjustR = 1 + random.uniform(-pressure, 0) # Pressure changes of churning water, shaking.
else:
adjustR = 1 # Stay full radius.
cutoffupdate += 1
if cutoffupdate > cutoffcalctime:
cutoffupdate = 0
groupupdate += 1
if groupupdate > groupcalctime:
groupupdate = 0
update_grouplist()
avcutoffs = update_cutofflist()
ljmaxwas = ljmax
ljavr, ljmin, ljmin_i, ljmin_j, ljmax, ljmax_i, atmaxforce = update_positions()
if atmaxforce != 0:
pressureRampDirection = 1
else:
pressureRampDirection = -1
if ljmax < 10: # If activity is low then increase Lennard Jones force.
dt *= 1.1
if dt >= .5:
dt = .5
else: # else decrease.
dt *= .95
if dt <= .001:
dt = .001
# if ljavr < .5: # If activity is low then increase Lennard Jones force.
# dt *= 1.1
# if dt >= .5:
# dt = .5
# else: # else decrease.
# dt *= .95
# if dt <= .001:
# dt = .001
update_connections()
update_angles()
cutoffcalctime = int(1000 / ljmax) * 5 # Timesteps before recalculating distances.
# cutoffcalctime = int(1 / (np.sqrt(abs(ljavr)) / 5) + 1) # Timesteps before recalculating distances.
# cutoffcalctime = int(2 / ljavr) * 5 # Timesteps before recalculating distances.
if cutoffcalctime > 200:
cutoffcalctime = 200
print(timestepnumber, " LJavr=", ljavr, " LJmin=", ljmin, " LJmax=", ljmax,
" dt=", dt, " r=", vesicleR, "of", screenradius)
# Popping out under pressure.
indexmin = ljsum.index(min(ljsum))
ljdiffwas = ljdiff
ljdiff = abs(ljmax - ljmaxwas)
if ljdiff < .00002:
print(ljdiff)
popwait += 1
if ljmax > 10: # If there is enough force on any one point.
pressureRampDirection = -1
if popwait > popwaitcount: # If there has been enough time since last swap
# Pop out the particle with most force around it
popwait = -1 # start new count for waiting again before swapping again
indexmax = ljsum.index(max(ljsum))
xpop = xp[indexmax]
ypop = yp[indexmax]
zpop = zp[indexmax]
azm, ele, rad = vectorsum(indexmin, cutofflist[indexmin])
xo, yo, zo = sphere2cart(azm, ele, rad)
xnew = xp[indexmin] + xo
ynew = yp[indexmin] + yo
znew = zp[indexmin] + zo
xp[indexmax] = xnew
yp[indexmax] = ynew
zp[indexmax] = znew
print(timestepnumber, '>>>>>> One popped out and was sent to one can pop in. <<<<<<')
myfont = pygame.font.SysFont('arial', 26)
textsurface = myfont.render('Pop-Out, 1 Moved', False, White)
screen.blit(textsurface, (height * .8, height * .84)) # height is half the width of screen
myfont = pygame.font.SysFont('arial', 22)
textsurface = myfont.render('From Blue to White Circle', False, White)
screen.blit(textsurface, (height * .76, height * .90)) # height is half the width of screen
pops =+ pops
cutoffupdate = cutoffcalctime # set to next update list of points in cutoff range.
groupupdate = groupcalctime # and set to next update list of points in larger groups.
# Propagate signals.
timestepsignals()
# Label screen:
myfont = pygame.font.SysFont('arial', 27)
textsurface = myfont.render('Front', False, White)
screen.blit(textsurface, (height * .83, height/10)) # height is half the width of screen
textsurface = myfont.render('Back', False, White)
screen.blit(textsurface, (height * 1.06, height/10))
myfont = pygame.font.SysFont('arial', 18)
textsurface = myfont.render(str(points) + ' mol/cell/column', False, White)
screen.blit(textsurface, (5, 0))
textsurface = myfont.render('LJmax: ' + str(round(ljmax, 4)), False, White)
screen.blit(textsurface, (5, height * .04))
textsurface = myfont.render('On: ' + str(signalactivity), False, White)
screen.blit(textsurface, (5, height * .08))
myfont = pygame.font.SysFont('arial', 20)
textsurface = myfont.render('TimeStep ' + str(timestepnumber), False, White)
screen.blit(textsurface, (height * .75, 0))
myfont = pygame.font.SysFont('arial', 18)
textsurface = myfont.render("LJ-Forces=" + str(forces) + " HoldRadius=" + str(holdradius) + " Pressure=" + str(pressureChanges), False, White)
screen.blit(textsurface, (5, height - 24))
if drawcolors == 2:
textsurface = myfont.render('Connections: Red=2 Orng=3 Yell=4 Grn=5 Blu=6 Viol=7 Gray=8', False, White)
if drawcolors == 1:
textsurface = myfont.render('Lennard Jones Potential, White = Extreme Pressure', False, White)
if drawcolors == 0:
textsurface = myfont.render('Shown as Points', False, White)
screen.blit(textsurface, (height, height - 24))
myfont = pygame.font.SysFont('arial', 18)
textsurface = myfont.render('Gary Gaulin, 2023', False, White)
screen.blit(textsurface, (height * 1.7, 4)) # Height is half the width of screen
render() # Draw to screen.
pygame.display.update # Display screen that was being drawn, to monitor display.
# pygame.image.save(screen, 'Images/' + str(timestepnumber) + '.png') # Save image file for making gif video.
fps.tick(50)
screen.fill(Black) # Screen can now be cleared for next draw.
# Rotate sphere a small amount for next draw.
rotation += rotationstep # Add screen rotation amount, for next rendering.
if rotation > radian: # Keep rotation amount in range from 0 to 1 Radian.
rotation = 0
timestepnumber += 1
Wednesday, October 25, 2023
Create Vesicle using Lennard Jones Potential, and propagate Traveling Waves over surface
Subscribe to:
Posts (Atom)
