Merge e8ced01 into 75bb7f8

src/r.jl

@@ -1,7 +1,9 @@
 export 
     rav2xf,
-    rotate,
-    recrad
+    reclat,
+    recpgr,
+    recrad,
+    rotate
 
 """
     rav2xf(rot, av)
@@ -30,6 +32,74 @@ function rav2xf(rot, av)
     xform
 end
 
+
+
+
+"""
+    reclat(rectan)
+
+    Convert from rectangular coordinates to latitudinal coordinates.
+
+### Arguments ###
+
+- `rectan`: Rectangular coordinates of a point.
+
+### Output ###
+
+- `rad`: Distance of the point from the origin
+- `lon`: Planetographic longitude of the point (radians).
+- `lat`: Planetographic latitude of the point (radians).
+
+### References ###
+
+- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/reclat_c.html)
+"""
+function reclat(rectan)
+    lon = Ref{Cdouble}()
+    lat = Ref{Cdouble}()
+    rad = Ref{Cdouble}()
+    ccall((:reclat_c, libcspice), Cvoid, (Ptr{SpiceDouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}), rectan, rad, lon, lat)
+    return rad[], lon[], lat[]
+end
+
+
+
+
+
+
+"""
+    recpgr(body, rectan, re, f)
+
+    Convert rectangular coordinates to planetographic coordinates.
+
+### Arguments ###
+
+- `body`: Body with which coordinate system is associated
+- `rectan`: Rectangular coordinates of a point
+- `re`: Equatorial radius of the reference spheroid
+- `f`: flattening coefficient
+
+### Output ###
+
+- `lon`: Planetographic longitude of the point (radians).
+- `lat`: Planetographic 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/recpgr_c.html)
+"""
+function recpgr(body::AbstractString, rectan, re::Float64, f::Float64)
+    lon = Ref{Cdouble}()
+    lat = Ref{Cdouble}()
+    alt = Ref{Cdouble}()
+    ccall((:recpgr_c, libcspice), Cvoid, (Cstring, Ptr{SpiceDouble}, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}), body, rectan, re, f, lon, lat, alt)
+    handleerror()
+    return lon[], lat[], alt[]
+end
+
+
+
 """
     rotate(angle, iaxis)
 

src/s.jl

@@ -9,6 +9,8 @@ export
     spkpos,
     spkw13,
     str2et,
+    subpnt,
+    subslr,
     swpool,
     sxform
 
@@ -128,6 +130,90 @@ function spd()
     ccall((:spd_c, libcspice), Cdouble, ())
 end
 
+
+
+
+
+"""
+    subpnt(method, target, et, fixref, obsrvr, abcorr)
+
+Compute the rectangular coordinates of the sub-observer point on 
+a target body at a specified epoch, optionally corrected for 
+light time and stellar aberration. 
+
+### Arguments ###
+
+- `method`: Computation method. 
+- `target`: Name of target body. 
+- `et`: Epoch in ephemeris seconds past J2000 TDB. 
+- `fixref`: Body-fixed, body-centered target body frame.  
+- `obsrvr`: Name of observing body. 
+- `abcorr`: Aberration correction.
+
+### Output ###
+
+- `spoint`: Sub-solar point on the target body. 
+- `trgepc`: Sub-solar point epoch. 
+- `srfvec`: Vector from observer to sub-solar point.
+
+Returns `cell` with its cardinality set to `card`.
+
+### References ###
+
+- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/subpnt_c.html)
+"""
+function subpnt(method, target, et, fixref, obsrvr; abcorr="NONE")
+    spoint = Array{SpiceDouble}(undef, 3)
+    trgepc = Ref{SpiceDouble}(0)
+    srfvec = Array{SpiceDouble}(undef, 3)
+    ccall((:subpnt_c, libcspice), Cvoid,
+          (Cstring, Cstring, SpiceDouble, Cstring, Cstring, Cstring, Ptr{SpiceDouble}, Ref{SpiceDouble}, Ptr{SpiceDouble}),
+          method, target, et, fixref, abcorr, obsrvr, spoint, trgepc, srfvec) 
+    handleerror()
+    spoint, trgepc[], srfvec
+end
+
+
+"""
+    subslr(method, target, et, fixref, obsrvr, abcorr)
+
+Compute the rectangular coordinates of the sub-solar point on 
+a target body at a specified epoch, optionally corrected for 
+light time and stellar aberration.
+
+### Arguments ###
+
+- `method`: Computation method. 
+- `target`: Name of target body. 
+- `et`: Epoch in ephemeris seconds past J2000 TDB. 
+- `fixref`: Body-fixed, body-centered target body frame. 
+- `obsrvr`: Name of observing body. 
+- `abcorr`: Aberration correction. 
+
+### Output ###
+
+- `spoint`: Sub-solar point on the target body. 
+- `trgepc`: Sub-solar point epoch. 
+- `srfvec`: Vector from observer to sub-solar point.
+
+Returns `cell` with its cardinality set to `card`.
+
+### References ###
+
+- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/subslr_c.html)
+"""
+function subslr(method::AbstractString, target::AbstractString, et::Float64, fixref::AbstractString, obsrvr::AbstractString; abcorr::AbstractString="NONE")
+    spoint = Array{Cdouble}(undef, 3)
+    trgepc = Ref{Cdouble}(0)
+    srfvec = Array{Cdouble}(undef, 3)
+    ccall((:subslr_c, libcspice), Cvoid, (Cstring, Cstring, Cdouble, Cstring, Cstring, Cstring, Ptr{Cdouble}, Ref{Cdouble}, Ptr{Cdouble}), method, target, et, fixref, abcorr, obsrvr, spoint, trgepc, srfvec) 
+    handleerror()
+    return spoint, trgepc[], srfvec
+end
+
+
+
+
 """
     swpool(agent, names)
 

test/r.jl

@@ -25,4 +25,30 @@
             @test act3[i] ≈ exp3[i]
         end
     end
+
+    kclear()
+    furnsh(
+        path(CORE, :lsk),
+        path(CORE, :pck),
+        path(CORE, :spk))
+        radii = bodvrd("MARS", "RADII", 3)
+        flat = (radii[1] - radii[3])/ radii[1]
+        x = [0.0, -2620.678914818178, 2592.408908856967]
+        lon, lat, alt = recpgr("mars", x, radii[1], flat)
+        actual = [rad2deg(lon), rad2deg(lat), alt]
+        expected = [90., 45, 300]
+        @test actual ≈ expected
+    kclear()  
+
+    kclear()
+    furnsh(
+        path(CORE, :lsk),
+        path(CORE, :pck),
+        path(CORE, :spk))
+        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]
 end

test/s.jl

@@ -25,5 +25,64 @@ using LinearAlgebra: I
     @test spkpos("EARTH", 0.0, "J2000", 399) == ([0.0, 0.0, 0.0], 0.0)
     @test spkpos(399, 0.0, "J2000", 399) == ([0.0, 0.0, 0.0], 0.0)
     kclear()
+
+
+    @testset "subslr" begin
+        try
+            furnsh(
+                path(CORE, :lsk),
+                path(CORE, :pck),
+                path(CORE, :spk),
+            ) 
+            et = str2et("2008 aug 11 00:00:00")
+            re, _, rp = bodvrd("MARS", "RADII", 3)
+            f = (re - rp) / re
+            methods = ["INTERCEPT/ELLIPSOID", "NEAR POINT/ELLIPSOID"]
+            expecteds = [
+                [
+                    0.0,
+                    175.8106755102322,
+                    23.668550281477703,
+                    -175.81067551023222,
+                    23.420819936106213,
+                    175.810721536362,
+                    23.42082337182491,
+                    -175.810721536362,
+                    23.42081994605096,
+                ],
+                [
+                    0.0,
+                    175.8106754100492,
+                    23.420823361866685,
+                    -175.81067551023222,
+                    23.175085577910583,
+                    175.81072152220804,
+                    23.420823371828,
+                    -175.81072152220804,
+                    23.420819946054046,
+                ]
+            ]
+            for (expected, method) in zip(expecteds, methods)
+                spoint, trgepc, srfvec = subslr(method, "Mars", et, "IAU_MARS", "Earth", abcorr="LT+S")
+                spglon, spglat, spgalt = recpgr("mars", spoint, re, f)
+
+                @test spgalt ≈ expected[1]
+                @test rad2deg(spglon) ≈ expected[2]
+                @test rad2deg(spglat) ≈ expected[3]
+                spcrad, spclon, spclat = reclat(spoint)
+                @test rad2deg(spclon) ≈ expected[4]
+                @test rad2deg(spclat) ≈ expected[5]
+                sunpos, sunlt = spkpos("sun", trgepc, "iau_mars", "mars", abcorr="LT+S")
+                supgln, supglt, supgal = recpgr("mars", sunpos, re, f)
+                @test rad2deg(supgln) ≈ expected[6]
+                @test rad2deg(supglt) ≈ expected[7]
+                supcrd, supcln, supclt = reclat(sunpos)
+                @test rad2deg(supcln) ≈ expected[8]
+                @test rad2deg(supclt) ≈ expected[9]
+            end
+        finally
+            kclear()
+        end
+    end
 end