Skip to content

Instantly share code, notes, and snippets.

@jstults
Created February 13, 2012 02:04
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save jstults/1812752 to your computer and use it in GitHub Desktop.
Script to generate a Sears-Haack body for a mini-Estes motor
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment