Script to compare coordinate

[rinoa79]

Hi folk!
I have created a script to compare ISS's real coordinates (RA, Dec, r) found in this site:
http://www.n2yo.com/
with ISS's coordinates calculated using Celestia.
But they didn't correspond.
Please, can someone tell me why?
I think I made some mistakes but I don't understand where!

This is my code:

Code: Select all

earth = celestia:find("Sol/Earth")
pos_earth = terra:getposition()
sist_geo_cart = celestia:newframe("equatorial", earth)

ISS = celestia:find("Sol/Earth/ISS")
pos_ISS = ISS:getposition()

-- ISS's universal coordinates
coordinataX = tostring(pos_ISS.x)
coordinataY = tostring(pos_ISS.y)
coordinataZ = tostring(pos_ISS.z)

-- ISS's cartesian equatorial coordinates in km
pos_ISS_geoc = sist_geo_cart:to(pos_ISS)
coord_geoc_X = tostring(pos_ISS_geoc.x * 9466411.842)
coord_geoc_Y = tostring(pos_ISS_geoc.y * 9466411.842)
coord_geoc_Z = tostring(pos_ISS_geoc.z * 9466411.842)
   
-- coordinate (alpha, delta, r) in deg
RA = math.deg(math.atan(coord_geoc_Z /coord_geoc_X))
DEC = math.deg(math.atan(coord_geoc_Y /(math.sqrt(coord_geoc_X ^2 + coord_geoc_Z ^2))))
r = math.sqrt(coord_geoc_X ^2+coord_geoc_Y ^2+coord_geoc_Z ^2)

celestia:print(RA.." "..DEC.." "..r,30)
wait(30)

Thanks in advance,
Mary.

[selden]

Celestia's orbit for the ISS is not kept up to date. It changes too often. This is explained in the FAQ.

[Vincent]

Mary,

On top of the issue raised by Selden, your script does not display the right Celestia values for R, Ra and Dec.
Though I don't have the time to investigate what is wrong with your script, you can use the following script that
displays R, RA, Dec for the selected object. RA/Dec are displayed respectively in hour/min/sec and in deg/min/arcsec.
You should be able to easily edit the script to make them displayed in decimal degrees.

Code: Select all

-- Title: Display current RA/Dec for selected object

KM_PER_LY = 9460730472580.8
KM_PER_AU = 149597870.7
PI = math.pi
degToRad = PI / 180;
J2000Obliquity = 23.4392911 * degToRad

LOOK  = celestia:newvector(0, 0, -1)
earth = celestia:find("Sol/Earth")

-- Convert coordinates from cartesian to polar:
xyz2rtp = function (x, y, z)
   local r = math.sqrt(x * x + y * y + z * z)
   local phi = math.atan2(y, x)
   local theta = math.atan2(math.sqrt(x * x + y * y), z)

   return r, theta, phi
end

-- Return current distance of object from Earth:
getR = function (obj)
   return earth:getposition():distanceto(obj:getposition())
end

-- Return current RA, Dec for object:
getRADec = function (obj)
   local base_rot = celestia:newrotation(celestia:newvector(1, 0, 0), -J2000Obliquity)
   local rot = earth:getposition():orientationto(obj:getposition(), LOOK) * base_rot
   local look = rot:transform(LOOK):normalize()
   local r, theta, phi = xyz2rtp(look.x, look.z, look.y)
   local phi = math.mod(720 - math.deg(phi), 360)
   local theta = math.deg(theta)
   if theta > 0 then
      theta = 90 - theta
   else
      theta = (-90 - theta)
   end
   return phi, theta
end

km2Unit =
   function(km)
      local sign, value, units
      if km < 0 then sign = -1 else sign = 1 end
      km = math.abs(km)
      if km > 1e12 then
         value = km / KM_PER_LY
         units = "ly"
      elseif km >= 1e8 then
         value = km/KM_PER_AU
         units = "AU"
      else
         value = km
         units = "km"
      end;
      return string.format("%.2f", sign * value).." "..units
   end

deg2dms = function(deg)
   local   a = math.abs(deg)
  local   d = math.floor(a)
   local   r = (a - d) * 60
   local   m = math.floor(r)
   local   s = (r - m) * 60
   if deg < 0 then d = -d end
   return string.format("%0.0fd %02.0f' %2.0f''",d,m,s)
end

deg2hms =   function(deg)
   local   a = math.abs(deg / 15)
  local   d = math.floor(a)
   local   r = (a - d) * 60
   local   m = math.floor(r)
   local   s = (r - m) * 60
   return string.format("%0.0fh %02.0fm %2.0fs", d, m, s)
end

while true do
   obj = celestia:getselection()
   objR = getR(obj)
   objRA, objDec = getRADec(obj)
   if objR >= 0 then
      -- Display geocentric coordinates for object:
      objRStr = km2Unit(objR)
      objRAStr = deg2hms(objRA)
      objDecStr = deg2dms(objDec)
      celestia:print("Geocentric coordinates for "..obj:name()..":\nR: "..objRStr.."\nRA: "..objRAStr.."\nDec: "..objDecStr, 1, -1, -1, 1, 7)
   end
   wait(0)
end

[rinoa79]

Mary,

On top of the issue raised by Selden, your script does not display the right Celestia values for R, Ra and Dec.
Though I don't have the time to investigate what is wrong with your script, you can use the following script that
displays R, RA, Dec for the selected object. RA/Dec are displayed respectively in hour/min/sec and in deg/min/arcsec.
You should be able to easily edit the script to make them displayed in decimal degrees.

Code: Select all

-- Title: Display current RA/Dec for selected object

KM_PER_LY = 9460730472580.8
KM_PER_AU = 149597870.7
PI = math.pi
degToRad = PI / 180;
J2000Obliquity = 23.4392911 * degToRad

LOOK  = celestia:newvector(0, 0, -1)
earth = celestia:find("Sol/Earth")

-- Convert coordinates from cartesian to polar:
xyz2rtp = function (x, y, z)
   local r = math.sqrt(x * x + y * y + z * z)
   local phi = math.atan2(y, x)
   local theta = math.atan2(math.sqrt(x * x + y * y), z)

   return r, theta, phi
end

-- Return current distance of object from Earth:
getR = function (obj)
   return earth:getposition():distanceto(obj:getposition())
end

-- Return current RA, Dec for object:
getRADec = function (obj)
   local base_rot = celestia:newrotation(celestia:newvector(1, 0, 0), -J2000Obliquity)
   local rot = earth:getposition():orientationto(obj:getposition(), LOOK) * base_rot
   local look = rot:transform(LOOK):normalize()
   local r, theta, phi = xyz2rtp(look.x, look.z, look.y)
   local phi = math.mod(720 - math.deg(phi), 360)
   local theta = math.deg(theta)
   if theta > 0 then
      theta = 90 - theta
   else
      theta = (-90 - theta)
   end
   return phi, theta
end

km2Unit =
   function(km)
      local sign, value, units
      if km < 0 then sign = -1 else sign = 1 end
      km = math.abs(km)
      if km > 1e12 then
         value = km / KM_PER_LY
         units = "ly"
      elseif km >= 1e8 then
         value = km/KM_PER_AU
         units = "AU"
      else
         value = km
         units = "km"
      end;
      return string.format("%.2f", sign * value).." "..units
   end

deg2dms = function(deg)
   local   a = math.abs(deg)
  local   d = math.floor(a)
   local   r = (a - d) * 60
   local   m = math.floor(r)
   local   s = (r - m) * 60
   if deg < 0 then d = -d end
   return string.format("%0.0fd %02.0f' %2.0f''",d,m,s)
end

deg2hms =   function(deg)
   local   a = math.abs(deg / 15)
  local   d = math.floor(a)
   local   r = (a - d) * 60
   local   m = math.floor(r)
   local   s = (r - m) * 60
   return string.format("%0.0fh %02.0fm %2.0fs", d, m, s)
end

while true do
   obj = celestia:getselection()
   objR = getR(obj)
   objRA, objDec = getRADec(obj)
   if objR >= 0 then
      -- Display geocentric coordinates for object:
      objRStr = km2Unit(objR)
      objRAStr = deg2hms(objRA)
      objDecStr = deg2dms(objDec)
      celestia:print("Geocentric coordinates for "..obj:name()..":\nR: "..objRStr.."\nRA: "..objRAStr.."\nDec: "..objDecStr, 1, -1, -1, 1, 7)
   end
   wait(0)
end

Hi Vincent,
I have tested your script.
I have choosen a satellite: SEASAT 1 and have created a file .ssc with its orbital parameters.
I have calculated RA, DEC and I have compared them with coordinates calculate using ORBITRON like reference.
But they don't corrispond...
Mary.

[selden]

Mary,

What values do you get?

[rinoa79]

Mary,

What values do you get?

RA = 13h 33min...
DEC = -57° 35' ...
but real values are around:
RA = 3h...
DEC = -9°...

[rinoa79]

I forgot to tell you that I have changed only the first line with this:

Code: Select all

obj = celestia:find("Sol/Earth/SEASAT1")

I maybe made a mess...

[Vincent]

Mary,

Actually, my script only gives accurate RA/Dec values for distant objects (stars, DSO's) and
when the observer is in the neighbourhood of Earth. If you want to get RA/Dec for spacecraft,
please replace the getRADec function with the following one that uses an observer centered
coordinate system. Then, the displayed values should match the Orbitron values, and also the values
obtained from the projection onto the Celestial grid.

Code: Select all

-- Return current RA, Dec for object:
getRADec = function (obj)
   obs = celestia:getobserver()
   local base_rot = celestia:newrotation(celestia:newvector(1, 0, 0), -J2000Obliquity)
   local rot = obs:getposition():orientationto(obj:getposition(), LOOK) * base_rot
   local look = rot:transform(LOOK):normalize()
   local r, theta, phi = xyz2rtp(look.x, look.z, look.y)
   local phi = math.mod(720 - math.deg(phi), 360)
   local theta = math.deg(theta)
   if theta > 0 then
      theta = 90 - theta
   else
      theta = (-90 - theta)
   end
   return phi, theta
end

[rinoa79]

Mary,

Actually, my script only gives accurate RA/Dec values for distant objects (stars, DSO's) and
when the observer is in the neighbourhood of Earth. If you want to get RA/Dec for spacecraft,
please replace the getRADec function with the following one that uses an observer centered
coordinate system. Then, the displayed values should match the Orbitron values, and also the values
obtained from the projection onto the Celestial grid.

Code: Select all

-- Return current RA, Dec for object:
getRADec = function (obj)
   obs = celestia:getobserver()
   local base_rot = celestia:newrotation(celestia:newvector(1, 0, 0), -J2000Obliquity)
   local rot = obs:getposition():orientationto(obj:getposition(), LOOK) * base_rot
   local look = rot:transform(LOOK):normalize()
   local r, theta, phi = xyz2rtp(look.x, look.z, look.y)
   local phi = math.mod(720 - math.deg(phi), 360)
   local theta = math.deg(theta)
   if theta > 0 then
      theta = 90 - theta
   else
      theta = (-90 - theta)
   end
   return phi, theta
end

Hi Vincent,
I tried your scripts with these parameters that I think they are correct:

"Hubble" "Sol/Earth"
{
EllipticalOrbit
{
Period 0.0666458
SemiMajorAxis 6936.51
Eccentricity 0.0003677
Inclination 28.4703
AscendingNode 240.866
ArgOfPericenter 245.866
MeanAnomaly 114.3777
Epoch "2008 10 18 4:51:45"
}
}

The original code gives RA enough correct (1h difference) but DEC is completly wrong. The values must be negative (in this instant) and in your code is positive.
When I have changed the piece of code how you told me, RA is wrong too...

[rinoa79]

Now original script gives me RA wrong too...