The program repeatedly goes from a squeeze cycle to align points into symmetric as possible geometry, and expansion cycle:
I’m still new to Python, so please forgive where it may be obvious that I need to figure out how to most simply improve array structure and other things. Coding suggestions are appreciated:
# Form stable vesicle membrane around sphere, Gary Gaulin, 8/4/2019
# Uses Lennard Jones Potential for attraction and repulsion behavior,
# representative of polar molecules, bubbles, and cortical cell migration.
# Starts with almost enough points to fill sphere.
# There are then large holes and unfilled 5 sided pores.
# For cortical stem cells these pores are a space for them to fill, not be any.
# Radius is then slowly decreased until pressure causes chain reaction of realignments,
# towards most stable configuration, ideal symmetry, for given number of points.
import sys
import pygame
import random
import numpy as np
Pi = 3.14159265358979323846264338327950288412
HalfPi = Pi / 2
Radian = Pi * 2
ThreeQtrRadian = Radian * 3/4
# Startup Parameters
# Example of 392 is the dimple count of a vintage Hogan golfball.
# It has an unusually uniform dimple uniformity, as though made by packing real balls around it.
# Other examples are:
# points = 12 # https://en.wikipedia.org/wiki/Regular_icosahedron
# points = 32 # https://en.wikipedia.org/wiki/Pentakis_dodecahedron
points = 392 # Total number of (cells, molecules, etc.) points.
screensize = Width, Height = (640, 640)
sr = Height*.38 # Screen radius, with added space around sphere.
sas = 4*Pi*(sr*sr) # Area of sphere at this given screen radius.
sar = sas/points # Divide for area of one Lennard-Jones radius, between points.
r = np.sqrt(sar/Pi)*2 # Radius between points, given area.
spacing = r-((1/points)*r) # Radius between points, ideal spacing.
spheremaxradius = int(sr * 1.2) # Maximum radius of sphere, points are located around.
pressuremax = 20 # Maximum amount of average Lennard-Jones force to apply.
pressurestep = 1 # Amount of radius change per time step, in screen pixels.
pressuretime = 200 # Time step duration to hold under pressure before releasing.
dt = .01 # Timestep size, how (slower but) precise it moves.
maxvelocity = spacing/10 # Velocity limit for when enough force to fly off sphere.
cutoffdist = spacing*3 # Distance two points are too far away to exert force.
connectdist = spacing*1.15 # Distance two points are touching, bumping each other.
cutoffcalctime = 50 # Time steps between updates of i,j points within cut-off range.
rotationstep = .05 # Amount of rotation per timestep.
rotation = 0 # Initialize rotation angle, to draw all points on screen.
timestepnumber = 0
epsilon = 1 # Lennard Jones variable, remains 1, normalized.
sigma = 1 # Lennard Jones variable, remains 1, normalized.
Blue = (0, 0, 150)
Black = (0, 0, 0)
Green = (0, 100, 0)
Yellow = (120, 120, 0)
# Initialize x,y,z arrays for points on sphere.
xp = [0 for _ in range(points)]
yp = [0 for _ in range(points)]
zp = [0 for _ in range(points)]
# Initialize Lennard-Jones sum array.
ljsum = [0 for _ in range(points)]
# X,Y of points on screen.
xyzs = [0 for _ in range(points)]
# Periodically updated list of points within cutoff range.
near = [0 for _ in range(points*100)]
# Periodically updated list of connections between points.
connect = [0 for _ in range(points*20)]
# Cartesian to Spherical coordinates.
def cart2sphere(x, y, z):
hxy = np.hypot(x, y)
el = np.arctan2(z, hxy)
az = np.arctan2(y, x)
return az, el, 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
def drawaxisline(x, y, z, c):
az, el, r = cart2sphere(x, y, z) # Cartesian to Spherical (az, el, r)
x1, y1, z1 = sphere2cart(az + rotation, el, r) # Spherical to Cartesian.
pygame.draw.line(screen, c, (xc, yc), (x1+xc, z1+yc), 1)
def findnearby():
ni = 0
ci = 0
for i in range(points): # For all the points in array
for j in range(points): # Except itself
if i != j:
dx = xp[j] - xp[i] # Calculate x,y,z distances.
dy = yp[j] - yp[i]
dz = zp[j] - zp[i]
# Use the three distances for calculating an always positive 3D distance.
# This is the very time consuming command to less occasionally use.
dist = np.sqrt((dx * dx) + (dy * dy) + (dz * dz))
# If distance is within cutoff range then save 'j' point, at this 'i'.
if dist < cutoffdist:
near[ni] = i, j
ni += 1
# If distance is within touching range then save 'j' point, at this 'i'.
if dist < connectdist:
connect[ci] = i, j
ci += 1
print(ni, ci)
return ni, ci
def update_positions():
# 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 , and maximum
ljt = 0 # Initialize variable for totaling Lennard-Jones force.
ljmax = -10000 # Initialize maximum variable, negative when points attract.
# For all the points in array.
for n in range(ni):
i,j = near[n]
# 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 ljr=1 force is 0.
if ljr > 0: # Prevent divide by zero error when LJ radius is zero.
lj = 4 * epsilon * ((sigma/ljr)**12 - (sigma/ljr)**6)
else:
lj = 0
# Update monitor variables for totaling average and maximum Lennard-Jones reading..
if lj > ljmax:
ljmax = lj
ljt += lj # Add to sum for average of all points.
ljsum[i] += lj # Used for brightness to draw on screen.
# Reduce timestep for smooth movement, fewer points wildly bouncing off each other.
lj = lj * dt
# Limit amount of displacement/velocity as though through water, or cell migration, not frictionless vacuum.
if lj > maxvelocity:
lj = maxvelocity
if lj < -maxvelocity:
lj = -maxvelocity
# Use distance for calculating spherical azimuth and elevation angle r, as in cart2sphere.
el = np.arctan2(dz, dist)
az = np.arctan2(dy, dx)
# Use Lennard Jones amount to calculate displacement amount, for given force.
x, y, z = sphere2cart(az, el, lj)
# 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
# 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]
return ljt/ni, ljmax
def render():
screen.fill(Black)
spherewidth = sphereradius * 2
# 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)
x, y, z = sphere2cart(az + rotation, el, r) # Spherical to Cartesian.
# Vary brightness by amount of force, positive increasingly orange, negative blue.
poscolor = 0
negcolor = 0
if ljsum[i] > 0:
poscolor = int(ljsum[i])
if poscolor > 255:
poscolor = 255
if ljsum[i] < 0:
negcolor = int(-ljsum[i] * 50)
if negcolor > 255:
negcolor = 255
c = (poscolor, (poscolor + negcolor) * .6, negcolor)
xyzs[i] = x, y, z, c
# To properly draw on top of each other sort the list according to what is furthest away, [1] = by Y axis.
xyzs.sort(key=lambda x: x[1])
# Draw all points to screen.
for i in range(points):
# Calculate x,y and radius then draw point to screen buffer.
x, y, z, c = xyzs[i]
xs = int(x + xc)
ys = int(z + yc)
d = ((y+sphereradius)/spherewidth*.4) + .6 # Fractional amount, to 1, according to how far away, depth.
if d > 1:
d = 1
cr = int(d * spacing * .45) # Circle radius.
pygame.draw.circle(screen, c, [xs, ys], cr, 0)
pygame.draw.circle(screen, (255, 255, 255), [xs, ys], cr, 1)
# Draw axis lines.
ra = sphereradius * 1.1 # Distance of axis line, from 0.
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)
pygame.display.update() # Display everything drawn into screen buffer.
fps.tick(60)
# Fill Environment arrays with randomly placed points.
# They will become located on surface of sphere.
for i in range(points):
xp[i] = random.randrange(-100, 100) # Fill positive half of cube.
yp[i] = random.randrange(-100, 100)
zp[i] = random.randrange(-100, 100)
fnd = cutoffcalctime # Less than full amount when starting off.
sphereradius = spheremaxradius # Start with maximum size sphere.
dialog = "Squeezing"
holdcount = pressuretime
# PyGame Initialization
pygame.init()
screen = pygame.display.set_mode(screensize)
xc = int(screensize[0]/2)
yc = int(screensize[1]/2)
pygame.display.set_caption("Spherical Membrane. Positive Axis: x=green, y=blue, z=yellow")
fps = pygame.time.Clock()
paused = False
while True:
for event in pygame.event.get():
if event.type == pygame.QUIT:
pygame.quit()
sys.exit()
if event.type == pygame.KEYUP:
if event.key == pygame.K_SPACE:
paused = not paused
if not paused:
timestepnumber += 1
# Keep points on sphere by converting to Spherical then (using sphere radius) to Cartesian again.
for i in range(points): # 1D Radial distance between x points.
az, el, r = cart2sphere(xp[i], yp[i], zp[i]) # Cartesian to Spherical.
xp[i], yp[i], zp[i] = sphere2cart(az, el, sphereradius) # Spherical to Cartesian again.
# Update positions of points, and rotation to drawn on screen.
fnd += 1
if fnd > cutoffcalctime:
fnd = 0
ni, ci = findnearby()
ljsum = [0 for _ in range(points)]
ljavr, ljmax = update_positions()
render()
rotation += rotationstep # Add to screen rotation amount.
if rotation > Radian: # Keep in range from 0 to 1 Radian
rotation = 0
# Pressurize and unpressurize by changing radius of sphere, by average L-J force amount.
moreless = "" # Clear optional "More" or "Less" dialog.
if dialog == "Squeezing":
holdcount -= 1 # Decrement squeezing time counter.
if holdcount < 0: # Time to recalculate points within cutoff range?
dialog = "Begin Expansion LJ=0" # Changing 'dialog' ends Squeezing.
if ljavr <pressuremax:
sphereradius -= pressurestep # Decrease sphere radius, increase pressure.
moreless = "More, " # Add to dialog as "Squeezing More".
else:
sphereradius += pressurestep # Increase sphere radius, decrease pressure.
moreless = "Less, " # Add to dialog as "Squeezing Less".
else:
if ljavr < 0: # Fully expanded?
# Save file here.
dialog = "Squeezing" # Changing dialog to "Squeezing" ends expansion.
holdcount = pressuretime # Start new countdown to next expansion cycle.
else:
dialog = "Expanding" # Otherwise use "Expanding" for dialog instead.
sphereradius += pressurestep # Increase sphere radius, decreases pressure.
print(timestepnumber, " ", dialog, moreless, "Radius =", int(sphereradius), " LJaverage =", ljavr, " LJmax =", ljmax)My next challenge is to use this to articulate a network, to attempt propagating traveling waves in the usual simple way I use, around a 3D brain ball or surface with some curvature as in biological membranes.
And credit to my wife Laurie, for counting all the dimples in ten’s using a magic marker, while I worked on the code to recreate the pattern. It’s sometimes hard to force out the influence of a small amount of space, golf ball dimple mould engineers could have in other ways squeezed out too.
