;+ ; NAME: ; QTEULER ; ; AUTHOR: ; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770 ; craigm@lheamail.gsfc.nasa.gov ; UPDATED VERSIONs can be found on my WEB PAGE: ; http://cow.physics.wisc.edu/~craigm/idl/idl.html ; ; PURPOSE: ; Compose a series of euler-type rotations into a single quaternion ; ; MAJOR TOPICS: ; Geometry ; ; CALLING SEQUENCE: ; Q = QTEULER(AXES, ANG0, ANG1, ... ) ; ; DESCRIPTION: ; ; The function QTEULER composes a series of Euler-type rotations into ; a single set of quaternion representations. ; ; The user specifies a set of axes, and the angles to rotation about ; those axes, and QTEULER constructs the corresponding quaternion. ; ; There must be a one-to-one correspondence between the elements of ; AXES and the number of rotations. AXES specifies the rotation axes ; as an string, which must be one of 'X', 'Y', or 'Z'. Other axes ; are invalid. For example, the following call: ; ; QTEULER(['X','Z'], THETA, PHI) ; ; will rotate first about the *Z* axis by the angle PHI, and then ; around the *resulting X* axis by angle THETA. ; ; Several things are worth noting here. First, rotations are applied ; first from the right, not the left. This conforms to the usual ; matrix notation for applying rotations to a vector on the right ; hand side. For example, in matrix notation, ; ; XNEW = A3 A2 A1 XOLD ; ; applies first A1, then A2 and finally A3 to the XOLD vector, ; resulting in the new vector XNEW. The same semantics apply here. ; ; A second thing to bear in mind is that the axes themselves change ; during the rotations. Thus, the coordinates specified in AXES ; should be considered attached to the "body" and not the inertial ; frame. ; ; ; INPUTS: ; ; AXES - a string array, specifies the rotation axes. Rotations are ; applied last element first. Each element of AXES must be ; one of 'X', 'Y' or 'Z'. ; ; ANG0, ..., ANGi - the successive rotation angles. [radians] ; Angle ANGi corresponds to axis AXES(i). ; ; If ANGi is a scalar, then it will be promoted to a vector ; the same size as the other rotation angles being performed. ; Otherwise, if the angles ANGi are vectors, then they must ; all be of the same size. ; ; RETURNS: ; ; The resulting quaternion (or, if ANGi are vectors, array of ; quaternions), which represent the requested rotations. ; ; ; KEYWORD PARAMETERS: ; ; NONE ; ; EXAMPLE: ; ; ;; Precession Nutation ; qtot = qteuler(['z','y','z', 'x','z','x' ], $ ; -zeta, +theta, -z, +eps0, -dpsi, -eps) ; ; Applies a series of rotations to correct for earth nutation and ; precession. The order of rotations on a vector would be ; X-Z-X-Z-Y-Z (i.e., the reverse order printed). ; ; SEE ALSO ; QTANG, QTAXIS, QTCOMPOSE, QTERP, QTEXP, QTFIND, QTINV, QTLOG, ; QTMAT, QTMULT, QTPOW, QTVROT ; ; MODIFICATION HISTORY: ; Written, 27 Jan 2002, CM ; More error checking, 03 Mar 2002, CM ; ; $Id: qteuler.pro,v 1.5 2012/09/27 23:53:08 cmarkwar Exp $ ; ;- ; Copyright (C) 2002, Craig Markwardt ; This software is provided as is without any warranty whatsoever. ; Permission to use, copy, modify, and distribute modified or ; unmodified copies is granted, provided this copyright and disclaimer ; are included unchanged. ;- ;; Extract axis ei and angle angi pro qteuler_extract, ax, i, ei, angi, $ ang0, ang1, ang2, ang3, ang4, $ ang5, ang6, ang7, ang8, ang9, $ status=status, errmsg=errmsg status = 0 zero = ang0(0)*0 ex = [1,zero,zero] & ey = [zero,1,zero] & ez = [zero,zero,1] ei = [0D, 0D, 0D] if execute('ei = e'+ax(i)+' & angi = ang'+strtrim(i,2)) NE 1 then begin stop errmsg = 'Invalid axis specification' return endif status = 1 return end function qteuler, axes, block=block, $ ang0, ang1, ang2, ang3, ang4, ang5, ang6, ang7, ang8, ang9, $ ang10, ang11, ang12, ang13, ang14, ang15 if n_params() EQ 0 then begin info = 1 USAGE_ERR: message, 'USAGE: Q = QTEULER(AXES, ANG0, ...)', /info message, ' AXES = ["X",...] ("X" or "Y" or "Z")', /info message, ' ANGn = rotation angle (radians)', info=info return, 0 endif if n_elements(axes) LT 1 OR n_elements(ang0) LT 1 then $ goto, USAGE_ERR nang = n_params()-1 ;; Check to be sure each axis label is 'X' 'Y' or 'Z' ax = strupcase(strmid(strtrim(axes,2),0,1)) wh = where(ax NE 'X' AND ax NE 'Y' AND ax NE 'Z', ct) if ct GT 0 then begin errmsg = 'AXES must be one of "X", "Y" or "Z"' goto, BAD_AXIS endif if n_elements(ax) NE nang then begin errmsg = 'Number of AXES and rotations ANGi must agree' goto, BAD_AXIS endif qteuler_extract, ax, 0, ev, angv, status=status, errmsg=errmsg, $ ang0, ang1, ang2, ang3, ang4, ang5, ang6, ang7, ang8, ang9 if status EQ 0 then begin BAD_AXIS: message, 'ERROR: '+errmsg, /info goto, USAGE_ERR endif qq = qtcompose(ev, angv) for i = 1, nang-1 do begin qteuler_extract, ax, i, ev, angv, status=status, errmsg=errmsg, $ ang0, ang1, ang2, ang3, ang4, ang5, ang6, ang7, ang8, ang9 if status EQ 0 then goto, BAD_AXIS qq = qtmult(qq, qtcompose(ev, angv)) endfor return, qq end