;+
; NAME:
; QTPOW
;
; 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:
; Raise quaternion Q to the "power" POW
;
; MAJOR TOPICS:
; Geometry
;
; CALLING SEQUENCE:
; QNEW = QTPOW(Q, POW)
;
; DESCRIPTION:
;
; The function QTPOW raises a quaterion Q to the power P. The
; operation
;
; QNEW = QTPOW(Q, POW)
;
; is equivalent to
;
; QNEW = QTEXP( POW * QTLOG(Q))
;
; which is the same as the definition of raising a real number to
; any power (however, QTPOW is faster than using QTLOG and QTEXP).
;
; For integer values of POW, this form of exponentiation is also
; directly equivalent to the multiplication of that many Q's
; together.
;
; Geometrically, raising Q to any power between 0 and 1 realizes a
; rotation that smoothly interpolates between the identity
; quaternion and Q. Thus, QTPOW is useful for interpolation of
; quaternions or SLERPing (spherical linear interpolation).
;
; When raising more than one quaternion to a power at a time, the
; number of quaternions and powers must be equal.
;
; Conventions for storing quaternions vary in the literature and from
; library to library. This library uses the convention that the
; first three components of each quaternion are the 3-vector axis of
; rotation, and the 4th component is the rotation angle. Expressed
; in formulae, a single quaternion is given by:
;
; Q(0:2) = [VX, VY, VZ]*SIN(PHI/2)
; Q(3) = COS(PHI/2)
;
; where PHI is the rotation angle, and VAXIS = [VX, VY, VZ] is the
; rotation eigen axis expressed as a unit vector. This library
; accepts quaternions of both signs, but by preference returns
; quaternions with a positive 4th component.
;
;
; INPUTS:
;
; Q - array of one or more unit quaternions. For a single
; quaternion, Q should be a 4-vector. For N quaternions, Q
; should be a 4xN array.
;
; POW - array of N powers, where N is the number of quaternions.
;
;
; RETURNS:
;
; The resulting exponentiated unit quaternions. For a single
; inputs, returns a 4-vector. For N input quaternions, returns N
; quaternions as a 4xN array.
;
;
; KEYWORD PARAMETERS:
;
; NONE
;
; EXAMPLE:
;
; ;; Form a rotation quaternion of 45 degrees about the X axis
; Q = qtcompose([1,0,0], !dpi/4)
;
; ;; Make an array of 1001 values smoothly varying from 0 to 1
; P = dindgen(1001)/1000d
;
; ;; Perform spherical linear interpolation
; QNEW = QTERP(Q, P)
;
;
; SEE ALSO
; QTANG, QTAXIS, QTCOMPOSE, QTERP, QTEXP, QTFIND, QTINV, QTLOG,
; QTMAT, QTMULT, QTPOW, QTVROT
;
; MODIFICATION HISTORY:
; Written, July 2001, CM
; Documented, Dec 2001, CM
; Usage message, error checking, 15 Mar 2002, CM
;
; $Id: qtpow.pro,v 1.5 2002/05/09 23:03:27 craigm Exp $
;
;-
; Copyright (C) 2001, 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.
;-
function qtpow, q, pow
if n_params() EQ 0 then begin
info = 1
USAGE:
message, 'USAGE:', /info
message, 'QNEW = QTPOW(Q, POW)', info=info
return, 0
endif
nq = n_elements(q)/4
np = n_elements(pow)
if nq LT 1 OR np LT 1 then goto, USAGE
v = q(0:2,*)
sinth = sqrt(total(v^2,1))
th = atan(sinth, q(3,*))
rat = th*0
wh = where(sinth NE 0, ct)
if ct GT 0 then rat(wh) = (sin(pow*th))(wh)/sinth(wh)
q1 = q
q1(3,*) = cos(th*pow)
q1(0:2,*) = rebin(reform(rat,1,nq),3,nq)*temporary(v)
return, q1
end