;+ ;Function: qdecompose,q ; ;Purpose: decompose quaternions into axes and angeles ; ;Inputs: q: a 4 element quaternion or an Nx4 element array of quaternions ; ;Returns: a 4 element array with a[0] = angle, and a[1:3] = axis, or ;an Nx4 element array or -1L on failure ; ;;Notes: Implementation largely copied from the euve c library for ;quaternions ;Represention has q[0] = scalar component ; q[1] = vector x ; q[2] = vector y ; q[3] = vector z ; ;The vector component of the quaternion can also be thought of as ;an eigenvalue of the rotation the quaterion performs ; ;As per the euve implementation, if q[0] is outside of the range of ;acos...[-1,1] the value of the quaternion will be turned into an ;identity quaternion...in other words clipped, this seems suspect, ;a better solution may be to wrap the value back into range using ;modular arithmatic, future modifiers of this routine should consider ;adding this. ; ; ;Written by: Patrick Cruce(pcruce@igpp.ucla.edu) ; ; $LastChangedBy: pcruce $ ; $LastChangedDate: 2007-11-11 17:12:08 -0800 (Sun, 11 Nov 2007) $ ; $LastChangedRevision: 2027 $ ; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/ssl_general/trunk/cotrans/cotrans.pro $ ;- function qdecompose,q EPSILON = 1.0e-20 ;Where sin(theta) is close enough to theta compile_opt idl2 ;this is to avoid mutating the input variable qi = q ;check to make sure input has the correct dimensions qi = qvalidate(qi,'q','qdecompose') if(size(qi,/n_dim) eq 0 && qi[0] eq -1) then return,qi qdims = size(qi,/dimensions) aout = make_array(qdims,/double) ;the following code will clip into range idx = where(qi[*,0] ge 1.0D) if(idx[0] ne -1) then begin aout[idx,0] = 0.0D aout[idx,1] = 1.0D aout[idx,2:3] = 0.0D endif idx = where(qi[*,0] le -1.0D) if(idx[0] ne -1) then begin aout[idx,0] = 2*!DPI aout[idx,1] = 1.0D aout[idx,2:3] = 0.0D endif idx = where(qi[*,0] gt -1.0D and qi[*,0] lt 1.0D) if(idx[0] ne -1) then begin theta2 = acos(qi[idx,0]) aout[idx,0] = 2 * theta2 idx2 = where(theta2 lt EPSILON) if(idx2[0] ne -1) then begin aout[idx[idx2],1] = 1.0D aout[idx[idx2],2:3] = 0.0D endif idx2 = where(theta2 ge EPSILON) if(idx2[0] ne -1) then begin aout[idx[idx2],1] = qi[idx[idx2],1] / sin(theta2[idx2]) aout[idx[idx2],2] = qi[idx[idx2],2] / sin(theta2[idx2]) aout[idx[idx2],3] = qi[idx[idx2],3] / sin(theta2[idx2]) endif endif return,reform(aout) end