Haven't gotten to go through the answers yet, but to be clear there's only ever 1 gravity source per scene. I like the idea of using iterators and will probably try that approach over the weekend and see how it performs

So your gravity-style force is actually proportional to the inverse [not-squared] distance. Neither of my differential-equation experts had thought about this situation before.

I managed to find an energy function, which is something like E = 0.5*velocity^2 - log(distance*5000). Could you see if that quantity is preserved (up to integration error) in your trajectories, so I can see if I am on the right track?

So your gravity-style force is actually proportional to the inverse [not-squared] distance.

Doh, totally missed that...

So your gravity-style force is actually proportional to the inverse [not-squared] distance.

Doh, totally missed that...

But it IS proportional to the inverse squared distance, like real gravity:
newForce = (newForce / newForce.sqrMagnitude) * 5000.0f;

Or is it me missing something? Just make sure before spending time...

Notice he computes force via distance / |distance|^2, which boils down to 1/|distance| * normalized direction. Square laws use distance / |distance|^3 or equivalently 1/distance^2 * normalized direction.

Ah damn, ... not the first time that mistake happened to me

Some time ago Alvaro helped me to solve a similar question: http://www.gamedev.net/topic/672589-equations-of-motion-for-position-dependent-acceleration/
Trying to apply this practice here i get this formula for acceleration:

acc = 5000 / -(pos - gravityCenter)

Now giving this to Wolfram Alpha i get this: http://www.wolframalpha.com/input/?i=x%27%27+-+k%2F-x+%3D+0
It shows a formula to get position at a given time, nice!
This could be used to speed up a searching approach, or better, maybe it can be rearranged to get the time for a questioned target position.

The problem is i don't know what is 'inverse error function'?
My math background is too weak.

The differential equation that JoeJ plugged into Wolfram Alpha describes how an object would fall into the star, if its velocity didn't have a horizontal component (i.e., if it's moving directly towards the star or away from it).

In inverse error function is the inverse of the error function. It is very closely related to probability of a getting number less than some threshold from a random variable that follows a normal distribution.

I don't think differential equations wouldn't be too useful here. The purpose of a differential is to model the rate of change to a higher order function at some point x or y. It can't really be used to find a predictable path alone, but it can be used to find the changes in velocity, acceleration, and position.

Most if not all movement equations in physics are differential equations...besides, changes in position integrated over time is the path.

y(t) = 1/2At^2 + Ct + k, is now a parabolic curve that computes the position of (x, y, or z) in space to a "ground place" being effected by a constant force. C and K are constants of anything. Which makes this equation an "ok" general solution for this linear equation. But not a perfect one.

If you pose your problem correctly, i.e. prescribe position and velocity, C and K are not "anything", but uniquely defined. This is then the only solution to this equation.

By the way, most if not all problems in mechanics and engineering can be reduced to solving (partial or ordinary) differential equations. You grossly understate their utility.

Of course OPs problem can be formulated as a ODE using Newton's law as a well-enough approximation:

our gravitational potentials are of the form V_i(r) = - G * M / |r - r_i| (central potential), so the force at point r can be written (if I'm not mistaken) as F_i(r) = grad V_i(r) = -G * M (r - r_i) / |r - r_i|^3.

Given an initial position r_0 and initial velocity v_0, you have enough data to compute the exact path. Plug this into m*a(t) = F(t) and you obtain an ODE for your path. I don't think solving this exactly is easy (maybe not even possible). You can certainly solve this numerically with ODE step solvers, but OP is looking at the inverse problem, which is hard in itself already. I don't know if this can be solved without some heuristic "aim, guess resulting path, adjust aim" algorithm, unless the ODE can be solved exactly. Maybe I'm off on this one, who knows

Issue still is that you won't be able to generate a full path ^-^; you can get the vector at any given point, but your graph will be for that point in space only. The forces of gravity is not constant. The forces of gravity are a function of your distance relative to each body of mass. Which does not create an ODE. But I think a partial differential equation might be able to solve it. But you will have more than two independents.

Your current location in the x and y plane, then all of the gravitational forces acting on you in relation to your distance to them. You're probably going to do better with looking into the future, than deriving a general equation for X number of celestial bodies.

Issue still is that you won't be able to generate a full path ^-^; you can get the vector at any given point, but your graph will be for that point in space only. The forces of gravity is not constant. The forces of gravity are a function of your distance relative to each body of mass. Which does not create an ODE. But I think a partial differential equation might be able to solve it. But you will have more than two independents.

This is simply not true. Not only are the equations of motion under a force that is described by a vector field and example of ODEs, but they are actually the first example of ODEs historically.

https://en.wikipedia.org/wiki/Equations_of_motion says:

This is equivalent to saying an equation of motion in r is a second order ordinary differential equation (ODE) in r,[...]