Page 1 of 1

Gravity Script thrust issue

Posted: 17.04.2007, 22:28
by Hungry4info
I have downloaded the gravity script from the motherlode, it's quite nice, but doesn't seem to be completely functional. From what I undersand, pressing Ctrl + plus will provide foreward thrust, and Ctrl + minus provides backwards thrust. I've tried and such has never worked for me. I've even raised the thrust variable quite a bit, and yet no change. Here's what I've got: (my apologies, as it's somewhat large)

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?