Skip to content

Commit

Permalink
Add some r and fix some a
Browse files Browse the repository at this point in the history
  • Loading branch information
helgee committed Feb 9, 2019
1 parent 76ae306 commit d71fb40
Show file tree
Hide file tree
Showing 3 changed files with 216 additions and 46 deletions.
2 changes: 1 addition & 1 deletion src/a.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
165 changes: 163 additions & 2 deletions src/r.jl
Original file line number Diff line number Diff line change
@@ -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)
Expand Down Expand Up @@ -32,14 +64,109 @@ 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)
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 ###
Expand All @@ -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}()
Expand Down Expand Up @@ -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)
Expand Down
95 changes: 52 additions & 43 deletions test/r.jl
Original file line number Diff line number Diff line change
@@ -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]
Expand All @@ -19,45 +18,53 @@
@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])
act3 = reclat((-1.0, 0.0, 0.0))
@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
Expand All @@ -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
Expand All @@ -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"] =#
Expand Down

0 comments on commit d71fb40

Please sign in to comment.