path: root/throw.py

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

### trig approximations
cos = taylor_cos
sin = taylor_sin

### radii:
global r0
r0 = 1. #meters -- the length of cable in 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*(1 - phi/(2*phi_max))*phi # cable let out
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

### cable data:
global dia, l, max_tension
#diameter in m: [l in kg/m, breaking point in kN]
cables = {.006: [.023,   41.2], .008: [.039,   65.7],
          .009: [.042,   82.4], .010: [.059,  105.9],
          .011: [.063,  133.4], .012: [.095,  161.9],
          .014: [.128,  215.8], .016: [.160,  269.8],
          .018: [.208,  343.3], .020: [.255,  407.1],
          .022: [.305,  490.5], .024: [.358,  569.0],
          .026: [.410,  647.4], .028: [.465,  725.9],
          .030: [.520,  799.5], .032: [.570,  868.1],
          .034: [.625,  941.2], .036: [.680, 1020.2],
          .038: [.740, 1098.7], .040: [.840, 1245.8],
          .042: [.930, 1373.4], .044: [1.02, 1491.1],
          .046: [1.11, 1618.6], .048: [1.21, 1755.9],
          .050: [1.31, 1893.3], .052: [1.41, 2020.8],
          .056: [1.63, 2315.1], .060: [1.75, 2472.0],
          .064: [2.00, 2766.3], .068: [2.26, 3099.9],
          .072: [2.54, 3413.8], .080: [3.13, 4139.7],
          .088: [3.79, 4934.3], .096: [4.51, 5768.1]}
dia = .026
cable_choices = numpy.array(cables.keys())
nearest = lambda dia: cable_choices[numpy.abs(dia - cable_choices).argmin()]
l, max_tension = cables[nearest(dia)]

### masses:
global m, Mass
m = 39.7 #kilograms -- mass of the satellite cage (overestimate)
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)

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: cos(numpy.arctan(tan(phi))) # cos(theta)
    tension = (m/cos_t(phi))*(phi1**2 * r(phi) - g(tau/4))/1000 #kiloNewtons
    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)
    
    return phi, phi1, v_tan, cos_t, tension, moment, momentum, max_arg

def plot(time, soln, graphs=['all'], i=7):
    phi, phi1, v_tan, cos_t, tension, moment, momentum, max_arg = \
        interpret(soln)
    phi_max, phi1_max, v_tan_max, tension_max, moment_max, momentum_max = \
        map(lambda x: numpy.ndarray.__getitem__(x, max_arg),
            [phi, phi1, v_tan, tension, moment, momentum])
    
    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, m, Mass, l, dia, max_tension
    r0, m, Mass, dia = guesses
    l, max_tension = cables[nearest(dia)]
    
    time, soln = solve_fun(func, init_vars)
    
    phi, phi1, v_tan, cos_t, tension, moment, momentum, max_arg = \
        interpret(soln)
    phi_max, phi1_max, v_tan_max, tension_max, moment_max, momentum_max = \
        map(lambda x: numpy.ndarray.__getitem__(x, max_arg),
            [phi, phi1, v_tan, tension, moment, momentum])
    
    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)

    return 2/v_tan_max + Mass/600 + phi1.max() / tau + cable*50

v0 = numpy.sqrt(g(0) * r(0))
w0 = v0/r0
inits = [0, w0]

#outfile = 'C:\\Users\\dave\\Desktop\\out'
outfile = '../out'
def optimize():
    fil = file(outfile, 'w')
    fil.write('r0, sat. mass, Total Mass, cable diameter, '
            + '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 + '()')