From 6bf468b1d146f5a8282c99d9deee6e43df546236 Mon Sep 17 00:00:00 2001 From: Joe Anderson Date: Sun, 30 Oct 2011 15:19:14 -0500 Subject: Initial State --- throw.py | 326 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 326 insertions(+) create mode 100644 throw.py diff --git a/throw.py b/throw.py new file mode 100644 index 0000000..646dea9 --- /dev/null +++ b/throw.py @@ -0,0 +1,326 @@ + +from time import time as current_time +import sys +import numpy +from scipy.interpolate import interp1d +from scipy.integrate import odeint +from scipy.optimize import fmin +from matplotlib import pyplot + +global time0 +time0 = current_time() + +tau = 2*numpy.pi +def taylor_sin(a): + x = a % tau + r = (x + -x**3/6 + +x**5/120 + -x**7/5040 + +x**9/362880 + -x**11/39916800 + +x**13/6227020800 + -x**15/1307674368000 + +x**17/355687428096000 + -x**19/121645100408832000 + +x**21/51090942171709440000 + -x**23/25852016738884976640000 + +x**25/15511210043330985984000000 + -x**27/10888869450418352160768000000 + #+x**29/8841761993739701954543616000000 + #-x**31/8222838654177922817725562880000000 + #+x**33/8683317618811886495518194401280000000 + #-x**35/10333147966386144929666651337523200000000 + #+x**37/13763753091226345046315979581580902400000000 + #-x**39/20397882081197443358640281739902897356800000000 + ) + return r +def taylor_cos(a): + x = a % tau + r = (1 + -x**2/2 + +x**4/24 + -x**6/720 + +x**8/40320 + -x**10/3628800 + +x**12/479001600 + -x**14/87178291200 + +x**16/20922789888000 + -x**18/6402373705728000 + +x**20/2432902008176640000 + -x**22/1124000727777607680000 + +x**24/620448401733239439360000 + -x**26/403291461126605635584000000 + #+x**28/304888344611713860501504000000 + #-x**30/265252859812191058636308480000000 + #+x**32/263130836933693530167218012160000000 + #-x**34/295232799039604140847618609643520000000 + #+x**36/371993326789901217467999448150835200000000 + #-x**38/523022617466601111760007224100074291200000000 + ) + return r + +#### +# The interpolation is slower, uses more RAM, +# and is two orders of magnitude less accurate. +# +#scale = 300. +#angles = numpy.linspace(0, tau, scale) +#_cos = numpy.cos(angles) +#def floor_cos(a): +# return _cos[int((a % tau)/tau * scale)] +#def floor_sin(a): +# return floor_cos(a - tau/4) +# +#interpolate_cos = interp1d(angles, numpy.cos(angles)) +#interpolate_sin = interp1d(angles, numpy.sin(angles)) +#def interp_cos(a): return interpolate_cos(a % tau) +#def interp_sin(a): return interpolate_sin(a % tau) + +### trig approximations +cos = taylor_cos +sin = taylor_sin + +### radii: +global r0 +r0 = 1. #meters -- the radius of the first full-circle rotation +r_anchor = .2 #meters -- the radius of the anchor point +r_winch = .2 #meters -- the max radius of the roll of cable +length = 1000 #meters -- the length of the cable +phi_max = 2 * (length - r0) / r_winch #rad +#r = lambda phi: r0 + r_winch*phi # the radius of the satellite +r = lambda phi: r0 + r_winch*(phi - phi**2/(2*phi_max)) + +### gravitational acceleration: +### the deprecated approximation g = 9.71 m/s^2 +### factoring in the r-phi-dependence +geopotential = 3.987 * 10**14 #cubic meters per second^2 +earth_radius = 6.3675 * 10**6 #meters +balloon_height = 40*10**3 +height = lambda phi: balloon_height + r(phi) * cos(phi) +g = lambda phi: geopotential / (earth_radius + height(phi))**2 + +### linear friction can be neglected due to: +### low viscosity +### high velocity => high Reynold's number +### the formulation below neglets pressure differences. +# for standard air: +#visc0 = .01827 # reference viscosity +#T0 = 291.15 # reference temperature +#suth_C = 120 # Sutherland's constant +# for pure diatomic hydrogen: +#visc0 = .00876 +#T0 = 293.85 +#suth_C = 72 +#a = 0.555 * T0 + suth_C +#b = lambda T: 0.555 * T + suth_C +#visc = lambda T: visc0 * (a/b(T)) * (T/T0)**(3/2) +#visc = visc(300) # standard temperature -- inaccurate +#lin_co = 3 * numpy.pi * visc * Diameter # for a sphere +# fric_lin = lin_co * r * phi1 + +### quadratic friction: +#molar_mass = .02897 #kg/mol -- for air up to +#ideal_gas_const = 8.314 #m^3 Pa / (mol K) +#pressure = 500 #Pa +#Temp = 300 #K +#density = pressure * molar_mass / (ideal_gas_const * Temp) #kg/m^3 +density = .003 #kg/m^3 +kappa = .25 # for a sphere +radius = .1 #meters -- radius of the satellite cage +Area = numpy.pi * radius **2 # for a sphere +quad_co = kappa * density * Area +# fric_quad = quad_co * (r * phi1 * vcm_co)**2 + +### mass: +global m, Mass +m = 39.7 #kilograms -- mass of the satellite cage (overestimate) +#l = 0.2 #kilograms/meter -- density of the cable (20mm: 41.9 kN) +#l = 0.52 #kilograms/meter -- density of the cable (32mm: 105.0 kN) +l = 0.16 #kilograms/meter -- density of Super Max cable, 16mm: 269.8kN +max_tension = 269.8 #kN, given by choice of cable above. +Mass = 670 #kilograms -- mass of /everything/ +M = lambda phi: Mass - m - l*r(phi) +v_cm_co = lambda phi: (.5 * l*r(phi) + M(phi)) / (m + l*r(phi) + M(phi)) + +tan = lambda phi: r_anchor / r(phi) # tan(theta) +def throw(phi, time): + """The diff eqs for the windup.""" + phi, phi1 = phi + prime_vars = [ + # phi' = phi' + phi1, + # phi'' = + ( tan(phi) - r_winch/r(phi) ) * phi1**2 + + g(phi) * ( sin(phi) - cos(phi)*tan(phi) ) / r(phi) + - quad_co * r(phi) *(v_cm_co(phi) * phi1 )**2 / m + ] + return prime_vars + +#good final times: 20 185 1874 +def solve_fun(func, init_vars, space=(0, 1874, 10**5)): + time0 = current_time() + time = numpy.linspace(*space) + soln = odeint(func, init_vars, time) + print "took %s seconds to solve" % (current_time() - time0) + return time, soln + +def interpret(soln, max_tension=max_tension): + phi = soln[:, 0] + phi1 = soln[:, 1] + v_tan = v_cm_co(phi) * phi1 * r(phi) / 1000 #kilometers/second + cos_t = lambda phi: numpy.cos(numpy.arctan(tan(phi))) # cos(theta) + tension = (m/cos_t(phi))*(phi1**2 * r(phi) - g(tau/4))/1000 #kiloNewtons + + try: + maxindex = numpy.argmax(tension >= max_tension) + max_arg = numpy.argmax(v_tan[:maxindex]) + except ValueError: max_arg = numpy.argmax(v_tan) + + return phi, phi1, v_tan, cos_t, tension, max_arg + +def plot(time, soln, graphs=['all'], i=5): + phi, phi1, v_tan, cos_t, tension, max_arg = interpret(soln) + phi_max, phi1_max, v_tan_max, tension_max = \ + map(lambda x: numpy.ndarray.__getitem__(x, max_arg), + [phi, phi1, v_tan, tension]) + + def test(*args): + b = 0 + for string in args: + b += string in graphs or \ + ('all' in graphs and "!" + string not in graphs) + return bool(b) + + i -= 1 + if test('gravity'): + pyplot.figure(i) + pyplot.plot(time, g(phi)) + pyplot.grid(True) + pyplot.ylabel("gravitational acceleration (m/s^2)") + pyplot.xlabel("time (seconds)") + + i -= 1 + if test('phi', 'r', 'radius'): + pyplot.figure(i) + pyplot.plot(time, phi / tau, label="phi (cycles)") + pyplot.plot(time, r(phi), label="r (m)") + pyplot.grid(True) + pyplot.legend(loc="upper left") + pyplot.title("phi and radius vs. time") + pyplot.ylabel("phi (cycles) and radius (meters)") + pyplot.xlabel("time (seconds)") + + i -= 1 + if test("phi'", "omega", "phi1"): + pyplot.figure(i) + pyplot.plot(time, phi1 / tau) + pyplot.grid(True) + pyplot.ylabel("phi' (Hertz)") + pyplot.xlabel("time (seconds)") + + i -= 1 + if test('v_tan', 'tension'): + pyplot.figure(i) + rect = [.1, .1, .8, .8] + + ax1 = pyplot.axes(rect) + ax1.yaxis.tick_left() + ax1.xaxis.grid(True) + pyplot.plot(phi / tau, v_tan, marker='+', markevery=len(phi)/50) + pyplot.ylabel("+: tangential velocity (km/s)") + pyplot.xlabel("phi (cycles)") + pyplot.ylim(ymin=0) + + ax2 = pyplot.axes(rect, frameon=False, sharex=ax1) + ax2.yaxis.tick_right() + ax2.yaxis.set_label_position('right') + pyplot.setp(ax2.get_xticklabels(), visible=False) + pyplot.plot(phi / tau, tension, marker='x', markevery=len(phi)/50) + pyplot.ylabel("x: Tension (kN)") + pyplot.ylim(ymin=0) + + pyplot.axvline(x=phi_max/tau) + + #import pdb; pdb.set_trace() + +def opt_fun(guesses, f, func, init_vars): + """Solve, plot, and compare the diff eq.""" + global r0, m, Mass + r0, m, Mass = guesses + + time, soln = solve_fun(func, init_vars) + + phi, phi1, v_tan, cos_t, tension, max_arg = interpret(soln) + phi_max, phi1_max, v_tan_max, tension_max = \ + map(lambda x: numpy.ndarray.__getitem__(x, max_arg), + [phi, phi1, v_tan, tension]) + + f.write("%s, " % r0) + f.write("%s, " % m) + f.write("%s, " % Mass) + f.write("%s, " % phi_max) + f.write("%s, " % phi1_max) + f.write("%s, " % tension_max) + f.write("%s\n" % v_tan_max) + + # If you want to add limits, you can do so by changing the return to: + #return 1/v_tan + Mass/600 + phi1.max() + # or + #return 1/v_tan + 1000*(Mass > Mass_limit) + 1000*(phi1.max() > phi1_limit) + return 2/v_tan_max + Mass/600 + phi1.max() / tau + +v0 = numpy.sqrt(g(0) * r(0)) +w0 = v0/r0 +inits = [0, w0] + +#outfile = 'C:\\Users\\dave\\Desktop\\out' +outfile = '/home/jandew/Documents/Impossible Challenges/microsatellite/out' +def optimize(): + fil = file(outfile, 'w') + fil.write('r0, sat. mass, Total Mass, phi, omega, Tension, Max Vel.\n') + guess = (r0, m, Mass) + #try: print opt_fun(guess, fil, throw, inits) + try: print fmin(opt_fun, guess, args=(fil, throw, inits)) + finally: fil.close() + +def solve(): + time, soln = solve_fun(throw, inits) + plot(time, soln, ['all', '!gravity']) + pyplot.show() + +def compare(): + global m, Mass + m = 39.7 + Mass = 650 + time, soln = solve_fun(throw, inits) + plot(time, soln, ['tension'], 8) + pyplot.title("Mass = %i" % Mass) + + m = 39.7 + Mass = 670 + time, soln = solve_fun(throw, inits) + plot(time, soln, ['tension'], 7) + pyplot.title("Mass = %i" % Mass) + + m = 39.7 + Mass = 690 + time, soln = solve_fun(throw, inits) + plot(time, soln, ['tension'], 6) + pyplot.title("Mass = %i" % Mass) + + m = 41.81320863 + Mass = 4172 + time, soln = solve_fun(throw, inits) + plot(time, soln, ['tension'], 5) + pyplot.title("Mass = %i" % Mass) + + pyplot.show() + +if len(sys.argv) == 2: + func_name = sys.argv[1] + if func_name in dir(): + exec(func_name + '()') + +optimize() + -- cgit v1.2.3-54-g00ecf