Modeling and animating the evolution of a cinder cone in Python

Numerically modeling of the evolution of a cinder cone where individual rock particles or "bombs" are launched from a central point with "random" trajectories from a central location and deposition/accumulation. Each process is explicit, eventually this process develops the shape of a classic cinder cone. 

I wrote this as part of a graduate course and added comments throughout so it is easy to follow what is happening. Also included is an example of how to save images of the process and combine them into a gif animation. 

Scientific and numerical modeling in Python is best accomplished with the Numpy library. This book is a great resource if you are interested in improving your numerical methods in Python using Numpy, SciPy, and matplotlib. 

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import imageio
from pathlib import Path
import re

Define physical parameters of bomb size and trajectories, determined from literature and trial and error.

## Bomb velocities, m/s
vel_min = 32
vel_max = 76 
## trajectory angles, degrees
ang_min = 345
ang_max = 375 
## diameter range, in meters
dia_min = .1 
dia_max = 0.9

Functions for creating, launching, and depositing bombs and defining the spatial domain.

######
##  bomb creation
######
def makeBomb(vMin,vMax,aMin,aMax,dMin,dMax,num=1):
    """Inputs are velocity (m/s), angle (degrees), and diameter (m) min and max values 
    outputs are particle resultant velocity, angle(rad), and area (sq m)"""
    vel = (vMax - vMin)*np.random.random(num) + vMin  ## velocity in m/s 
    ang = (aMax - aMin)*np.random.random(num) + aMin
    ang *= np.pi/180  ## return angle in radians
    dia = (dMax - dMin)*np.random.random(num) + dMin
    are = np.pi*(dia/2.)**2  ## return area in square meters
    return {'velocity':vel,'angle':ang,'area':are} ## returned as dict for readability

## x spatial grid for simulation and plotting
x_grid = np.linspace(-1000,1000,2001) # unit meter grid

######
##  bomb launching 
######
def launchBomb(vel,ang,z_0=0,x=x_grid):
    """ Inputs: velocity (vel), angle of launch in rads (ang), launch elev (z_0), 
    and x grid that is defaulted to use a predifined grid in main program
    return x and z position vectors """ 
    g = 9.81 ## m/s^2
    x_vel = vel * np.sin(ang)
    z_vel = np.abs( vel * np.cos(ang))
    z = np.zeros_like(x)
    z = (z_vel * x)/x_vel - (1/2.) * g * (x/x_vel)**2 + z_0
    return  z  ## elevation position array
    
######
##  bomb depsiting
######
def depositBomb(x_traj,z_traj,elev_profile,bombArea):
    """Inputs: x distance grid, z (elevation) coordinates for a trajectory, previous 
    land elevation profile and bomb area. Finds where bomb hits ground and 
    adds accumulated mass to the elevation profile and returns it"""
    vent_ind = np.where(x_traj == 0)[0][0] ## get index of vent (middle of x grid)
    
    ## check if lauch went to the right
    if np.where(z_traj == np.max(z_traj)) > vent_ind:
        x = np.where(z_traj[vent_ind:] <= elev_profile[vent_ind:])[0]  ## where trajectory below elevation
        deposit_ind = np.min(x[1:]) + vent_ind -1  ## get the minimum index of occurence (first intersection with land)
        ## redeposit on cliffs 
        while elev_profile[deposit_ind] - elev_profile[deposit_ind + 1] > np.sin(np.pi/6):  ## forward roll
            deposit_ind += 1
        while elev_profile[deposit_ind] - elev_profile[deposit_ind - 1] > np.sin(np.pi/6):  ## backward fall
            deposit_ind -= 1  
            ## prevent indefinate rolling if steep in both directions 
            if elev_profile[deposit_ind] - elev_profile[deposit_ind+1] > bombArea:
                break        
    ## check if launch went to the left of vent    
    elif np.where(z_traj == np.max(z_traj)) < vent_ind:
        deposit_ind = np.max(np.where(z_traj[:vent_ind] <= elev_profile[:vent_ind])[0]) ## 
        ##redeposit particles on cliffs
        while elev_profile[deposit_ind] - elev_profile[deposit_ind - 1] > np.sin(np.pi/6): # forward
            deposit_ind -= 1  
        while elev_profile[deposit_ind] - elev_profile[deposit_ind + 1] > np.sin(np.pi/6): # backwards
            deposit_ind += 1
            ## Prevent infinite movement if steep on both directions
            if elev_profile[deposit_ind] - elev_profile[deposit_ind-1] > bombArea:
                break                
    ## check if bomb landed in the vent itself
    else:
        deposit_ind = vent_ind ## need to fix redistribution in the vent, prevent huge buildup 
        if elev[deposit_ind] > elev[deposit_ind+1] or elev[deposit_ind] > elev[deposit_ind-1]:
            deposit_ind += np.random.choice(np.arange(-15,15,1))  ## create some chaos to redistribute vent deps
            elev_profile[vent_ind-5:vent_ind+5] -= 0.1 ## cause meltdown of vent occasionally to prevent filling in
    ## increment elevation profile at the correct index by area of bomb
    elev_profile[deposit_ind] += bombArea
    return elev_profile

Combining semi-random bomb creation, ejection, and deposition in a loop. This is also an example of how you can plot and save a series of images and then combine them to make an animation of the physical simulation.

######
## Launch Zone
######

## initialize elevation profile using zero grid as input
bomb = makeBomb(vel_min,vel_max,ang_min,ang_max,dia_min,dia_max)
traj = launchBomb(bomb['velocity'],bomb['angle'])
elev = depositBomb(x_grid,traj,np.zeros_like(x_grid),bomb['area'])

## number of eruptions and how often to show
ndiv = 5000 
nerupts = 400000

## for final figure, not movie
fig = plt.figure(figsize=(12,9))
ax = fig.add_subplot(111)

## loop through number of eruptions simulating cinder cone development
for i in range(nerupts):
    bomb = makeBomb(vel_min,vel_max,ang_min,ang_max,dia_min,dia_max)
    traj = launchBomb(bomb['velocity'],bomb['angle'],z_0=elev[np.where(x_grid == 0)[0]])
    elev = depositBomb(x_grid,traj,elev,bomb['area'])
    if (i % ndiv == 0):
        plt.plot(x_grid,elev) ## plot updated elevation profile
        plt.xlim(-200,200)
        plt.show()
        plt.pause(0.2)
        ## will plot current trajectory 
        plt.plot(x_grid[traj>elev],traj[traj>elev], 'r.', ms = 1.5)  
        plt.savefig(f'fig{i}.png') # for animation

## plot final elevation profile with every ndiv profile and trajectory overlain
plt.plot(x_grid,elev) 
plt.xlim(-400,400)
plt.ylim(0,650)
plt.ylabel('Height (m)')
plt.title('Cinder cone evolution after %i eruptions,\n deposits and trajectory shown every %i th eruption' % (nerupts,ndiv), fontsize = 16)
plt.savefig('Cinder cone evolution after {} eruptions.png'.format(nerupts), dpi = 600)

Gather the saved images and use a natural (human) sorting to order them from start to finish then save as a gif animation with the imageio Python library. 

def natural_key(string_):
    """See https://blog.codinghorror.com/sorting-for-humans-natural-sort-order/"""
    return [int(s) if s.isdigit() else s for s in re.split(r'(\d+)', string_)]

images = sorted(list(map(str, Path().glob('*.png'))), key=natural_key)
images = [x for x in images if not 'Cinder' in x]
images = [imageio.imread(i) for i in images]
imageio.mimsave('eruption.gif',images,fps=3)

And now for the fun part, click on the video below to watch the resulting animation of cinder cone evolution and bomb trajectories plotted every 5000 eruptions: