-- Phase test
-- Test the new timeline features of Celestia

km = 9460730.4725808


sel = celestia:getselection()

function date_to_string(tdb)
    if tdb == math.huge then
	return "+Inf"
    elseif tdb == -math.huge then
	return "-Inf"
    else
       t = celestia:tdbtoutc(tdb)
       return string.format("%02d:%02d:%02d %d %d %d", t.hour, t.minute, t.seconds, t.year, t.month, t.day)
    end
end


function position_to_string(p)
   return string.format("%f %f %f", p.x * km, -p.z * km, p.y * km)
end


function rotation_to_euler(q)
   local qx = q.x
   local qy = q.y
   local qz = q.z
   local qw = q.w

   local Nq = qw * qw + qx * qx + qy * qy + qz * qz
   local s = 2.0 / Nq
   
   local xs = qx * s
   local ys = qy * s
   local zs = qz * s

   local wx = qw * xs
   local wy = qw * ys
   local wz = qw * zs
   local xx = qx * xs
   local xy = qx * ys
   local xz = qx * zs
   local yy = qy * ys
   local yz = qy * zs
   local zz = qz * zs

   local m00 = 1 - (yy + zz)
   local m11 = 1 - (xx + zz)
   local m22 = 1 - (xx + yy)

   local m01 = xy - wz
   local m10 = xy + wz
   local m02 = xz + wy
   local m20 = xz - wy
   local m12 = yz - wx
   local m21 = yz + wx

--[[
   -- XYX
   local sy = math.sqrt(m01*m01+m02*m02)
   if (sy > 1e-10) then
      ex = math.atan2(m01, m02)
      ey = math.atan2(sy, m00)
      ez = math.atan2(m10, -m20)
   else
      ex = math.atan2(-m12, m11)
      ey = math.atan2(sy, m00)
      ez = 0
   end
]]


--[[
   -- YZY
   local sy = math.sqrt(m12*m12+m10*m10)
   if (sy > 1e-10) then
      ex = math.atan2(m12, m10)
      ey = math.atan2(sy, m11)
      ez = math.atan2(m21, -m01)
   else
      ex = math.atan2(-m20, m22)
      ey = math.atan2(sy, m11)
      ez = 0
   end
]]


--[[
   -- ZXZ
   local sy = math.sqrt(m20*m20+m21*m21)
   if (sy > 1e-10) then
      ex = math.atan2(m20, m21)
      ey = math.atan2(sy, m22)
      ez = math.atan2(m02, -m12)
   else
      ex = math.atan2(-m01, m00)
      ey = math.atan2(sy, m22)
      ez = 0
   end
]]


   -- YXY
   local sy = math.sqrt(m10*m10+m12*m12)
   if (sy > 1e-10) then
      ex = math.atan2(m10, m12)
      ey = math.atan2(sy, m11)
      ez = math.atan2(m01, -m21)
   else
      ex = math.atan2(-m02, m00)
      ey = math.atan2(sy, m11)
      ez = 0
   end

   return ex, ey, ez
end


function weld(p0, p1, n)
   -- t == p0.ending == p1.ending
   t = p1:timespan()

   pos_univ = p1:orbitframe():from(p1:getposition(t), t)
   pos_p0 = p0:orbitframe():to(pos_univ, t)
   rot_univ = p1:bodyframe():from(p1:getorientation(t), t)
   rot_p0 = p0:bodyframe():to(rot_univ, t)
   
   pos_univ = p0:orbitframe():from(p0:getposition(t), t)
   pos_p1 = p1:orbitframe():to(pos_univ, t)
   rot_univ = p0:bodyframe():from(p0:getorientation(t), t)
   rot_p1 = p1:bodyframe():to(rot_univ, t)

   local node, inc, meridian

   celestia:log(date_to_string(t))

   celestia:log(string.format("Phase %d in %d: ", n, n + 1))
   celestia:log(position_to_string(pos_p1))
   node, inc, meridian = rotation_to_euler(rot_p1)
   celestia:log(string.format("AscendingNode %f  Inclination %f  MeridianAngle %f", math.deg(node), math.deg(inc), math.deg(meridian)))

   celestia:log(string.format("Phase %d in %d: ", n + 1, n))
   celestia:log(position_to_string(pos_p0))
   node, inc, meridian = rotation_to_euler(rot_p0)
   celestia:log(string.format("AscendingNode %f  Inclination %f  MeridianAngle %f", math.deg(node), math.deg(inc), math.deg(meridian)))
end


lastphase = nil
count = 0
for p in sel:phases() do
   if lastphase then
      weld(lastphase, p, count)
   end

   lastphase = p
   count = count + 1
end

rot = celestia:newrotation(1, 0, 0, 0)
rot:setaxisangle(celestia:newvector(1, 0, 0), -math.pi / 2)

t = celestia:gettime()
q = lastphase:getorientation(t)
-- q = rot * q
ex, ey, ez = rotation_to_euler(q)

celestia:log(string.format("%f %f %f", math.deg(ex), math.deg(ey), math.deg(ez)))