2.5.B.sagews

from sympy.physics.units import newton, m
from sympy import *
from sympy import N as Num

from mpmath import degrees, radians

print("------ User input ----------------------")

symbolic = True
symbolic = False

if symbolic:
    # quantity = symbol:
    M = var("M")
    G = var("G")
    l = var("l")
    r1, r2 = var("r1, r2")
else:
    # quantity =  factor times unit:
    M = 10 *newton*m
    G = 79*10**9 *newton/m/m
    l = 1 *m
    r1 = S(5)/1000 *m
    r2 = S(10)/1000 *m

print("------ Radius --------------------------")
x = var("x")

# Radius as a function of x:
rx = r2 - (r2-r1)*x/l

# polar second moment of area:
Ip = pi/2*rx**4

print("------ Integration ---------------------")
# integrand:
i = 1/Ip.simplify()
I = integrate(i, (x, 0, l))
pprint(["I", I.simplify()])

# check manual integration:
# derivative of antiderivative of
# integrand must be equal to integrand:
term = -l/(r2-r1)
anti = -2/pi/3*rx**(-3)*term
i_check = diff(anti, x).simplify()
print("Should be \"True\":", i == i_check)

theta_l = M/G*I.simplify()
pprint(["theta(l)", theta_l])
pprint(["theta(l) / deg", Num(degrees(theta_l),2)])

------ User input ---------------------- ------ Radius -------------------------- ------ Integration --------------------- 2800000000 [I, ----------] 3 3*pi*m ('Should be "True":', True) 28 [theta(l), ------] 237*pi [theta(l) / deg, 2.2]