Transform torque from Euler angles to infinitesimal Cartesian rotations For a certain pair of rigid bodies, I have the gradient of energy in terms of Euler angles. I want to transform this gradient to the gradient of energy in terms of rotations about the $x, y, z$ axes (the very useful gradient known as torque). 
Performing the analogous transformation between, say, spherical coordinates and Cartesian is easy because both form a full basis. One simply expresses the spherical coordinates in terms of Cartesian and finds the partial derivatives for each $r, \theta, \phi$ with respect to $x, y, z$ (the Jacobian). Then matrix multiplication gives the answer. 
But the Cartesian rotations are not a basis for full rotations! They are only a basis for infinitesimal rotations. Therefore it is not possible to take this simple chain-rule approach, because one cannot represent an arbitrary Euler angle in terms of fixed-axis $x, y, z$ rotations. 
I strongly suspect that this Jacobian exists, especially since this Phys.SE post seems to have found it (without stating the Euler convention or any derivation or sources, making it useless). How can I find this Jacobian?
 A: I think this is a very messy problem no matter how you approach it. Let me show you how to at least eliminate the calculus so I can leave you a linear algebraic problem to solve.
I don't know which Euler angle convention you're using. I'll guess:
\begin{align}
R(\phi,\theta,\psi) &=& R_z(\phi)R_y(\theta)R_z(\psi)
\end{align}
Let's rewrite this using quaternions:
\begin{align}R(\phi,\theta,\psi) &=& \exp(k\phi/2)\exp(j\theta/2)\exp(k\psi/2)
\end{align}
We're interested in how this rotation changes when we multiply by rotations through infinitesimal $x$ round the $x$-axis, $y$ round the $y$-axis and $z$ around the $z$-axis. (Note that order doesn't matter because we're asking about infinitesimals and computing first derivatives.)
\begin{align}
R+\delta R &=& \exp(ix/2)\exp(jy/2)\exp(kz/2)\exp(k\phi/2)\exp(j\theta/2)\exp(k\psi/2) \\
\end{align}
Using the fact that $x$, $y$ and $z$ are infinitesimal we get:
\begin{align}
R+\delta R &=& (1+ix/2+jy/2+kz/2)\exp(k\phi/2)\exp(j\theta/2)\exp(k\psi/2) \\
\delta R &=& (ix/2+jy/2+kz/2)\exp(k\phi/2)\exp(j\theta/2)\exp(k\psi/2)
\end{align}
But we also have:
\begin{align}
R+\delta R &=& \exp(k(\phi+\delta\phi)/2)\exp(j(\theta+\delta\theta)/2)\exp(k(\psi+\delta\psi)/2)\\
\delta R &=& (\delta\phi/2)k\exp(k\phi/2)\exp(j\theta/2)\exp(k\psi/2)\\
&&+(\delta\theta/2)\exp(k\phi/2)j\exp(j\theta/2)\exp(k\psi/2)\\
&&+(\delta\psi/2)\exp(k\phi/2)\exp(j\theta/2)k\exp(k\psi/2)
\end{align}
We want to use the fact that $\delta\psi={\partial\psi\over\partial x}dx+{\partial\psi\over\partial y}dy+{\partial\psi\over\partial z}dz$ etc. to get all of the partial derivatives. You can do this by simplifying the expressions above using basic quaternion algebra and $\exp(\theta n)=cos\theta+n\sin\theta$ for any unit quaternion $n$ and then equating coefficients of $1$, $i$, $j$ and $k$. Start by multiplying on the right by $R^{-1}$ because $\delta R R^{-1}$ has zero real part and you're then left with three linear equations in three variables. Shouldn't be too bad in a good computer algebra package. I'll let you finish it. By linearity you can drop all the divisions by two. That'll save some ink.
This sidesteps the horrible problem of converting to Euler coordinates and back and differentiating the resulting mess.
