5-var optimize, Started writing solver()

1 files changed, 59 insertions, 34 deletions

diff --git a/throw.py b/throw.py
index 2f8d36f..e36b151 100644
--- a/throw.py
+++ b/throw.py
@@ -65,13 +65,14 @@ cos = taylor_cos
 sin = taylor_sin
 
 ### radii:
-global r0
-r0 = 1. #meters -- the length of cable in the first full-circle rotation
+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_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_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:
@@ -100,7 +101,6 @@ 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
@@ -108,26 +108,27 @@ def throw(phi, time):
       # 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
+        ( (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, space=(0, 1874, 10**5)):
+def solve_fun(func, init_vars, timespace):
     time0 = current_time()
-    time = numpy.linspace(*space)
+    time = numpy.linspace(*timespace)
     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):
+
+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
-    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
+    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
@@ -137,14 +138,32 @@ def interpret(soln, max_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
+    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):
+    for n in itertools.count(1):
+        time, soln = solve_fun(func, init_vars, timespace)
+        
+        try: phi, phi1, v_tan, tension, moment, momentum, max_arg = \
+                interpret(soln, max_tension)
+        except TimespaceTooSmall:
+            timespace = map(lambda x: x*.2*1.5**n, timespace)
+            continue
+        break
+    
+    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])
 
 def plot(time, soln, graphs=['all'], i=7):
-    phi, phi1, v_tan, cos_t, tension, moment, momentum, max_arg = \
+    phi, phi1, v_tan, 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])
+    (phi_max, phi1_max, v_tan_max, tension_max, moment_max, momentum_max,
+     time_max) = map(lambda x: numpy.ndarray.__getitem__(x, max_arg),
+        [phi, phi1, v_tan, tension, moment, momentum, time])
     
     def test(*args):
         b = 0
@@ -223,28 +242,33 @@ def plot(time, soln, graphs=['all'], i=7):
 
 def opt_fun(guesses, f, func, init_vars):
     """Solve, plot, and compare the diff eq."""
-    global r0, m, Mass
-    r0, m, Mass = guesses
+    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, cos_t, tension, moment, momentum, max_arg = \
+    phi, phi1, v_tan, 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])
+    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])
     
+    # Parameter values.
+    f.write("%s, " % r_roll)
+    f.write("%s, " % length)
     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)
+    
+    # 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_max + Mass/600 + phi1.max() / tau
+    return 2/v_tan[-1] + Mass/600 + phi1[-1] / tau
 
 v0 = numpy.sqrt(g(0) * r(0))
 w0 = v0/r0
@@ -254,8 +278,9 @@ inits = [0, w0]
 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)
+    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()