import sympy as sm
import sympy.physics.mechanics as me
me.init_vprinting()

theta, beta, alpha, l = me.dynamicsymbols('theta, beta, alpha, l')
A, B, C, D = sm.symbols('A, B, C, D', cls=me.ReferenceFrame)

B.orient_axis(A, theta, A.z)
C.orient_axis(B, beta, -B.x)
D.orient_axis(C, alpha, C.z)

r = l *D.x
print(r)

print(r.express(C))

print(r.express(B))

print(r.express(A))

print(r.dot(C.x).diff(alpha))

print(r.dot(C.y).diff(alpha))

r_alpha_C = r.dot(C.x).diff(alpha) * C.x + r.dot(C.y).diff(alpha) * C.y
print(r_alpha_C)

print(r.diff(alpha, C).express(C))

t = me.dynamicsymbols._t
print(theta.diff(t))