Coordinate conversion functions

[chris]

There has been a lot of confusion lately about how to get the position of an object in a celx script. I've written some well-documented functions to make this task easier. I will be expanding on this utility library, but for now, there are four functions:

position_ecl_j2000 -- get the position of an object in the J2000 ecliptic coordinate system
position_eq_j2000 -- get the position of an object in the J2000 equatorial coordinate system
rect_to_spherical -- convert rectangular coordinates to spherical coordinates
spherical_to_rect -- convert spherical coordinates to rectangular coordinates

Other people have their own versions of these functions. I've written mine with clarity as the main objective, so that you can see exactly the transformations required when dealing with coordinates in Celestia.

Here is the code for the four functions:

Code: Select all

-- Constants required for the coordinate conversion functions
obliquity_j2000 = math.rad(-23.4392911)
cos_obliquity = math.cos(obliquity_j2000)
sin_obliquity = math.sin(obliquity_j2000)

-- position_ecl_j2000
--
-- Get the position of the object obj relative to the center object
-- at time t. The position is the J2000 ecliptic coordinate system,
-- which is the basis of Celestia's native coordinate system.
--
-- Return value is a Celestia vector with the coordinates of obj
-- in kilometers.
--
-- Parameters:
--   obj    : Celestia object
--   center : Celestia object
--   t      : time (TDB)
--
function position_ecl_j2000(obj, center, t)
    local uly_to_km = 9460730.4725808
    local p0 = center:getposition(t)
    local p = obj:getposition(t)

    -- Convert to kilometers
    local v = (p - p0) * uly_to_km

    -- Convert from Celestia's internal coordinate system
    return celestia:newvector(v.x, -v.z, v.y)
end


-- position_eq_j2000
--
-- Get the position of the object obj relative to the center object
-- at time t. The position is the J2000 Earth equatorial coordinate
-- system.
--
-- Parameters:
--   obj    : Celestia object
--   center : Celestia object
--   t      : time (TDB)
--
-- Return value is a Celestia vector with the coordinates of obj
-- in kilometers.
--
function position_eq_j2000(obj, center, t)
    -- Get the object position in ecliptic coordinates
    local v = position_ecl_j2000(obj, center, t)

    -- Rotate about the x-axis by the Earth's obliquity at J2000.0
    return celestia:newvector(v.x,
                              v.y * cos_obliquity + v.z * sin_obliquity,
                              v.y * -sin_obliquity + v.z * cos_obliquity)
end


-- rect_to_spherical
--
-- Convert rectangular coordinates to spherical coordinates.
--
-- Parameters:
--   v :     Celestia vector
--
-- Return values:
--   altitude, azimuth, distance
--   altitude and azimuth are both in radians. Altitude will be between
--   +pi/2 and -pi/2; azimuth will be the range -pi to +pi.
--
function rect_to_spherical(v)
    local r = v:length()
    local altitude = math.asin(v.z / r)
    local azimuth = math.atan2(v.y, v.x)
    return altitude, azimuth, r
end


-- spherical_to_rect
--
-- Convert spherical coordinates to rectangular coordinates.
--
-- Parameters:
--   altitude : altitude angle in radians
--   azimuth  : azimuth angle in radians
--   distance
--
-- Return value:
--   A Celestia vector with rectangular coordinates.
--
function spherical_to_rect(altitude, azimuth, distance)
    local s = math.cos(altitude)
    local x = math.cos(azimuth) * s * distance
    local y = math.sin(azimuth) * s * distance
    local z = math.sin(altitude) * distance
    return celestia:newvector(x, y, z)
end


And here is a very simple example usage that shows the geocentric right ascension and declination of the currently selected object:

Code: Select all

-- Test the functions; show the geocentric RA and Declination of the
-- currently selected object.
sun = celestia:find("Sol")
earth = celestia:find("Sol/Earth")

while true do
    local obj = celestia:getselection()
    if obj then
        local pos = position_eq_j2000(obj, earth, celestia:gettime())
        local alt, azimuth, r = rect_to_spherical(pos)

        -- Convert the spherical coordinates in radians to hours (RA)
        -- and degrees (Declination)
        local ra = math.deg(azimuth) / 15
        if ra < 0 then ra = ra + 24 end
        local dec = math.deg(alt)

        celestia:flash(string.format("RA: %f Dec: %f", ra, dec))
    end
   
    wait(0)
end


--Chris

[BobHegwood]

As always, Chris, thanks very much from the Brain-Dead.
Any clarifications for CELX are much appreciated here.

[chris]

It is often useful to be able to compute the location of one object projected onto the surface of another. This point on the surface is known as the subpoint. For example, we may want to know the latitude and longitude of the point directly beneath an Earth orbiting spacecraft. The following code defines a function called subpoint which is used to show the current latitude and longitude of the observer:

Code: Select all

-- subpoint
--
-- Return the longitude   and latitude of   the point p projected onto the surface
-- of the solar   system body obj
--
-- Parameters:
--  obj      : a Celestia   object which must be a solar system body (i.e. the
--            function doesn't work for stars and deep sky objects.)
--  p      : a Celestia   position in universal coordinates
--
-- Return values:
--  latitude, longitude   : in radians
--
-- Note:
--  This function only returns the true   subpoint for spherical objects
--
function subpoint(obj, p)
    local f = obj:bodyfixedframe()
    if f then
        local p1 = f:to(p)

        -- Convert to standard coordinates
        local v = celestia:newvector(p1.x, -p1.z, p1.y)
        local lat, long = rect_to_spherical(v)
        return lat, long
    else
        return 0, 0
    end
end


Code: Select all

earth = celestia:find("Sol/Earth")

while true do
    local obs = celestia:getobserver()
    local long, lat = subpoint(earth, obs:getposition())
    celestia:flash(string.format("Lat: %f Long: %f", math.deg(lat), math.deg(long)))
    wait(0)
end


The code for subpoint shows how using reference frames properly can make some code very elegant. The 'magic' of the function is this code that converts a point in Celestia's universal coordinate system to the body fixed coordinate system of obj:

Code: Select all

local f = obj:bodyfixedframe()
local p1 = f:to(p)


After this step, we get the latitude and longitude simply by converting from rectangular to spherical coordinates.

Two caveats:
- The function subpoint only returns the true subpoint for spherical objects. For oblate objects, it is necessary to compute the closest point on an ellipsoid, which is much more complex than computing the closest point on a sphere.
- subpoint will only work on Celestia 1.6.0 because it uses the new function bodyfixedframe.

--Chris

[eric_cannes]

As a beginner in Celestia development (scriptedorbit and scriptedrotation), I have some questions about the used frame and the coordinate conversion functions.

For visualising my satellite around earth, I used J2000 equatorial x,y,z coordinate, and also a quaternion orientation from J2000 eq to to the Satellite body frame.
I tried first the scriptedrotation method with a unitary quaternion (1, 0, 0, 0), it means, no angular rotations around X, Y or Z.
I was surprised to observe that the Satellite "Body frame" is not aligned on the "Frame axes" displayed in celestia view (no rotation).

I found out that X=-X, Y=-Y and Z=Z.

It sounds like the ssc local frame (in my case J2000 eq frame) is different than the displayed "frame axes" of my object satellite.

As I only used J2000 frame and Satellite frame, I would like to see only these 2 frames. I dont know if it is possible.

what is the "frame axes" displayed in celestia ? universal coordinate ? celestia internal coordinate ?
what is the "body axes" displayed in celestia ? J2000 eq frame ?

When using the scriptedrotation, what are the frame axes used ? it is not the J2000 eq frame..
Scriptedrotation is about rotation in which frame ? Celestia's internal coordinate system ?

[chris]

As a beginner in Celestia development (scriptedorbit and scriptedrotation), I have some questions about the used frame and the coordinate conversion functions.

...

what is the "frame axes" displayed in celestia ? universal coordinate ? celestia internal coordinate ?
what is the "body axes" displayed in celestia ? J2000 eq frame ?

The frame axes and body axes are both shown in the Celestia's 'universal coordinate system', which is based on the J2000 ecliptic (not equatorial frame.) The body axes are the frame axes transformed by the objects rotation model (UniformRotation, ScriptedRotation, etc.)

When using the scriptedrotation, what are the frame axes used ? it is not the J2000 eq frame..
Scriptedrotation is about rotation in which frame ? Celestia's internal coordinate system ?

ScriptedRotation is a rotation in the object's body frame; a rotation of (1, 0, 0, 0)--identity--*should* cause the an objects body axes to be aligned with the frame axes.

--Chris

[eric_cannes]

Thank's Chris for your answer but its indeed more confusing than I though.
I've made a very simple test of ScriptedRotation and I observed the result for 4 simple quaternions (id and 3 rotations around X,Y and Z).
Here is my file "rtfifo.lua"
--------------------------------------------
rtAttProto = { }

-- constructor method
function rtAttProto:new(o)
o = o or {}
setmetatable (o, self)
self.__index = self
o.period = 1
return o
end

-- orientation function
function rtAttProto:orientation(tjd)
local t = tjd
-- return a quaternion w, x, y, z
return 1, 0, 0, 0
end

function rtAtt(sscvals)
return rtAttProto:new(sscvals)
end
--------------

If scriptedrotation returns 1,0,0,0 => I observe a rotation of 180 deg around Zaxis.
If scriptedrotation returns 0,1,0,0 => I observe a rotation of 180 deg around Yaxis.
If scriptedrotation returns 0,0,1,0 => I observe no rotation (body axes and frame axes are aligned).
If scriptedrotation returns 0,0,0,1 => I observe a rotation of 180 deg around Xaxis.

Some of my colleagues correct the function (SampleRotation or ScriptedRotation) by adding some extra rotations.
For example, in my example, we observe good results and correct simulations in regards to what we expected.
Here the added lines :

pkg = {
Xaxis = celestia:newvector(1,0,0),
Yaxis = celestia:newvector(0,1,0),
Zaxis = celestia:newvector(0,0,1)
}
pkg.zmPi = celestia:newrotation( pkg.Zaxis, -math.pi )
pkg.xmPi2 = celestia:newrotation( pkg.Xaxis, -math.pi/2.0 )

function convert_quaternion( w, x, y, z)
n = math.sqrt(w*w + x*x + y*y + z*z)
q_sc_i = celestia:newrotation( w/n, x/n, y/n, z/n )
return pkg.xmPi2 * q_sc_i * pkg.xmPi2 * pkg.zmPi
end

-- orientation function
function rtAttProto:orientation(tjd)
local t = tjd
-- return a quaternion
w, x, y, z = 1, 0, 0, 0
q = convert_quaternion ( w, x, y, z)
return q.w, -q.x, -q.y, -q.z
end

I think that the problem came from our cmod file. As a confirmation, we had to correct the mesh orientation by adding a line in ssc file.
The body axes were incorrect with the mesh definition. We added in ssc file :
Orientation [ 180 0 0.7017 0.7017 ]

Eric

[Fenerit]

It is often useful to be able to compute the location of one object projected onto the surface of another. This point on the surface is known as the subpoint. For example, we may want to know the latitude and longitude of the point directly beneath an Earth orbiting spacecraft. The following code defines a function called subpoint which is used to show the current latitude and longitude of the observer:

.
.
.

--Chris

This is very interesting. Reflecting on it I wonder whether would possible to do the same for the mouse pointer, so to have directly the Earth coordinates, if only numerically visible just close to it.

EDIT LATER

Never mind, it's a foolishly, what have in mind cannot be realized in this way; if I move the pointer toward the pole, how it could to know the subpoint? It wouldn't be perpendicular!

[Fenerit]

Well, I'll try to explain what I've in mind. Suppose that at start, at prompt, Celestia shows the US centered on Denver. The crosshair will be on Denver and by LUA EDU Tools I can see its coordinates. Now, with LUA EDU TOOLS I can see the coordinate for all places that are "centered" on screen, nor whether the crosshair is moved everywhere. Let's go back to example: once centered on Denver, I wonder whether I move the crosshair on Boston (without toggle the Earth) there would be a way to know also the Boston's coordinates and so on, for all place on Earth in which insist the crosshair; to say a way in which, by internal computation, Boston or whatever else location were as in center of screen (reference) and then, as visible output, showing it's
coordinates either within LUATOOLS, somewhere on screen or near the crosshair.

[duds26]

This are good functions but the subpoint is somewhat too general.
A better name would be something like: surfaceprojecttedpoint,
surfacepoint, bodiepoint,...

subpoint leaves, in my opinion, too much space open for other explanations, confusion.