Angular Momentum graphs

1 files changed, 37 insertions, 14 deletions

diff --git a/throw.py b/throw.py
index deaea5b..3fcfdd8 100644
--- a/throw.py
+++ b/throw.py
@@ -66,13 +66,13 @@ sin = taylor_sin
 
 ### radii:
 global r0
-r0 = 1. #meters -- the radius of the first full-circle rotation
+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*phi # the radius of the satellite
-r = lambda phi: r0 + r_winch*(phi - phi**2/(2*phi_max))
+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
@@ -93,10 +93,12 @@ quad_co = kappa * density * Area
 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
-max_tension = 269.8 #kN, given by choice of cable above.
+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)
-v_cm_co = lambda phi: (.5 * l*r(phi) + M(phi)) / (m + l*r(phi) + M(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):
@@ -124,21 +126,25 @@ 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)
+    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, max_arg
+    return phi, phi1, v_tan, cos_t, tension, moment, momentum, 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 = \
+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])
+            [phi, phi1, v_tan, tension, moment, momentum])
     
     def test(*args):
         b = 0
@@ -197,6 +203,22 @@ def plot(time, soln, graphs=['all'], i=5):
         
         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):
@@ -206,10 +228,11 @@ def opt_fun(guesses, f, func, init_vars):
     
     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 = \
+    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])
+            [phi, phi1, v_tan, tension, moment, momentum])
     
     f.write("%s, " % r0)
     f.write("%s, " % m)