Contradiction in Euler Angle Rotation Matrix Calculation

------------ The Setup --------------------------

Assume that I have an object in 3D space. Its orientation is defined by the transform below:

| 1 0 0 |
| 0 1 0 |
| 0 0 1 |

As can be seen, the object’s local axes align with the world’s axes, so the object’s transform is an identity matrix.

-----------The First Rotation --------------------------------

Now, I will first rotate the object around the X axis by multiply the object’s transform with the rotation matrix below:

| 1 0 0 |
| 0 Cos(ψ) -Sin(ψ) |
| 0 Sin(ψ) Cos(ψ) |

As can be seen, the matrix presented above is a rotation matrix around the X axis, and Ψ is the amount of rotation around the X axis. Here, I ask which X axis will the first rotation rotates the object around: the object’s local X axis or the world X axis? The answer is that since the object’s local axes align with the world’s axes, so the first rotation is a rotation around both the object’s local X axis and the world X axis.

-----------The Second Rotation -----------------------------------

Now, I will do another rotation operation by multiplying the object’s transform with another rotation matrix presented below:

| Cos(θ) 0 Sin(θ) |
| 0 1 0 |
| -Sin(θ) 0 Cos(θ) |

As can be seen, the matrix presented above is a rotation matrix around the Y axis, and θ is the amount of rotation around the Y axis. Here, I ask again which Y axis will this second rotation rotates the object around: the object’s local Y axis or the world Y axis? Keep in mind that after the first rotation, the object’s local Y axis and the world’s Y axis are no longer align, so which Y axis is the second rotation be around?

According to many textbooks, papers and Youtube math tutorials, the answer is that the second rotation will be around the object’s local Y axis. However, my calculation below conflicts with this answer.

M = Ry(θ)Ry(ψ)O where θ = 18, ψ = 22.2, O is the starting transform of the object, M is the final result

| Cos(18) 0 Sin(18) | | 1 0 0 | | 1 0 0 |
M = | 0 1 0 | | 0 Cos(22.2) -Sin(22.2) | | 0 1 0 |
| -Sin(18) 0 Cos(18) | | 0 Sin(22.2) Cos(22.2) | | 0 0 1 |

| 0.951 0.116759 0.2861 |
M = | 0 0.92587 -0.37784 |
| -0.3090 0.359347 0.880555 |

The fact that you see zero in M21(second row first column) suggests that the second rotation rotates the object around the world’s Y axis, NOT the object’s local Y axis. As can be seen, my calculation contradicts the math in many textbook and papers.

Using Godot, I had verified the fact that if the second rotation is around the object’s local Y axis, the correct result should be as shown below:

| 0.951057 0 0.309017 |
M = | 0.116759 0.925871 -0.359348 |
| -0.28611 0.377841 0.880555 |

Please explain why this contradiction happens.