OV-GSO-DC / dustem-wrapper_idl

Blame view

src/idl_extern/CMTotal_for_Dustemwrap/qtexp.pro 3.31 KB
Edit Raw Normal View History
517b8f98   Annie Hughes   first commit 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
 
;+
; NAME:
;   QTEXP
;
; 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:
;   Compute "exponentiation" of a non-unit quaternion
;
; MAJOR TOPICS:
;   Geometry
;
; CALLING SEQUENCE:
;   Q = QTEXP(QLOG)
;
; DESCRIPTION:
;
;   The function QTEXP computes the "exponentiation" of a quaternion.
;
;   The exponential is only defined for a non-unit quaternion with a
;   *zero* rotation angle.  Specifically, the expression
;
;      QTEXP([VAXIS * PHI/2, 0])    
;
;   becomes
; 
;      [VAXIS*SIN(PHI/2), COS(PHI/2)]
;
;   where VAXIS is the unit vector rotation axis and PHI is the
;   rotation angle.  Note that since VAXIS is a unit vector, the
;   product VAXIS*PHI can have an arbitrary direction and magnitude.
;
;   Typically the input to QTEXP is found by taking the logarithm of a
;   unit quaternion using QTLOG, and the identity QTEXP(QTLOG(Q)) is
;   the same as Q.
;   
;  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:
;
;   QLOG - a non-unit quaternion of the form [VX, VY, VZ, 0]; or,
;          N quaternions of the same form, as a 4xN array.
;
;
; RETURNS:
;
;   The exponentiated unit quaternion(s).  For a single input
;   quaternion, returns a 4-vector; for N input quaternions, returns a
;   4xN array.
;
;
; KEYWORD PARAMETERS:
;
;   NONE
;
; EXAMPLE:
;
;   IDL> print, qtlog(qtcompose([0d,1,0], !dpi/4))
;         0.0000000      0.39269908       0.0000000       0.0000000
;
;   Prints the logarithm of the quaternion composed of a rotation of
;   !dpi/4 radians around the axis [0,1,0]
;
;
; SEE ALSO
;   QTANG, QTAXIS, QTCOMPOSE, QTERP, QTEXP, QTFIND, QTINV, QTLOG,
;   QTMAT, QTMULT, QTPOW, QTVROT
;
; MODIFICATION HISTORY:
;   Written, July 2001, CM
;   Documented, Dec 2001, CM
;   Documentation corrected, 27 Jan 2002, CM
;   Usage message, error checking, 15 Mar 2002, CM
;
;  $Id: qtexp.pro,v 1.7 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 qtexp, q

  if n_params() EQ 0 then begin
      info = 1
      USAGE:
      message, 'USAGE:', /info
      message, 'QNEW = QTEXP(Q)', info=info
      return, 0
  endif
  nq = n_elements(q)/4
  if nq LT 1 then goto, USAGE

  v = q(0:2,*)
  th = sqrt(total(v^2,1))
  wh = where(th NE 0, ct)

  if ct GT 0 then v(*,wh) = v(*,wh)/rebin(reform(th(wh),1,ct),3,ct)

  q1 = q

  q1(3,*) = cos(th)
  q1(0:2,*) = rebin(reform([sin(th)],1,nq),3,nq)*v

  return, q1
end