ang_vel_in() could give a simpler result for body fixed angular velocities

[moorepants]

If you set up a body fixed rotation and then call .ang_vel_in() I would expect the simplest form of the angular velocity to be returned, because we know apriori what the angular velocity is for any body fixed rotation.

import sympy as sm
import sympy.physics.mechanics as me

psi, theta, phi = me.dynamicsymbols('psi, theta, varphi')

A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')
B.orient_body_fixed(A, (psi, theta, phi), 'ZXZ')
B.ang_vel_in(A)

gives:

⎛                         d                              d            ⎞       ⎛                         d                              d
⎜sin(θ(t))⋅sin(varphi(t))⋅──(ψ(t))        cos(varphi(t))⋅──(θ(t))     ⎟       ⎜sin(θ(t))⋅cos(varphi(t))⋅──(ψ(t))        sin(varphi(t))⋅──(
⎜                         dt                             dt           ⎟       ⎜                         dt                             dt
⎜───────────────────────────────── + ─────────────────────────────────⎟ b_x + ⎜───────────────────────────────── - ───────────────────────
⎜   2                 2                 2                 2           ⎟       ⎜   2                 2                 2                 2
⎝sin (varphi(t)) + cos (varphi(t))   sin (varphi(t)) + cos (varphi(t))⎠       ⎝sin (varphi(t)) + cos (varphi(t))   sin (varphi(t)) + cos (

          ⎞
θ(t))     ⎟
          ⎟       ⎛          d          d            ⎞
──────────⎟ b_y + ⎜cos(θ(t))⋅──(ψ(t)) + ──(varphi(t))⎟ b_z
          ⎟       ⎝          dt         dt           ⎠
varphi(t))⎠

We know the answer to this and I though that the code was supposed to return the simpler form for a lookup table. The simple form is:

>>> B.ang_vel_in(A).simplify()
⎛                         d                         d       ⎞       ⎛                         d                         d       ⎞       ⎛
⎜sin(θ(t))⋅sin(varphi(t))⋅──(ψ(t)) + cos(varphi(t))⋅──(θ(t))⎟ b_x + ⎜sin(θ(t))⋅cos(varphi(t))⋅──(ψ(t)) - sin(varphi(t))⋅──(θ(t))⎟ b_y + ⎜c
⎝                         dt                        dt      ⎠       ⎝                         dt                        dt      ⎠       ⎝

         d          d            ⎞
os(θ(t))⋅──(ψ(t)) + ──(varphi(t))⎟ b_z
         dt         dt           ⎠

The manual calculation is:

mnx = me.dot(B.y.express(A).dt(A), B.z)
mny = me.dot(B.z.express(A).dt(A), B.x)
mnz = me.dot(B.x.express(A).dt(A), B.y)
A_w_B = mnx*B.x + mny*B.y + mnz*B.z
A_w_B.simplify()

We have this function:
https://github.com/sympy/sympy/blob/master/sympy/physics/vector/functions.py#L227
which has this big lookup table but I don't think it is used.

[moorepants]

I guess the issue is that direction cosine matrices have no knowledge of how they were created, but when B.orient_body_fixed(A) is called, the simplest form of B.ang_vel_in(A) should be computed and stored.

[soma2000-lang]

@moorepants I am working on this issue

[moorepants]

I think that when orien_* is called we should also set the angular velocity to be a simple form expressed in the reference frame that is being rotated. I think what happens now is that we wait until calling ang_vel_in to calculate the angular velocity and then it has to use a generic formula. It is also possible that the generic formula can get simpler results too.

[moorepants]

@soma2000-lang if you want to work on this, the best way to signal that is to open a draft pull request where we can observe your work on it. That lays "claim" a bit more than just a comment in this issue. I can also give feedback there as you work on it.

[moorepants]

This actually looks like a regression. In SymPy 1.0, I get the preferred answer:

In [1]: import sympy as sm

In [2]: sm.__version__
Out[2]: '1.0'

In [3]: %paste
import sympy as sm
import sympy.physics.mechanics as me

psi, theta, phi = me.dynamicsymbols('psi, theta, varphi')

## -- End pasted text --

In [4]: %paste
A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')

## -- End pasted text --

In [5]: me.init_vprinting()

In [6]: B.orient(A, 'Body', (psi, theta, phi), 'ZXZ')

In [7]: B.ang_vel_in(A)
Out[7]: (sin(θ)⋅sin(varphi)⋅ψ̇ + cos(varphi)⋅θ̇) b_x + (sin(θ)⋅cos(varphi)⋅ψ̇ - sin(varphi)⋅θ̇) b_y + (cos(θ)⋅ψ̇ + varṗhi) b_z

[moorepants]

Git bisect tells me that this is the commit that changed the result:

f0b7996

[moorepants]

This is the PR in which the commit appears: #18844

[moorepants]

This is the code that generates the angular velocity and it does use kinematic_equations. Interesting that there is now some a polyerrors reference there. I don't know what that would have to do with this module.

sympy/sympy/physics/vector/frame.py

Lines 882 to 898 in 1de2e07

	try:
	from sympy.polys.polyerrors import CoercionFailed
	from sympy.physics.vector.functions import kinematic_equations
	q1, q2, q3 = amounts
	u1, u2, u3 = symbols('u1, u2, u3', cls=Dummy)
	templist = kinematic_equations([u1, u2, u3], [q1, q2, q3],
	'body', rot_order)
	templist = [expand(i) for i in templist]
	td = solve(templist, [u1, u2, u3])
	u1 = expand(td[u1])
	u2 = expand(td[u2])
	u3 = expand(td[u3])
	wvec = u1 * self.x + u2 * self.y + u3 * self.z
	except (CoercionFailed, AssertionError):
	wvec = self._w_diff_dcm(parent)
	self._ang_vel_dict.update({parent: wvec})
	parent._ang_vel_dict.update({self: -wvec})

[moorepants]

The line 890 td = solve(templist, [u1, u2, u3]) is where this new form comes from. So solve() is doing something worse for these input equations.

[moorepants]

Ok I opened #23140 for the deep issue, but we can also improve orient_body_fixed(). Right now the method extracts the q' = expr from kinematic_equations() and solves the set of three equations for the u's. We could simply hard code the simplest form of the u = expr also and not rely on solving the linear system. Also, it is a simple 3x3 linear system and we could write our own solve() function specifically for a 3x3 that returns the ideal result, swapping this for SymPy`'s solve.

[moorepants]

Also, I think this issue plays a role here: #23075

[moorepants]

A super quick fix would be to add trigsimp right after solve in orient_body_fixed(). That would restore the old behavior. In #23140 Oscar explained how simplify in linear solve was removed in SymPy 1.7. But it would be better if we hard coded the answers and returned those.