(RA in hour/min/sec, Dec in deg/arc min/arc sec)
Code: Select all
-- Title: Display current RA/Dec for observer
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 observer from Earth:
getR = function (obs)
return earth:getposition():distanceto(obs:getposition())
end
-- Return current RA, Dec for observer:
getRADec = function (obs)
local base_rot = celestia:newrotation(celestia:newvector(1, 0, 0), -J2000Obliquity)
local rot = earth:getposition():orientationto(obs: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
obs = celestia:getobserver()
obsR = getR(obs)
obsRA, obsDec = getRADec(obs)
celestia:print(obsR..' '..obsRA..' '..obsDec)
if obsR >= 0 then
-- Display geocentric coordinates for observer:
obsRStr = km2Unit(obsR)
obsRAStr = deg2hms(obsRA)
obsDecStr = deg2dms(obsDec)
celestia:print("Geocentric coordinates for observer:\nR: "..obsRStr.."\nRA: "..obsRAStr.."\nDec: "..obsDecStr, 1, -1, -1, 1, 7)
end
wait(0)
end
