Skip to content

Instantly share code, notes, and snippets.

# jstults/sears_haack.py

Created February 13, 2012 02:04
Star You must be signed in to star a gist
Script to generate a Sears-Haack body for a mini-Estes motor
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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()
to join this conversation on GitHub. Already have an account? Sign in to comment