I have been working on a simulation and have got a spaceship in orbit around a planet. Now I want to present some information useful for navigation, in particular for landing on the surface. So the speed the surface moving under you is important, but for the life of me I can not figure it out or find the right information on the web (or the math forum)
I have two hunches. The first is that I need to find the tangential plane to the surface of the planet and project the velocity vector onto this plane, which will match the vector you find in this plane for the velocity of the rotation of the body. This is not a simple orthogonal projection though, so there I get stuck.
My second hunch is that there's a chance converting the vector represention to Kepler elements might be better, scaling the distance between two moments in time to the circle the ship orbits, but perhaps this is only easier in the idealised circular, non inclined orbit. (Isn't everything?)
Some details:
This simulation is running on a 1mhz Commodore 64, so it's highly likely the code to calc the speed, will take as long as updating the physics, I might have to forego considering the rotation speed, but basically the point is I am looking for a method that is both computationally cheap, and memory friendly. Approximations are okay, because of the resoloution of the dimensions, you'll be landing at unrealistic speeds anyway. I have slow and accurate trig routines, and fast and approximate ones, so with lots of fast trig the errors compound.
Any pointers are appreciated. Nice to be on this forum, and hope to have anyone interested flying an 8-bit spaceship soon.
Hi, 'malcontent'
I've been following this thread off and on for the last few weeks. Here's an approach to your problem that you may want to consider.
As I understand it, you have developed an integration engine that reproduces a normal Keplerian (i.e., circular, elliptical, parabolic or hyperbolic orbit) around a planet. For the most part, these integration schemes tend to be set up in an inertial reference frame. For many problems, the choice of an inertial reference frame makes ensuing calculations simpler. However, for some problems - and possibly for your problem - using an integration scheme in inertial space makes life harder.
It is always possible, though, to set up a numerical integration scheme in a reference frame that co-rotates with the planet such that all points on the planet have a fixed position. In the past (and largely form my own amusement), I've developed a number of explicit symplectic integration schemes that work in these rotating reference frames. Since, the position of all points on the surface of the planet are fixed, determining the speed of your vessel with respect to a (now stationary) point directly underneath your vessel can be calculated using straightforward projective geometry. At worst, this calculation involves an inverse square root calculation which is something that you can already calculate reasonably cheaply.
As far as the symplectic integrator is concerned only simple arithmetic calculations are required (i.e., no sine or cosine evaluations); and can be constructed as a second order, fourth order or higher order integrator as desired.
If this is of interest to you, I'll post the second-order symplectic integrator (in 3D, rotating coordinates) so that you can play around with it and see if it meets your requirements. If you need higher order integration schemes, I can show you how to combine the second-order scheme to construct higher-order schemes.
P.S. If you also require the integrated trajectory in inertial coordinates (as well as rotating coordinates) then rather convert between inertial and rotating reference frames, the easier thing to do may be to perform two parallel integrations - one in inertial coordinates; and one in rotating coordinates.
---------- Post added 02-28-17 at 03:03 AM ---------- Previous post was 02-27-17 at 06:54 AM ----------
For anyone interested, here is a Lua script for a numerical integration of an orbit around a planet (in this case, the Earth) in rotating coordinates and calculating the grounds speed of the vessel without using sines and cosines.
The inputs to the integration are the position and the velocity of the vessel in rotating coordinates. In the specific example given below, the vessel is in a 300 x 300 km polar orbit starting from the equator and heading North towards the pole. In inertial coordinates, this means that it should be orbiting at around 7725.839 m/s, but in rotating coordinates it also picks up a perpendicular component due to the rotation of the Earth of -486.967 m/s.
The integration engine here is a simple second-order (symplectic) integrator. Higher order integrators can easily be constructed from this using the scheme given here
http://www.orbiter-forum.com/showthread.php?t=35489&highlight=symplectic
The integrator requires one inverse square root calculation per time-step; and to calculate the ground speed, one further inverse square root calculation is required and one sort calculation. All other operations are simple arithmetic operations +, -, * , / . No sines or cosines are required.
Although written in Lua, it should be easy to convert to other programming languages.
It is should also be easy to calculate ground tracks using the same algorithm,
Code:
-- enter a few specific details for the planet, in this case 'the Earth'
GM = 3.986004418e14 -- standard gravitational parameter
R = 6378000.0 -- mean radius of the Earth
w = 2 * math.pi / (86400 * 0.99726958) -- Earth's sidereal rate of rotation
-- enter some numerical integration constants
dt = 10 -- integration time-step
phi = w * dt -- the product of rate of rotation and the time step
-- define the second-order symplectic integrator - for rotating coordinates
integrator = function (Q0, V0)
-- initialise some quantities used for internal calculation
local iphi = 0
local rho = 0
local rho3 = 0
local irho3= 0
local Qa = {x = 0, y = 0, z = 0}
local Q1 = {x = 0, y = 0, z = 0}
local V1 = {x = 0, y = 0, z = 0}
local F = {x = 0, y = 0, z = 0}
-- the first part of the leap-frog integration scheme
Qa.x = Q0.x + 0.5 * dt * V0.x
Qa.y = Q0.y + 0.5 * dt * V0.y
Qa.z = Q0.z + 0.5 * dt * V0.z
-- the second part of the leap-frog integration scheme
rho = 1 / math.sqrt( Qa.x * Qa.x + Qa.y * Qa.y + Qa.z * Qa.z)
irho3= rho * rho * rho
iphi = 1 / (1 + phi * phi)
F.x = - GM * Qa.x * irho3 + w * w * Qa.x -- inverse square law and centrifugal force
F.y = - GM * Qa.y * irho3 + w * w * Qa.y -- inverse square law and centrifugal force
F.z = - GM * Qa.z * irho3 -- inverse square law only
V1.x = ( (1 - phi * phi) * V0.x + 2 * phi * V0.y + dt * (F.x + phi * F.y)) * iphi
V1.y = (-2 * phi * V0.x + (1 - phi * phi) * V0.y + dt * (F.y - phi * F.x)) * iphi
V1.z = V0.z + dt * F.z
-- the third part of the leap-frog integration scheme
Q1.x = Qa.x + 0.5 * dt * V1.x
Q1.y = Qa.y + 0.5 * dt * V1.y
Q1.z = Qa.z + 0.5 * dt * V1.z
-- return the result of the integration step
return Q1, V1
end
groundSpeed = function (Q, V)
local gspd = 0.0
local rhat = {x = 0.0, y = 0.0, z = 0.0}
local Vw = {x = 0.0, y = 0.0, z = 0.0}
local Vrhat = 0.0
local inormQ = 0.0
inormQ = 1 / math.sqrt(Q.x * Q.x + Q.y * Q.y + Q.z * Q.z)
rhat.x = Q.x * inormQ
rhat.y = Q.y * inormQ
rhat.z = Q.z * inormQ
-- calculate the projection of the velocity vector onto the Earth's surface
Vrhat = V.x * rhat.x + V.y * rhat.y + V.z * rhat.z
Vw.x = (V.x - rhat.x * Vrhat) * inormQ * R
Vw.y = (V.y - rhat.y * Vrhat) * inormQ * R
Vw.z = (V.z - rhat.z * Vrhat) * inormQ * R
gspd = math.sqrt(Vw.x * Vw.x + Vw.y * Vw.y + Vw.z * Vw.z)
return gspd
end
-- set up the initial conditions for position and velocity in rotating coordinates
Q = {x = 6678000.0, y = 0.0, z = 0.0}
V = {x = 0.0, y = -486.96749, z = 7725.83948}
print( groundSpeed (Q, V))
for idx = 1, 400 do
Q, V = integrator(Q, V)
print( groundSpeed (Q, V))
end