Merge 7d49578 into 95511f6

.gitignore

@@ -16,3 +16,4 @@ deps/windows/usr/
 deps/windows/x64
 docs/build/
 test/kernels/
+deps/build.log

src/e.jl

@@ -0,0 +1,28 @@
+export eul2m
+
+"""
+    eul2m(angle3, angle2, angle1, axis3, axis2, axis1)
+
+Construct a rotation matrix from a set of Euler angles.
+
+### Arguments ###
+
+- `angle3`, `angle2`, `angle1`: Rotation angles about third, second, and first rotation axes (radians)
+- `axis3`, `axis2`, `axis1`: Axis numbers of third, second, and first rotation axes
+
+### Output ###
+
+A rotation matrix corresponding to the product of the 3 rotations.
+
+### References ###
+
+- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/eul2m_c.html)
+"""
+function eul2m(angle3, angle2, angle1, axis3, axis2, axis1)
+    r = Matrix{SpiceDouble}(undef, 3, 3)
+    ccall((:eul2m_c, libcspice), Cvoid,
+          (SpiceDouble, SpiceDouble, SpiceDouble, SpiceInt, SpiceInt, SpiceInt, Ptr{SpiceDouble}),
+          angle3, angle2, angle1, axis3, axis2, axis1, r)
+    handleerror()
+    permutedims(r)
+end

src/m.jl

@@ -1,6 +1,8 @@
-export mxvg
+export 
+    mxvg,
+    m2q
 
-function mxvg(m1, v2)
+function _mxvg(m1, v2)
     lm1, lm2 = size(m1)
     lv = length(v2)
     if lm2 != lv
@@ -12,3 +14,37 @@ function mxvg(m1, v2)
     handleerror()
     return vout
 end
+
+"""
+    mxvg(m1,v2)
+
+**Deprecated:** Use `m1 * v2` instead.
+"""
+mxvg
+
+@deprecate mxvg(m1, v2) m1 * v2
+
+"""
+    m2q(r)
+
+Find a unit quaternion corresponding to a specified rotation matrix.
+
+### Arguments ###
+
+- `r`: A rotation matrix
+
+### Output ###
+
+A unit quaternion representing `r'
+
+### References ###
+
+- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/m2q_c.html)
+"""
+function m2q(r)
+    q = Array{SpiceDouble}(undef, 4)
+    ccall((:m2q_c, libcspice), Cvoid, (Ptr{Float64}, Ptr{Float64}), permutedims(r), q)
+    handleerror()
+    q
+end
+

src/q.jl

@@ -0,0 +1,83 @@
+export 
+    q2m,
+    qxq,
+    qdq2av
+
+"""
+    q2m(q...)
+
+Find the rotation matrix corresponding to a specified unit quaternion.
+
+### Arguments ###
+
+- `q`: A unit quaternion (as any kind of iterable with four elements)
+
+### Output ###
+
+A rotation matrix corresponding to `q`.
+
+### References ###
+
+- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/q2m_c.html)
+"""
+q2m(q...) = q2m(collect(q))
+
+function q2m(q)
+    length(q) != 4 && throw(ArgumentError("`q` needs to be an iterable with four elements."))    
+    r = Matrix{SpiceDouble}(undef, 3, 3)
+    ccall((:q2m_c, libcspice), Cvoid, (Ptr{SpiceDouble}, Ptr{SpiceDouble}), collect(q), r)
+    permutedims(r)
+end
+
+"""
+    qxq(q1,q2)
+
+Multiply two quaternions. 
+
+### Arguments ###
+
+- `q1`: First SPICE quaternion factor (as any kind of iterable with four elements)
+- `q2`: Second SPICE quaternion factor (as any kind of iterable with four elements)
+
+### Output ###
+
+A quaternion corresponding to the product of `q1' and `q2'
+
+### References ###
+
+- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/qxq_c.html)
+"""
+function qxq(q1, q2)
+    length(q1) != 4 && throw(ArgumentError("`q1` needs to be an iterable with four elements.")) 
+    length(q2) != 4 && throw(ArgumentError("`q2` needs to be an iterable with four elements."))   
+    q = Array{SpiceDouble}(undef, 4)
+    ccall((:qxq_c, libcspice), Cvoid, (Ptr{SpiceDouble}, Ptr{SpiceDouble}, Ptr{SpiceDouble}), collect(q1), collect(q2), q)
+    q
+end
+
+"""
+    qdq2av(q,dq)
+
+Derive angular velocity from a unit quaternion and its derivative 
+with respect to time. 
+
+### Arguments ###
+
+- `q`: Unit SPICE quaternion (as any kind of iterable with four elements)
+- `dq`: Derivative of `q' with respect to time
+
+### Output ###
+
+Angular velocity vector defined by `q' and `dq'
+
+### References ###
+
+- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/qdq2av_c.html)
+"""
+function qdq2av(q,dq)
+    length(q) != 4 && throw(ArgumentError("`q` needs to be an iterable with four elements.")) 
+    length(dq) != 4 && throw(ArgumentError("`dq` needs to be an iterable with four elements."))   
+    av = Array{SpiceDouble}(undef, 3)
+    ccall((:qdq2av_c, libcspice), Cvoid, (Ptr{SpiceDouble}, Ptr{SpiceDouble}, Ptr{SpiceDouble}), collect(q), collect(dq), av)
+    av
+end

src/r.jl

@@ -1,4 +1,6 @@
-export rav2xf
+export 
+    rav2xf,
+    rotate
 
 """
     rav2xf(rot, av)
@@ -27,3 +29,29 @@ function rav2xf(rot, av)
     xform
 end
 
+"""
+    rotate(angle, iaxis)
+
+ Calculate the 3x3 rotation matrix generated by a rotation 
+of a specified angle about a specified axis. This rotation 
+is thought of as rotating the coordinate system. 
+ 
+### Arguments ###
+
+- `angle`: Angle of rotation (radians)
+- `iaxis`: Axis of rotation (X=1, Y=2, Z=3)
+
+### Output ###
+
+Returns rotation matrix associated with `angle` and `iaxis`.
+
+### References ###
+
+- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/rotate_c.html)
+"""
+function rotate(angle, iaxis)
+    r = Matrix{SpiceDouble}(undef, 3, 3)
+    ccall((:rotate_c, libcspice), Cvoid, (SpiceDouble, SpiceInt, Ptr{SpiceDouble}), angle, iaxis, r)
+    permutedims(r)
+end
+

test/e.jl

@@ -0,0 +1,12 @@
+@testset "E" begin
+    let 
+        act = eul2m(3, 2, 1, 3, 2, 1)
+        @test size(act) == (3, 3)
+        exp = [0.411982245665683 -0.6812427202564033 0.6051272472413688;
+               0.05872664492762098 -0.642872836134547 -0.7637183366502791;
+               0.9092974268256817 0.35017548837401463  -0.2248450953661529]
+        @testset for i in eachindex(act, exp)
+            @test act[i] ≈ exp[i]
+        end
+    end
+end

test/m.jl

@@ -1,6 +1,16 @@
 @testset "M" begin 
-    matrix = [1. 1. 1. 
+    let 
+        matrix = [1. 1. 1. 
                 2. 3. 4.]
-    vecgood = [1., 2., 3.]
-    @test mxvg(matrix, vecgood) == [6.,20.]
-end 
+        vecgood = [1., 2., 3.]
+        @test SPICE._mxvg(matrix, vecgood) == [6.,20.]
+    end
+    let 
+        r = rotate(π/2, 3)
+        act = m2q(r)
+        exp = [sqrt(2)/2.0, 0.0, 0.0, -sqrt(2)/2.0]
+        @testset for i in eachindex(act, exp)
+            @test act[i] ≈ exp[i]
+        end
+    end
+end

test/q.jl

@@ -0,0 +1,69 @@
+@testset "Q" begin
+    let exp = [0.607843137254902 0.27450980392156854 0.7450980392156862;
+                0.6666666666666666 0.33333333333333326 -0.6666666666666666;
+                -0.43137254901960775 0.9019607843137255 0.019607843137254832]
+        act = q2m(0.5, 0.4, 0.3, 0.1)
+        @testset for i in eachindex(act, exp)
+            @test act[i] ≈ exp[i]
+        end
+        q = [0.5, 0.4, 0.3, 0.1]
+        act = q2m(q)
+        @testset for i in eachindex(act, exp)
+            @test act[i] ≈ exp[i]
+        end
+        q = (0.5, 0.4, 0.3, 0.1)
+        act = q2m(q)
+        @testset for i in eachindex(act, exp)
+            @test act[i] ≈ exp[i]
+        end
+    end
+    let angle = deg2rad.([-20.0, 50.0, -60.0])
+        m = eul2m(angle[3], angle[2], angle[1], 3, 1, 3)
+        q = m2q(m)
+        exp = [1.0, 2.0, 3.0]
+        qav = [0.0, 1.0, 2.0, 3.0]
+        dq = -0.5*qxq(q, qav)
+        act = qdq2av(q, dq)
+        @testset for i in eachindex(act, exp)
+            @test act[i] ≈ exp[i]
+        end
+    end
+    let qID = [1.0, 0.0, 0.0, 0.0]
+        nqID = [-1.0, 0.0, 0.0, 0.0]
+        qI = [0.0, 1.0, 0.0, 0.0]
+        qJ = [0.0, 0.0, 1.0, 0.0]
+        qK = [0.0, 0.0, 0.0, 1.0]
+        qIJ = qxq(qI, qJ)
+        @testset for i in eachindex(qIJ, qK)
+            @test qIJ[i] ≈ qK[i]
+        end
+        qJK = qxq(qJ, qK)
+        @testset for i in eachindex(qJK, qI)
+            @test qJK[i] ≈ qI[i]
+        end
+        qKI = qxq(qK, qI)
+        @testset for i in eachindex(qKI, qJ)
+            @test qKI[i] ≈ qJ[i]
+        end
+        qII = qxq(qI, qI)
+        @testset for i in eachindex(qII, nqID)
+            @test qII[i] ≈ nqID[i]
+        end
+        qJJ = qxq(qJ, qJ)
+        @testset for i in eachindex(qJJ, nqID)
+            @test qJJ[i] ≈ nqID[i]
+        end
+        qKK = qxq(qK, qK)
+        @testset for i in eachindex(qKK, nqID)
+            @test qKK[i] ≈ nqID[i]
+        end
+        qIDI = qxq(qID, qI)
+        @testset for i in eachindex(qIDI, qI)
+            @test qIDI[i] ≈ qI[i]
+        end
+        qIID = qxq(qI, qID)
+        @testset for i in eachindex(qIID, qI)
+            @test qIID[i] ≈ qI[i]
+        end 
+    end
+end

test/r.jl

@@ -0,0 +1,11 @@
+@testset "R" begin 
+    let 
+        act = rotate(π/4, 3)
+        exp = [sqrt(2)/2.0 sqrt(2)/2.0 0.0;
+              -sqrt(2)/2.0 sqrt(2)/2.0 0.0;
+               0.0 0.0 1.0]
+        @testset for i in eachindex(act, exp)
+             @test act[i] ≈ exp[i]
+        end
+    end
+end