From bfe5324fff7ab35b97220d7fdf9b1f9fb1a3369a Mon Sep 17 00:00:00 2001 From: Joe Anderson Date: Sun, 13 Nov 2011 10:06:27 -0600 Subject: deleting old file --- throw.py | 307 --------------------------------------------------------------- 1 file changed, 307 deletions(-) delete mode 100644 throw.py (limited to 'throw.py') diff --git a/throw.py b/throw.py deleted file mode 100644 index 4c7d691..0000000 --- a/throw.py +++ /dev/null @@ -1,307 +0,0 @@ - -from time import time as current_time -import sys, itertools -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 - -### trig approximations -cos = taylor_cos -sin = taylor_sin - -### radii: -global r0, r_roll, length -r0 = 1.151467 #meters -- the length of cable in the first full-circle rotation -r_anchor = .2 #meters -- the radius of the anchor point -r_roll = .2045575 #meters -- the max radius of the roll of cable -length = 578.63682 #meters -- the length of the cable -phi_max = 2 * (length - r0) / r_roll #rad -r = lambda phi: r0 + r_roll*(1 - phi/(2*phi_max))*phi # cable length let out -r_dot_co = lambda phi: r_roll*(1 - phi/phi_max) #*phi1 -r_sat = lambda phi: numpy.sqrt(r(phi)**2 + r_anchor) # radius of satellite - -### gravitational acceleration: -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 - -### quadratic friction: -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.16 #kilograms/meter -- density of Super Max cable, 16mm: 269.8kN -l = 0.41 #kilograms/meter -- density of Super Max cable, 26mm: 647.4kN -max_tension = 647.4 #kN, given by choice of cable above. -Mass = 670 #kilograms -- mass of /everything/ -M = lambda phi: Mass - m - l*r(phi) -r_cm = lambda phi: (m*r(phi) + .5*l*r(phi)**2) / (m + l*r(phi) + M(phi)) -v_cm_co = lambda phi: 1 - r_cm(phi)/r(phi) - -def throw(phi, time): - """The diff eqs for the windup.""" - phi, phi1 = phi - prime_vars = [ - # phi' = phi' - phi1, - # phi'' = - ( (r_sat(phi)*r_anchor/r(phi) - r_dot_co(phi)) * phi1**2 - + g(phi) * ( sin(phi) - cos(phi)*r_anchor/r(phi) ) - - quad_co * ( v_cm_co(phi) * r(phi) * phi1 )**2 / m - ) / r_sat(phi) - ] - return prime_vars - -#good final times: 20 185 1874 -def solve_fun(func, init_vars, timespace): - time0 = current_time() - time = numpy.linspace(*timespace) - soln = odeint(func, init_vars, time) - print " :: took %s seconds" % (current_time() - time0) - return time, soln - -TimespaceTooSmall = Warning(" *** Requires larger timespace!!! *** ") -def interpret(soln, max_tension): - phi = soln[:, 0] - phi1 = soln[:, 1] - v_tan = v_cm_co(phi) * phi1 * r(phi) / 1000 #kilometers/second - tension = (m*r_sat(phi)/r(phi))*(phi1**2 * r(phi) - g(tau/4))/1000 #kN - cable_moment = l*((r(phi) - r_cm(phi))**3 - r_cm(phi)**3)/3 #kg m^2 - moment = m*(r(phi) - r_cm(phi))**2 + cable_moment - M(phi)*r_cm(phi)**2 - momentum = moment*phi1**2 #kg m^2/s^2 - - try: - maxindex = numpy.argmax(tension >= max_tension) - max_arg = numpy.argmax(v_tan[:maxindex]) - except ValueError: max_arg = numpy.argmax(v_tan) - - if phi[max_arg] > abs(phi[-1] - tau): - raise TimespaceTooSmall - - return phi, phi1, v_tan, tension, moment, momentum, max_arg - -def solver(func, init_vars, timespace=(0,600,10**5), max_tension=max_tension, - returns=['time', 'variables', 'slices', 'maxes']): - time0 = current_time() - for n in itertools.count(1): - time, soln = solve_fun(func, init_vars, timespace) - - try: variables, max_arg = interpret(soln, max_tension) - except TimespaceTooSmall: - timespace = map(lambda x: x*.2*1.5**n, timespace) - continue - break - - getslice = lambda x: numpy.ndarray.__getslice__(x, 0, max_arg+1) - getmax = lambda x: numpy.ndarray.__getitem__( x, max_arg ) - slices = map(getslice, (time,) + variables) - maxes = map(getmax, (time,) + variables) - - print ": took %s seconds to solve" % (current_time() - time0) - - out = [] - for item in returns: out.append(vars()[item]) - return tuple(out) - -def plot(solved, graphs=['all'], i=7): - time, variables, slices, maxes = solved - phi, phi1, v_tan, tension, moment, momentum = variables - (phi_max, phi1_max, v_tan_max, - tension_max, moment_max, momentum_max) = maxes - - 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) - - i -= 1 - if test('moment'): - pyplot.figure(i) - pyplot.plot(phi / tau, moment) - pyplot.grid(True) - pyplot.ylabel("Moment of Inertia (kg m^2)") - pyplot.xlabel("phi (cycles)") - - i -= 1 - if test('momentum'): - pyplot.figure(i) - pyplot.plot(phi / tau, momentum) - pyplot.grid(True) - pyplot.ylabel("Angular Momentum (kg m^2 s^-2)") - pyplot.xlabel("phi (cycles)") - - #import pdb; pdb.set_trace() - -def opt_fun(guesses, f, func, init_vars): - """Solve, plot, and compare the diff eq.""" - global r0, r_roll, length, m, Mass - r0, r_roll, length, m, Mass = guesses - - time, soln = solve_fun(func, init_vars) - - phi, phi1, v_tan, tension, moment, momentum, max_arg = \ - interpret(soln) - phi, phi1, v_tan, tension, moment, momentum, time = \ - map(lambda x: numpy.ndarray.__getslice__(x, 0, max_arg+1), - [phi, phi1, v_tan, tension, moment, momentum, time]) - - time, (phi, phi1, v_tan, tension, moment, momentum, max_arg) = \ - solver(func, init_vars, returns=['time', 'slices']) - - # Parameter values. - f.write("%s, " % r_roll) - f.write("%s, " % length) - f.write("%s, " % r0) - f.write("%s, " % m) - f.write("%s, " % Mass) - - # Values at the max velocity point. - f.write("%s, " % phi[-1]) - f.write("%s, " % phi1[-1]) - f.write("%s, " % tension[-1]) - f.write("%s\n" % v_tan[-1]) - - # If you want strict limits, you can change the return to: - #return 1/v_tan + 1000*(Mass > Mass_limit) + 1000*(phi1.max() > phi1_limit) - return 2/v_tan[-1] + Mass/600 + phi1[-1] / 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, r_roll, length of cable, m, Mass, phi, omega, Tension, Max Vel.\n') - guess = (r0, r_roll, length, 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(): - solved = solver(throw, inits) - plot(solved, ['all', '!gravity']) - pyplot.show() - -if len(sys.argv) == 2: - func_name = sys.argv[1] - if func_name in dir(): - exec(func_name + '()') - -- cgit v1.2.3-70-g09d2