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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
//
// GENERATED FILE
//
use super::*;
use crate::SpiceContext;
use f2rust_std::*;
/// Axis and angle to rotation
///
/// Construct a rotation matrix that rotates vectors by a specified
/// angle about a specified axis.
///
/// # Required Reading
///
/// * [ROTATION](crate::required_reading::rotation)
///
/// # Brief I/O
///
/// ```text
/// VARIABLE I/O DESCRIPTION
/// -------- --- --------------------------------------------------
/// AXIS I Rotation axis.
/// ANGLE I Rotation angle, in radians.
/// R O Rotation matrix corresponding to AXIS and ANGLE.
/// ```
///
/// # Detailed Input
///
/// ```text
/// AXIS,
/// ANGLE are, respectively, a rotation axis and a rotation
/// angle. AXIS and ANGLE determine a coordinate
/// transformation whose effect on any vector V is to
/// rotate V by ANGLE radians about the vector AXIS.
/// ```
///
/// # Detailed Output
///
/// ```text
/// R is a rotation matrix representing the coordinate
/// transformation determined by AXIS and ANGLE: for
/// each vector V, R*V is the vector resulting from
/// rotating V by ANGLE radians about AXIS.
/// ```
///
/// # Exceptions
///
/// ```text
/// Error free.
///
/// 1) If AXIS is the zero vector, the rotation generated is the
/// identity. This is consistent with the specification of VROTV.
/// ```
///
/// # Particulars
///
/// ```text
/// AXISAR can be thought of as a partial inverse of RAXISA. AXISAR
/// really is a `left inverse': the code fragment
///
/// CALL RAXISA ( R, AXIS, ANGLE )
/// CALL AXISAR ( AXIS, ANGLE, R )
///
/// preserves R, except for round-off error, as long as R is a
/// rotation matrix.
///
/// On the other hand, the code fragment
///
/// CALL AXISAR ( AXIS, ANGLE, R )
/// CALL RAXISA ( R, AXIS, ANGLE )
///
/// preserves AXIS and ANGLE, except for round-off error, only if
/// ANGLE is in the range (0, pi). So AXISAR is a right inverse
/// of RAXISA only over a limited domain.
/// ```
///
/// # Examples
///
/// ```text
/// The numerical results shown for these examples may differ across
/// platforms. The results depend on the SPICE kernels used as
/// input, the compiler and supporting libraries, and the machine
/// specific arithmetic implementation.
///
/// 1) Compute a matrix that rotates vectors by pi/2 radians about
/// the Z-axis, and compute the rotation axis and angle based on
/// that matrix.
///
///
/// Example code begins here.
///
///
/// PROGRAM AXISAR_EX1
/// IMPLICIT NONE
///
/// C
/// C SPICELIB functions.
/// C
/// DOUBLE PRECISION DPR
/// DOUBLE PRECISION HALFPI
///
/// C
/// C Local variables
/// C
/// DOUBLE PRECISION ANGLE
/// DOUBLE PRECISION ANGOUT
/// DOUBLE PRECISION AXIS ( 3 )
/// DOUBLE PRECISION AXOUT ( 3 )
/// DOUBLE PRECISION ROTMAT ( 3, 3 )
///
/// INTEGER I
/// INTEGER J
///
/// C
/// C Define an axis and an angle for rotation.
/// C
/// AXIS(1) = 0.D0
/// AXIS(2) = 0.D0
/// AXIS(3) = 1.D0
/// ANGLE = HALFPI()
///
/// C
/// C Determine the rotation matrix.
/// C
/// CALL AXISAR ( AXIS, ANGLE, ROTMAT )
///
/// C
/// C Now calculate the rotation axis and angle based on
/// C ROTMAT as input.
/// C
/// CALL RAXISA ( ROTMAT, AXOUT, ANGOUT )
///
/// C
/// C Display the results.
/// C
/// WRITE(*,'(A)') 'Rotation matrix:'
/// WRITE(*,*)
/// DO I = 1, 3
/// WRITE(*,'(3F10.5)') ( ROTMAT(I,J), J=1,3 )
/// END DO
/// WRITE(*,*)
/// WRITE(*,'(A,3F10.5)') 'Rotation axis :', AXOUT
/// WRITE(*,'(A,F10.5)') 'Rotation angle (deg):',
/// . ANGOUT * DPR()
///
/// END
///
///
/// When this program was executed on a Mac/Intel/gfortran/64-bit
/// platform, the output was:
///
///
/// Rotation matrix:
///
/// 0.00000 -1.00000 0.00000
/// 1.00000 0.00000 0.00000
/// 0.00000 0.00000 1.00000
///
/// Rotation axis : 0.00000 0.00000 1.00000
/// Rotation angle (deg): 90.00000
///
///
/// 2) Linear interpolation between two rotation matrices.
///
/// Let R(t) be a time-varying rotation matrix; R could be
/// a C-matrix describing the orientation of a spacecraft
/// structure. Given two points in time t1 and t2 at which
/// R(t) is known, and given a third time t3, where
///
/// t1 < t3 < t2,
///
/// we can estimate R(t3) by linear interpolation. In other
/// words, we approximate the motion of R by pretending that
/// R rotates about a fixed axis at a uniform angular rate
/// during the time interval [t1, t2]. More specifically, we
/// assume that each column vector of R rotates in this
/// fashion. This procedure will not work if R rotates through
/// an angle of pi radians or more during the time interval
/// [t1, t2]; an aliasing effect would occur in that case.
///
///
/// Example code begins here.
///
///
/// PROGRAM AXISAR_EX2
/// IMPLICIT NONE
///
/// C
/// C SPICELIB functions.
/// C
/// DOUBLE PRECISION DPR
/// DOUBLE PRECISION HALFPI
///
/// C
/// C Local variables
/// C
/// DOUBLE PRECISION ANGLE
/// DOUBLE PRECISION AXIS ( 3 )
/// DOUBLE PRECISION DELTA ( 3, 3 )
/// DOUBLE PRECISION FRAC
/// DOUBLE PRECISION Q ( 3, 3 )
/// DOUBLE PRECISION R1 ( 3, 3 )
/// DOUBLE PRECISION R2 ( 3, 3 )
/// DOUBLE PRECISION R3 ( 3, 3 )
/// DOUBLE PRECISION T1
/// DOUBLE PRECISION T2
/// DOUBLE PRECISION T3
///
/// INTEGER I
/// INTEGER J
///
/// C
/// C Lets assume that R(t) is the matrix that rotates
/// C vectors by pi/2 radians about the Z-axis every
/// C minute.
/// C
/// C Let
/// C
/// C R1 = R(t1), for t1 = 0", and
/// C R2 = R(t2), for t1 = 60".
/// C
/// C Define both matrices and times.
/// C
/// AXIS(1) = 0.D0
/// AXIS(2) = 0.D0
/// AXIS(3) = 1.D0
///
/// T1 = 0.D0
/// T2 = 60.D0
/// T3 = 30.D0
///
/// CALL IDENT ( R1 )
/// CALL AXISAR ( AXIS, HALFPI(), R2 )
///
/// C
/// C Lets compute
/// C
/// C -1
/// C Q = R2 * R1 ,
/// C
/// C The rotation axis and angle of Q define the rotation
/// C that each column of R(t) undergoes from time `t1' to
/// C time `t2'. Since R(t) is orthogonal, we can find Q
/// C using the transpose of R1. We find the rotation axis
/// C and angle via RAXISA.
///
/// CALL MXMT ( R2, R1, Q )
/// CALL RAXISA ( Q, AXIS, ANGLE )
///
/// C
/// C Find the fraction of the total rotation angle that R
/// C rotates through in the time interval [t1, t3].
/// C
/// FRAC = ( T3 - T1 ) / ( T2 - T1 )
///
/// C
/// C Finally, find the rotation DELTA that R(t) undergoes
/// C during the time interval [t1, t3], and apply that
/// C rotation to R1, yielding R(t3), which we'll call R3.
/// C
/// CALL AXISAR ( AXIS, FRAC * ANGLE, DELTA )
/// CALL MXM ( DELTA, R1, R3 )
///
/// C
/// C Display the results.
/// C
/// WRITE(*,'(A,F10.5)') 'Time (s) :', T3
/// WRITE(*,'(A,3F10.5)') 'Rotation axis :', AXIS
/// WRITE(*,'(A,F10.5)') 'Rotation angle (deg):',
/// . FRAC * ANGLE * DPR()
/// WRITE(*,'(A)') 'Rotation matrix :'
/// WRITE(*,*)
/// DO I = 1, 3
/// WRITE(*,'(3F10.5)') ( R3(I,J), J=1,3 )
/// END DO
///
///
/// END
///
///
/// When this program was executed on a Mac/Intel/gfortran/64-bit
/// platform, the output was:
///
///
/// Time (s) : 30.00000
/// Rotation axis : 0.00000 0.00000 1.00000
/// Rotation angle (deg): 45.00000
/// Rotation matrix :
///
/// 0.70711 -0.70711 0.00000
/// 0.70711 0.70711 0.00000
/// 0.00000 0.00000 1.00000
/// ```
///
/// # Author and Institution
///
/// ```text
/// N.J. Bachman (JPL)
/// J. Diaz del Rio (ODC Space)
/// W.L. Taber (JPL)
/// ```
///
/// # Version
///
/// ```text
/// - SPICELIB Version 1.2.0, 06-JUL-2021 (JDR)
///
/// Added IMPLICIT NONE statement.
///
/// Edited the header to comply with NAIF standard. Removed
/// unnecessary $Revisions section.
///
/// Added complete code examples based on existing code fragments.
///
/// - SPICELIB Version 1.1.0, 25-AUG-2005 (NJB)
///
/// Updated to remove non-standard use of duplicate arguments
/// in VROTV call.
///
/// Identity matrix is now obtained from IDENT.
///
/// - SPICELIB Version 1.0.1, 10-MAR-1992 (WLT)
///
/// Comment section for permuted index source lines was added
/// following the header.
///
/// - SPICELIB Version 1.0.0, 30-AUG-1990 (NJB)
/// ```
pub fn axisar(axis: &[f64; 3], angle: f64, r: &mut [[f64; 3]; 3]) {
AXISAR(axis, angle, r.as_flattened_mut());
}
//$Procedure AXISAR ( Axis and angle to rotation )
pub fn AXISAR(AXIS: &[f64], ANGLE: f64, R: &mut [f64]) {
let AXIS = DummyArray::new(AXIS, 1..=3);
let mut R = DummyArrayMut2D::new(R, 1..=3, 1..=3);
let mut VTEMP = StackArray::<f64, 3>::new(1..=3);
//
// Local variables
//
//
// First, set R equal to the identity.
//
IDENT(R.as_slice_mut());
//
// The matrix we want rotates EVERY vector by ANGLE about AXIS.
// In particular, it does so to our basis vectors. The columns
// of R are the images of the basis vectors under this rotation.
//
for I in 1..=3 {
VROTV(
R.subarray([1, I]),
AXIS.as_slice(),
ANGLE,
VTEMP.as_slice_mut(),
);
VEQU(VTEMP.as_slice(), R.subarray_mut([1, I]));
}
}