Code: Select all
gKM_MLY = 9466411.842
gSEC_DAY = 86400
gMACH_KM_S = 0.34029
gG_KM_S2 = 0.0098 -- km/s^2
gTHRUST_ENERGY = 10000000.05 -- no units for now (burn out time not considered)
gBIG_NUMBER = 1e31
gSPHEROID_BOUND = 0.0001
gOBS_MASS = 85 -- kg
gOBS_RADIUS = 1 -- observer is sphere-shaped with this radius in meters
-- this is an oversimplification: drag coefficient depends on speed, air density, air viscosity, air temperature, etc,
-- all of which are strong functions of altitude. Also, each planet has different gaseous environments, so this should be
-- recomputed accordingly:
gDRAG_COEFF = 0.5
-- let's work in Celestia's own units: mLY^3 / (kg * day^2):
gG_CONST = 6.673e-20 * math.pow(gSEC_DAY,2) / math.pow(gKM_MLY,3)
gWAIT_DELAY = 0.0
gPRINT_DELAY = 15
gZERO = celestia:newposition(0, 0, 0)
gEARTH_MASS = 5.972e24 -- kg
gBodies = {
{name = "47 UMa", mass = 356666.31 },
{name = "47 UMa/b", mass = 826.8 },
{name = "47 Uma/c", mass = 426.12 },
{name = "55 Cnc", mass = 343333.3299 },
{name = "55 Cnc/b", mass = 267.12 },
{name = "55 Cnc/c", mass = 66.78 },
{name = "55 Cnc/d", mass = 1287.9 },
{name = "55 Cnc/e", mass = 14.31 },
{name = "Cygnus X-1/GCVS 311357",mass = 2500000. },
{name = "Gliese 876", mass = 106666.6656 },
{name = "Gliese 876/b", mass = 615.33 },
{name = "Gliese 876/c", mass = 178.08 },
{name = "Gliese 876/d", mass = 7.314 },
{name = "HD 37124", mass = 303333.03 },
{name = "HD 37124/b", mass = 193.98 },
{name = "HD 37124/c", mass = 217.194 },
{name = "HD 37124/d", mass = 190.8 },
{name = "HD 69830", mass = 286666.6638 },
{name = "HD 69830/b", mass = 10.2 },
{name = "HD 69830/c", mass = 11.8 },
{name = "HD 69830/d", mass = 18.1 },
{name = "Mu Ara", mass = 359999.64 },
{name = "Mu Ara/b", mass = 531.06 },
{name = "Mu Ara/c", mass = 985.8 },
{name = "Mu Ara/d", mass = 13.992 },
{name = "Mu Ara/e", mass = 575.58 },
{name = "Sol", mass = 333333.33 },
{name = "Sol/Mercury", mass = 0.055 },
{name = "Sol/Venus", mass = 0.82 },
{name = "Sol/Earth", mass = 1 },
{name = "Sol/Earth/Moon", mass = 0.0123 },
{name = "Sol/Mars", mass = 0.11 },
{name = "Sol/Jupiter", mass = 318 },
{name = "Sol/Saturn", mass = 95 },
{name = "Sol/Saturn/Mimas", mass = 0.000006363027 },
{name = "Sol/Saturn/Enceladus", mass = 0.000012223711 },
{name = "Sol/Saturn/Tethys", mass = 0.000104152713 },
{name = "Sol/Saturn/Dione", mass = 0.000175820496 },
{name = "Sol/Saturn/Rhea", mass = 0.004035498995 },
{name = "Sol/Saturn/Titan", mass = 0.022605492297 },
{name = "Sol/Saturn/Iapetus", mass = 0.000314802411 },
{name = "Sol/Uranus", mass = 14.5 },
{name = "Sol/Neptune", mass = 17.1 },
{name = "Sol/Pluto", mass = 0.0026 }
}
--***********************************************************
-- Selected Object
--***********************************************************
Selection = {body = nil, position = gZERO, speedVector = gZERO-gZERO, accelVector = gZERO-gZERO}
--***********************************************************
-- Point-like Observer Object
--***********************************************************
-- Also store some data just to avoid computing it more than once per cycle (ugly but faster):
MassPointObserver = {observer, position, speedVector, accelVector, t, dt,
frontArea, mass, dragConstant,
closestBody, closestBodyRadius, closestBodyDistance, collided
}
function MassPointObserver:getLatitude(relPos)
return math.deg(math.atan(math.sqrt(math.pow(relPos.x,2) + math.pow(relPos.y,2) ) / relPos.z))
end
function MassPointObserver:getLongitude(relPos)
return math.deg(math.atan(relPos.y/relPos.y))
end
function MassPointObserver:computeDistanceToSurface(distanceToBodyCenter)
-- returns observer's distance to body's surface in km:
-- first, compute observer's latitude relative to body:
local bodyFrame = celestia:newframe("Planetographic", self.closestBody)
local obsRelPos = bodyFrame:to(self.position)
local obsLatitude = self:getLatitude(obsRelPos)
self.closestBodyRadius = self.closestBody:radius() * math.sqrt(1 - math.pow(0*math.sin(obsLatitude),2))
self.closestBodyDistance = distanceToBodyCenter - self.closestBodyRadius
end
function MassPointObserver:setFinalPosition()
-- given a body-relative position, places the observer in the closest acceptable point to the body's surface
local closestBodyPos = self.closestBody:getposition()
self.position = ((self.position-closestBodyPos):normalize()*
(gSPHEROID_BOUND+1)*(self.closestBodyRadius/gKM_MLY) + closestBodyPos)
end
function MassPointObserver:applyFrontThrust()
self.speedVector = self.speedVector + gTHRUST_ENERGY * self:getLookingDirection()
end
function MassPointObserver:applyBackThrust()
self.speedVector = self.speedVector - gTHRUST_ENERGY * self:getLookingDirection()
end
function MassPointObserver:getLookingDirection()
return self.observer:getorientation():transform(celestia:newvector(0, 0, -1))
end
function MassPointObserver:initializeOriginalConditions()
-- determines original position, speed vector, etc.:
self.observer = celestia:getobserver()
self.crossArea = math.pi * math.pow(gOBS_RADIUS/(1000 * gKM_MLY), 2)
self.mass = gOBS_MASS
self.dragConstant = 0.5 * self.crossArea * gDRAG_COEFF / self.mass
self.collided = false
local t0 = celestia:gettime()
local p0 = self.observer:getposition()
wait(gWAIT_DELAY)
self.t = celestia:gettime()
self.dt = self.t - t0
self.position = self.observer:getposition()
self.speedVector = (self.position-p0) * (1 / self.dt)
self.accelVector = gZERO-gZERO
end
function MassPointObserver:solve()
-- solves the ODE using a 4th order Runge-Kutta method (for numerical stability at high time speeds):
self.dt = celestia:gettime() - self.t
local dv1 = externalAccel(self.t, self.position, self.speedVector) * self.dt
local dp1 = self.speedVector * self.dt
local dv2 = externalAccel(self.t + self.dt*0.5, self.position + dp1*0.5, self.speedVector + dv1*0.5) * self.dt
local dp2 = (self.speedVector + dv1*0.5) * self.dt
local dv3 = externalAccel(self.t + self.dt*0.5, self.position + dp2*0.5, self.speedVector + dv2*0.5) * self.dt
local dp3 = (self.speedVector + dv2*0.5) * self.dt
local dv4 = externalAccel(self.t + self.dt, self.position + dp3, self.speedVector + dv3) * self.dt
local dp4 = (self.speedVector + dv3) * self.dt
self.t = self.t + self.dt
self.position = self.position + (1/6) * (dp1 + 2*dp2 + 2*dp3 + dp4)
self.speedVector = self.speedVector + (1/6) * (dv1 + 2*dv2 + 2*dv3 + dv4) + self.accelVector * self.dt
self.accelVector = externalAccel(self.t, self.position, self.speedVector) + self.accelVector
end
function MassPointObserver:checkForCollision()
-- this won't work with arbitrary shaped bodies. (Access to mesh vertex data is needed for that):
if self.closestBodyDistance <= gSPHEROID_BOUND * self.closestBodyRadius then
self:setFinalPosition()
self.collided = true
end
end
function MassPointObserver:hasCollided()
return self.collided
end
function MassPointObserver:updatePosition()
self.observer:setposition(self.position)
end
function MassPointObserver:updateSelection()
local sel = celestia:getselection()
if sel:type() ~= "null" then
-- these values are just for the screen dump, so no refinement is really needed:
Selection.body = sel
local newPosition = sel:getposition(self.t)
local newSpeedVector = (newPosition - Selection.position) * (1/self.dt)
local newAccelVector = (newSpeedVector - Selection.speedVector) * (1/self.dt)
Selection.position = newPosition
Selection.speedVector = newSpeedVector
Selection.accelVector = newAccelVector
else
Selection.body = nil
end
end
function MassPointObserver:printState()
C1 = gKM_MLY / gSEC_DAY
C2 = gKM_MLY / (gSEC_DAY * gSEC_DAY)
local currSpeed = self.speedVector:length() * C1
local currAccel = self.accelVector:length() * C2
local selSpeed = (Selection.speedVector-self.speedVector):length() * C1
local selAccel = (Selection.speedVector-self.accelVector):length() * C2
local str = string.format("[Relative to UCS]\n Speed: %.4f km/s (%.3f Mach)\n Acceleration: %.7f km/s?? (%.3f G)",
currSpeed, currSpeed/gMACH_KM_S, currAccel, currAccel/gG_KM_S2)
if Selection.body ~= nil then
str = string.format("%s\n\n[Relative to selection]\n Speed: %.4f km/s (%.3f Mach)\n Acceleration: %.4f km/s?? (%.3f G)",
str, selSpeed, selSpeed/gMACH_KM_S, selAccel, selAccel/gG_KM_S2)
end
celestia:print(str, 1, -1, 0)
end
function MassPointObserver:resetForces()
self.accelVector = gZERO-gZERO
end
function MassPointObserver:stop()
-- sets observer's speed to 0 (for builtin speeds that are set with the keyboard function keys):
self.observer:setspeed(0)
end
--***********************************************************
-- Other functions
--***********************************************************
function externalAccel(tn, pn, sn)
-- the external forces accumulation function (this should be an internal method):
local minimumChord = gBIG_NUMBER
local accelVectorAC = gZERO-gZERO
for k, v in pairs(gBodies) do
local body = celestia:find(v.name)
local chord = body:getposition(tn) - pn
local chordLen = chord:length()
-- determine the closest body:
if chordLen < minimumChord then
minimumChord = chordLen
MassPointObserver.closestBody = body
end
accelVectorAC = accelVectorAC + ( v.mass * gEARTH_MASS / math.pow(chordLen, 3) ) * chord
end
accelVectorAC = gG_CONST * accelVectorAC
-- determine if observer is inside the atmosphere. If so, compute drag component:
MassPointObserver:computeDistanceToSurface(minimumChord * gKM_MLY)
local atmosphereHeight = MassPointObserver.closestBody:getinfo().atmosphereHeight
if atmosphereHeight ~= nil then
if MassPointObserver.closestBodyDistance <= atmosphereHeight then
--TODO: compute drag component more or less "seriously" (ie. consider at least a piecewise Reynolds "gradient", etc.)
end
end
return accelVectorAC
end
--***********************************************************
-- Callbacks
--***********************************************************
function celestia_keyboard_callback(input)
if input == "\044" then -- ctrl-plus
MassPointObserver:applyFrontThrust()
return true
elseif input == "\046" then -- ctrl-minus
MassPointObserver:applyBackThrust()
return true
else
return false
end
end
--***********************************************************
-- Main
--***********************************************************
celestia:requestkeyboard(true)
MassPointObserver:initializeOriginalConditions()
while not(MassPointObserver:hasCollided()) do
MassPointObserver:solve()
MassPointObserver:checkForCollision()
MassPointObserver:updatePosition()
MassPointObserver:updateSelection()
MassPointObserver:printState()
MassPointObserver:resetForces()
wait(gWAIT_DELAY)
end
MassPointObserver:stop()
celestia:print("TouchDown! End of script", gPRINT_DELAY)
Does anybody know how to get thrust to work?