From d71fb401ed11b1b607e64a5eea2378ab52528fa2 Mon Sep 17 00:00:00 2001 From: Helge Eichhorn Date: Sat, 9 Feb 2019 08:06:53 +0100 Subject: [PATCH] Add some r and fix some a --- src/a.jl | 2 +- src/r.jl | 165 +++++++++++++++++++++++++++++++++++++++++++++++++++++- test/r.jl | 95 +++++++++++++++++-------------- 3 files changed, 216 insertions(+), 46 deletions(-) diff --git a/src/a.jl b/src/a.jl index acc2a9b..70f8245 100644 --- a/src/a.jl +++ b/src/a.jl @@ -21,5 +21,5 @@ Rotation matrix corresponding to `axis` and `angle` function axisar(axis, angle) r = Matrix{Float64}(undef, 3,3) ccall((:axisar_c, libcspice), Cvoid, (Ptr{SpiceDouble}, SpiceDouble, Ptr{SpiceDouble}), axis, angle, r) - return r' + permutedims(r) end diff --git a/src/r.jl b/src/r.jl index 9a6589e..631d1db 100644 --- a/src/r.jl +++ b/src/r.jl @@ -1,10 +1,42 @@ export + radrec, rav2xf, + raxisa, + reccyl, + recgeo, reclat, recpgr, + recsph, recrad, rotate +""" + radrec(range, ra, dec) + +Convert from range, right ascension, and declination to rectangular coordinates. + +### Arguments ### + +range I Distance of a point from the origin. +ra I Right ascension of point in radians. +dec I Declination of point in radians. + +### Output ### + +Returns the rectangular coordinates of the point. + +### References ### + +- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/radrec_c.html) +""" +function radrec(range, ra, dec) + rectan = Array{SpiceDouble}(undef, 3) + ccall((:radrec_c, libcspice), Cvoid, + (SpiceDouble, SpiceDouble, SpiceDouble, Ptr{SpiceDouble}), + range, ra, dec, rectan) + rectan +end + """ rav2xf(rot, av) @@ -32,6 +64,101 @@ function rav2xf(rot, av) permutedims(xform) end +""" + raxisa(matrix) + +Compute the axis of the rotation given by an input matrix and the angle of the rotation +about that axis. + +### Arguments ### + +- `matrix`: A 3x3 rotation matrix + +### Output ### + +- `axis`: Axis of the rotation +- `angle`: Angle through which the rotation is performed + +### References ### + +- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/raxisa_c.html) +""" +function raxisa(matrix) + size(matrix) != (3, 3) && throw(ArgumentError("`matrix` must be a 3x3 matrix.")) + axis = Array{SpiceDouble}(undef, 3) + angle = Ref{SpiceDouble}() + ccall((:raxisa_c, libcspice), Cvoid, + (Ptr{SpiceDouble}, Ptr{SpiceDouble}, Ref{SpiceDouble}), + permutedims(matrix), axis, angle) + handleerror() + axis, angle[] +end + +""" + reccyl(rectan) + +Convert from rectangular to cylindrical coordinates. + +### Arguments ### + +- `rectan`: Rectangular coordinates of a point + +### Output ### + +- `r`: Distance of the point from the Z axis +- `lon`: Angle (radians) of the point from the XZ plane +- `z`: Height of the point above the XY plane + +### References ### + +- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/reccyl_c.html) +""" +function reccyl(rectan) + length(rectan) != 3 && throw(ArgumentError("Length of `rectan` must be 3.")) + r = Ref{SpiceDouble}() + lon = Ref{SpiceDouble}() + z = Ref{SpiceDouble}() + ccall((:reccyl_c, libcspice), Cvoid, + (Ptr{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}), + rectan, r, lon, z) + r[], lon[], z[] +end + + +""" + recgeo(rectan, re, f) + +Convert from rectangular coordinates to geodetic coordinates. + +### Arguments ### + +- `rectan`: Rectangular coordinates of a point +- `re`: Equatorial radius of the reference spheroid +- `f`: Flattening coefficient + +### Output ### + +- `lon`: Geodetic longitude of the point (radians) +- `lat`: Geodetic latitude of the point (radians) +- `alt`: Altitude of the point above reference spheroid + +### References ### + +- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/XXX_c.html) +""" +function recgeo(rectan, re, f) + length(rectan) != 3 && throw(ArgumentError("Length of `rectan` must be 3.")) + lon = Ref{SpiceDouble}() + lat = Ref{SpiceDouble}() + alt = Ref{SpiceDouble}() + ccall((:recgeo_c, libcspice), Cvoid, + (Ptr{SpiceDouble}, SpiceDouble, SpiceDouble, + Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}), + rectan, re, f, lon, lat, alt) + handleerror() + lon[], lat[], alt[] +end + """ reclat(rectan) @@ -39,7 +166,7 @@ Convert from rectangular coordinates to latitudinal coordinates. ### Arguments ### -- `rectan`: Rectangular coordinates of a point as an iterable with three elements. +- `rectan`: Rectangular coordinates of a point ### Output ### @@ -54,6 +181,7 @@ Returns a tuple consisting of: - [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/reclat_c.html) """ function reclat(rectan) + length(rectan) != 3 && throw(ArgumentError("Length of `rectan` must be 3.")) lon = Ref{SpiceDouble}() lat = Ref{SpiceDouble}() rad = Ref{SpiceDouble}() @@ -86,16 +214,49 @@ Convert rectangular coordinates to planetographic coordinates. - [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/recpgr_c.html) """ function recpgr(body, rectan, re, f) + length(rectan) != 3 && throw(ArgumentError("Length of `rectan` must be 3.")) lon = Ref{SpiceDouble}() lat = Ref{SpiceDouble}() alt = Ref{SpiceDouble}() ccall((:recpgr_c, libcspice), Cvoid, - (Cstring, Ptr{SpiceDouble}, SpiceDouble, SpiceDouble, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}), + (Cstring, Ptr{SpiceDouble}, SpiceDouble, SpiceDouble, + Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}), body, rectan, re, f, lon, lat, alt) handleerror() lon[], lat[], alt[] end + +""" + recsph(rectan) + +Convert from rectangular coordinates to spherical coordinates. + +### Arguments ### + +- `rectan`: Rectangular coordinates of a point + +### Output ### + +- `r `: Distance of the point from the origin +- `colat `: Angle of the point from the Z-axis in radian +- `lon `: Longitude of the point in radians + +### References ### + +- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/recsph_c.html) +""" +function recsph(rectan) + length(rectan) != 3 && throw(ArgumentError("Length of `rectan` must be 3.")) + r = Ref{SpiceDouble}() + colat = Ref{SpiceDouble}() + lon = Ref{SpiceDouble}() + ccall((:recsph_c, libcspice), Cvoid, + (Ptr{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}), + rectan, r, colat, lon) + r[], colat[], lon[] +end + """ rotate(angle, iaxis) diff --git a/test/r.jl b/test/r.jl index 9dc4429..a77f873 100644 --- a/test/r.jl +++ b/test/r.jl @@ -1,10 +1,9 @@ @testset "R" begin -#= @testset "radrec" begin =# -#= npt.assert_array_almost_equal([1.0, 0.0, 0.0], radrec(1.0, 0.0, 0.0)) =# -#= npt.assert_array_almost_equal([0.0, 1.0, 0.0], radrec(1.0, 90.0 * rpd(), 0.0)) =# -#= npt.assert_array_almost_equal([0.0, 0.0, 1.0], radrec(1.0, 0.0, 90.0 * rpd())) =# -#= =# -#= =# + @testset "radrec" begin + @test [1.0, 0.0, 0.0] ≈ radrec(1.0, 0.0, 0.0) + @test [0.0, 1.0, 0.0] ≈ radrec(1.0, π/2, 0.0) + @test [0.0, 0.0, 1.0] ≈ radrec(1.0, 0.0, π/2) + end @testset "rav2xf" begin e = [1.0, 0.0, 0.0] rz = [0.0 1.0 0.0; -1.0 0.0 0.0; 0.0 0.0 1.0] @@ -19,38 +18,45 @@ @test exp[i] ≈ act[i] end end -#= @testset "raxisa" begin =# -#= axis = [1.0, 2.0, 3.0] =# -#= angle = 0.1 * twopi() =# -#= rotate_matrix = axisar(axis, angle) =# -#= axout, angout = raxisa(rotate_matrix) =# -#= expectedAngout = [0.26726124, 0.53452248, 0.80178373] =# -#= npt.assert_approx_equal(angout, 0.62831853, significant=7) =# -#= npt.assert_array_almost_equal(axout, expectedAngout) =# -#= =# -#= =# -#= @testset "reccyl" begin =# -#= expected1 = array([0.0, 0.0, 0.0]) =# -#= expected2 = array([1.0, 90.0 * rpd(), 0.0]) =# -#= expected3 = array([1.0, 270.0 * rpd(), 0.0]) =# -#= npt.assert_array_almost_equal(expected1, reccyl([0.0, 0.0, 0.0]), decimal=7) =# -#= npt.assert_array_almost_equal(expected2, reccyl([0.0, 1.0, 0.0]), decimal=7) =# -#= npt.assert_array_almost_equal(expected3, reccyl([0.0, -1.0, 0.0]), decimal=7) =# -#= =# -#= =# -#= @testset "recgeo" begin =# -#= kclear() =# -#= furnsh(CoreKernels.testMetaKernel) =# -#= num_vals, radii = bodvrd("EARTH", "RADII", 3) =# -#= flat = (radii[0] - radii[2]) / radii[0] =# -#= x = [-2541.748162, 4780.333036, 3360.428190] =# -#= lon, lat, alt = recgeo(x, radii[0], flat) =# -#= actual = [lon * dpr(), lat * dpr(), alt] =# -#= expected = [118.000000, 32.000000, 0.001915518] =# -#= npt.assert_array_almost_equal(actual, expected, decimal=4) =# -#= kclear() =# -#= =# -#= =# + @testset "raxisa" begin + axis = [1.0, 2.0, 3.0] + angle = 0.2π + rotate_matrix = axisar(axis, angle) + axout, angout = raxisa(rotate_matrix) + expected_angout = [0.26726124, 0.53452248, 0.80178373] + @test angout ≈ 0.62831853 + @test axout ≈ expected_angout + @test_throws SpiceError raxisa(randn(3, 3)) + @test_throws ArgumentError raxisa(randn(2, 2)) + end + @testset "reccyl" begin + expected1 = (0.0, 0.0, 0.0) + expected2 = (1.0, π/2, 0.0) + expected3 = (1.0, 3/2 * π, 0.0) + @test all(reccyl([0.0, 0.0, 0.0]) .≈ expected1) + @test all(reccyl([0.0, 1.0, 0.0]) .≈ expected2) + @test all(reccyl([0.0, -1.0, 0.0]) .≈ expected3) + @test_throws ArgumentError reccyl([0.0, -1.0]) + end + @testset "recgeo" begin + try + furnsh(path(CORE, :pck)) + radii = bodvrd("EARTH", "RADII", 3) + flat = (radii[1] - radii[3]) / radii[1] + x = [-2541.748162, 4780.333036, 3360.428190] + lon, lat, alt = recgeo(x, radii[1], flat) + actual = [rad2deg(lon), rad2deg(lat), alt] + expected = [118.0, 32.0, 0.001915518] + @testset for i in eachindex(actual, expected) + @test actual[i] ≈ expected[i] rtol=1e-4 + end + @test_throws ArgumentError recgeo(x[1:2], radii[1], flat) + @test_throws SpiceError recgeo(x, -radii[1], flat) + @test_throws SpiceError recgeo(x, radii[1], 1.0) + finally + kclear() + end + end @testset "reclat" begin act1 = reclat([1.0, 0.0, 0.0]) act2 = reclat([0.0, 1.0, 0.0]) @@ -58,6 +64,7 @@ @test [act1[1], act1[2], act1[3]] ≈ [1.0, 0.0, 0.0] @test [act2[1], act2[2], act2[3]] ≈ [1.0, deg2rad(90.0), 0.0] @test [act3[1], act3[2], act3[3]] ≈ [1.0, deg2rad(180.0), 0.0] + @test_throws ArgumentError reclat([1.0, 2.0]) end @testset "recpgr" begin try @@ -69,6 +76,9 @@ actual = [rad2deg(lon), rad2deg(lat), alt] expected = [90., 45, 300] @test actual ≈ expected + @test_throws ArgumentError recpgr("MARS", x[1:2], radii[1], flat) + @test_throws SpiceError recpgr("MARS", x, -radii[1], flat) + @test_throws SpiceError recpgr("MARS", x, radii[1], 1.0) finally kclear() end @@ -90,11 +100,10 @@ @test act3[i] ≈ exp3[i] end end -#= @testset "recsph" begin =# -#= v1 = array([-1.0, 0.0, 0.0]) =# -#= assert recsph(v1) == (1.0, pi/2, pi) =# -#= =# -#= =# + @testset "recsph" begin + v1 = [-1.0, 0.0, 0.0] + @test all(recsph(v1) .≈ (1.0, π/2, π)) + end #= @testset "removc" begin =# #= cell = cell_char(10, 10) =# #= items = ["one", "two", "three", "four"] =#