Initial State

1 files changed, 326 insertions, 0 deletions

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