3.3.H.sagews

from sympy import *
from sympy import N as Num

# SI units:
m  = Symbol("m",  positive=True)
s  = Symbol("s",  positive=True)
kg = Symbol("kg", positive=True)


print("\n--- User input -----------------------")

symbolic = True
symbolic = False

if symbolic:
    # quantity = symbol:
    g, M, r, R = var("g, M, r, R")
else:
    # factors:
    g_val = 10
    M_val = S(150)/ 1000
    r_val = S(10) / 1000
    R_val = S(40) / 1000
    #
    # quantity =  factor times unit:
    g = g_val * m/s/s
    M = M_val * kg
    (r, R) = (r_val*m, R_val*m)

print("\n--- b: -------------------------------")
Theta = M*R*R/2
Theta_A = Theta + r*r*M

# x = x''
# p = phi''
S = var("S")
x = var("x")
p = var("p")

# Balance of momentum:
eq1 = Eq(M*x, M*g - S)

# Balance of moment of momentum:
eq2 = Eq(Theta*p, r*S)     # both equations
eq2 = Eq(Theta_A*p, r*M*g) # are correct

# Kin. constraint:
eq3 = Eq(x, r*p)

eqs = [eq1, eq2, eq3]
unknowns = [S, x, p]
sol = solve(eqs, unknowns)

for s in sol:
    pprint("\n")
    pprint(s)
    pprint(sol[s])

x=sol[x]

print("\n--- c: -------------------------------")
t1 = sqrt(2*m/x)
pprint(["t_1", Num(t1,3)])

print("\n--- d: -------------------------------")
v1 = g*t1/9
pprint(["v_1", Num(v1,3)])

print("\n--- e: -------------------------------")
w1 = v1/r
vB = w1 * (r+R)
pprint(["v_B", Num(vB,3)])

--- User input ----------------------- --- b: ------------------------------- S 4*kg*m ------ 2 3*s x 10*m ---- 2 9*s p 1000 ---- 2 9*s --- c: ------------------------------- [t_1, 1.34*s] --- d: ------------------------------- 1.49*m [v_1, ------] s --- e: ------------------------------- 7.45*m [v_B, ------] s