This is what I've managed to work out today:

__Thrust__
*T = r * (v_max - v) / v_max * t * e*
where:

**T** is the resulting thrust
**r** is atmospheric density (normalised to a max value of 1 at sea level)
**v** is the current velocity
**v_max** is the maximum possible velocity (1,000 m/s)
**t** is the throttle position (0..1)
**e** is engine's maximum output (120,000 N)

__Weight__
*W = m * g*
where:

**W** is the resulting weight
**m** is the current mass of the aircraft including fuel and weapons
**g** is 10 (close enough)

__Drag__
*D = Cd * r * v^2*
where:

**D** is the resulting drag
**Cd** is a constant
**r** is atmospheric density (normalised to a max value of 1 at sea level)
**v** is the current velocity

__Lift__
*L = Cl * r * v^2*
where:

**L** is the resulting lift
**Cl** is a constant
**r** is atmospheric density (normalised to a max value of 1 at sea level)
**v** is the current velocity

I need to calculate

**Cl** and

**Cd** from the other information I have. I don't know the correct way to do this, but I've thrown some calculations together based on a take off speed of 50 m/s and a max speed at sea level of 500 m/s and get

**Cl = 40** and

**Cd = 0.36**.

I have no idea if those numbers are right or even sensible, but it looks like they're my starting point.

I'm also very dubious about my thrust formula.

These formulae are all forces, so they all need to be divided by mass to get acceleration, which gets be a change in velocity per second. I will need to scale this value according to how long it's been since the last update.

My world coordinate system is in metres, shifted to the left by 8 bits, which gives a resolution of around 4mm. I don't know if this will give noticeable errors due to rounding, particularly at low velocity or acceleration. I may need to accumulate errors.