Script to generate a Sears-Haack body for a mini-Estes motor

sears_haack.py

import scipy as sp
from scipy.integrate import ode 
from scipy.optimize import fsolve 
from scipy.optimize import fmin 
from matplotlib import rc
rc('text', usetex=True)
tick_size = 'large'
label_size = 'xx-large'
rc('xtick', labelsize=8)
rc('ytick', labelsize=8)
rc('legend', fontsize=12)
rc('axes', labelsize=12)
fig_width_pt = 469.75502 
inches_per_pt = 1.0/72.27
golden_mean = (sp.sqrt(5.0)-1.0)/2.0
fig_width = fig_width_pt *inches_per_pt 
fig_height = fig_width * golden_mean
fig_size = [fig_width, fig_height] 
rc('figure', figsize=fig_size)

from matplotlib import pylab as plt

def sh_radius(rmax, x):
    return(rmax*(4*x*(1-x))**(3.0/4.0))
def sh_area(rmax,x):
    return(sp.pi*rmax**2*(4*x*(1-x))**(3.0/2.0))
def sh_vol(rmax,L):
    return(3*sp.pi**2*rmax**2*L/16.0)
def circ_area(d):
    return(sp.pi*d**2/4.0)

def r_given_l(rmax,D,eL,Lref):
    """function for zero finding the rmax that gives the r at 0.5+eL/2/Lref"""
    x = 0.5 + eL/2.0/Lref
    return(sh_radius(rmax,x)-D/2.0)

def find_min_v(L,D,eL):
    rmax = fsolve(r_given_l, D*2, (D,eL,L))
    return(sh_vol(rmax,L))

D = 13 # mm
eL = 44 # mm
Lref = sp.linspace(1.3*eL, 4*eL, 40) # mm 
nx = 1000
x = sp.linspace(0, 1, nx)

rmax = sp.zeros(Lref.shape)
for i,L in enumerate(Lref):
    rmax[i] = fsolve(r_given_l, 20, (D, eL, L))

vol = map(sh_vol, rmax, Lref)
minvi = sp.argmin(vol)

L_minv = fmin(find_min_v, Lref[minvi], (D,eL))
r_minv = fsolve(r_given_l, rmax[minvi], (D,eL,L_minv))
minv = sh_vol(r_minv,L_minv)

x1 = sp.array([2.250,-1.697,0.000])
x2 = sp.array([15.890,-1.697,0.000])

for i,L in enumerate(Lref):
    plt.plot(L*x, sh_radius(rmax[i],x))
plt.plot(sp.linspace(0,Lref.max(),50), D/2.0*sp.ones(50), 'k--')
plt.plot(L_minv*x, sh_radius(r_minv,x), 'k-', lw=4)
plt.xlabel(r'$x$')
plt.ylabel(r'$r$')
plt.savefig("sears_haack_variation.png")

plt.figure(figsize=(18,2))
plt.plot(L_minv*x, sh_radius(r_minv,x), 'k-', lw=4)
plt.axis('scaled') 
plt.savefig("shr.svg")

plt.figure()
plt.plot(Lref, vol)
plt.plot(L_minv, minv, 'bo')
plt.xlabel(r'length')
plt.ylabel(r'volume')
plt.savefig("sears_haack_length_volume.png")

plt.show()