1 /* 2 3 -Procedure raxisa_c ( Rotation axis of a matrix ) 4 5 -Abstract 6 7 Compute the axis of the rotation given by an input matrix 8 and the angle of the rotation about that axis. 9 10 -Disclaimer 11 12 THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE 13 CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. 14 GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE 15 ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE 16 PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" 17 TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY 18 WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A 19 PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC 20 SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE 21 SOFTWARE AND RELATED MATERIALS, HOWEVER USED. 22 23 IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA 24 BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT 25 LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, 26 INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, 27 REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE 28 REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. 29 30 RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF 31 THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY 32 CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE 33 ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. 34 35 -Required_Reading 36 37 ROTATION 38 39 -Keywords 40 41 ANGLE, MATRIX, ROTATION 42 43 */ 44 45 #include "SpiceUsr.h" 46 #include "SpiceZfc.h" 47 #undef raxisa_c 48 49 raxisa_c(ConstSpiceDouble matrix[3][3],SpiceDouble axis[3],SpiceDouble * angle)50 void raxisa_c ( ConstSpiceDouble matrix[3][3], 51 SpiceDouble axis [3], 52 SpiceDouble * angle ) 53 54 /* 55 56 -Brief_I/O 57 58 VARIABLE I/O DESCRIPTION 59 -------- --- -------------------------------------------------- 60 matrix I 3x3 rotation matrix in double precision. 61 axis O Axis of the rotation. 62 angle O Angle through which the rotation is performed. 63 64 -Detailed_Input 65 66 matrix is a 3x3 rotation matrix in double precision. 67 68 -Detailed_Output 69 70 axis is a unit vector pointing along the axis of the rotation. 71 In other words, `axis' is a unit eigenvector of the input 72 matrix, corresponding to the eigenvalue 1. If the input 73 matrix is the identity matrix, `axis' will be the vector 74 (0, 0, 1). If the input rotation is a rotation by pi 75 radians, both `axis' and -axis may be regarded as the 76 axis of the rotation. 77 78 angle is the angle between `v' and matrix*v for any non-zero 79 vector `v' orthogonal to `axis'. `angle' is given in 80 radians. The angle returned will be in the range from 0 81 to pi radians. 82 83 -Parameters 84 85 None. 86 87 -Exceptions 88 89 1) If the input matrix is not a rotation matrix (a fairly 90 loose tolerance is used to check this) a routine in the 91 call tree of this routine will signal an error indicating 92 the problem. 93 94 2) If the input matrix is the identity matrix, this routine 95 returns an angle of 0.0, and an axis of ( 0.0, 0.0, 1.0 ). 96 97 -Files 98 99 None. 100 101 -Particulars 102 103 Every rotation matrix has an axis `a' such any vector `v' 104 parallel to that axis satisfies the equation 105 106 v = matrix * v 107 108 This routine returns a unit vector `axis' parallel to the axis of 109 the input rotation matrix. Moreover for any vector `w' orthogonal 110 to the axis of the rotation, the two vectors 111 112 axis, 113 w x (matrix*w) 114 115 (where "x" denotes the cross product operation) 116 117 will be positive scalar multiples of one another (at least 118 to within the ability to make such computations with double 119 precision arithmetic, and under the assumption that `matrix' 120 does not represent a rotation by zero or pi radians). 121 122 The angle returned will be the angle between `w' and matrix*w 123 for any vector orthogonal to `axis'. 124 125 If the input matrix is a rotation by 0 or pi radians some 126 choice must be made for the axis returned. In the case of 127 a rotation by 0 radians, `axis' is along the positive z-axis. 128 In the case of a rotation by 180 degrees, two choices are 129 possible. The choice made this routine is unspecified. 130 131 -Examples 132 133 This routine can be used to numerically approximate the 134 instantaneous angular velocity vector of a rotating object. 135 136 Suppose that r(t) is the rotation matrix whose columns 137 represent the inertial pointing vectors of the bodyfixed 138 axes of an object at time t. 139 140 Then the angular velocity vector points along the vector 141 given by: 142 T 143 limit axis( r(t+h)r ) 144 h-->0 145 146 And the magnitude of the angular velocity at time t is given by: 147 148 T 149 d angle ( r(t+h)r(t) ) 150 ---------------------- at h = 0 151 dh 152 153 Thus to approximate the angular velocity vector the following 154 code fragment will do 155 156 [ Load t into the double precision variable t 157 Load h into the double precision variable h 158 Load r(t+h) into the 3 by 3 double precision array rth 159 Load r(t) into the 3 by 3 double precision array rt 160 . 161 . 162 . 163 ] 164 165 /. 166 T 167 Compute the infinitesimal rotation r(t+h)r(t) 168 ./ 169 mxmt_c ( rth, rt, infrot ); 170 171 /. 172 Compute the axis and angle of the infinitesimal rotation. 173 /. 174 raxisa_c ( infrot, axis, &angle ); 175 176 /. 177 Scale axis to get the angular velocity vector. 178 ./ 179 vscl_c ( angle/h, axis, angvel ); 180 181 182 -Restrictions 183 184 1) If the input matrix is not a rotation matrix but is close enough 185 to pass the tests this routine performs on it, no error will be 186 signaled, but the results may have poor accuracy. 187 188 2) The input matrix is taken to be an object that acts on (rotates) 189 vectors---it is not regarded as a coordinate transformation. To 190 find the axis and angle of a coordinate transformation, input 191 the transpose of that matrix to this routine. 192 193 -Author_and_Institution 194 195 N.J. Bachman (JPL) 196 W.L. Taber (JPL) 197 F.S. Turner (JPL) 198 199 -Literature_References 200 201 None. 202 203 -Version 204 205 -CSPICE Version 1.0.1, 05-JAN-2005 (NJB) (WLT) (FST) 206 207 Various header updates were made to reflect changes 208 made to the underlying SPICELIB Fortran code. 209 Miscellaneous header corrections were made as well. 210 211 -CSPICE Version 1.0.0, 31-MAY-1999 (WLT) (NJB) 212 213 -Index_Entries 214 215 rotation axis of a matrix 216 217 -& 218 */ 219 220 { /* Begin raxisa_c */ 221 222 /* 223 Local variables 224 */ 225 SpiceDouble tmpmat[3][3]; 226 227 228 /* 229 Error free: no error tracing. 230 */ 231 232 /* 233 Transpose the input matrix to put it in column-major order. 234 */ 235 236 xpose_c ( matrix, tmpmat ); 237 238 raxisa_ ( ( doublereal * ) tmpmat, 239 ( doublereal * ) axis, 240 ( doublereal * ) angle ); 241 242 } /* End raxisa_c */ 243 244