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,  = by Y axis. xyzs.sort(key=lambda x: x) # 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/2) yc = int(screensize/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.