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()