### 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))

## 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: