Calculate spaceship thrust intensity (only particle system) based on direction and velocity

I’I have a spaceship which can lead to any place of the 3D space and this ship has a thrust.

What I’d like to do is to adjust the thrust intensity based on the facing direction of the ship and its velocity. Please note that I’m talking about the visuals, that is the particle system of the ship, and not the physics itself. That part is done already. Velocity can be found in the root RigidBody and facing direction can be calculated.

The reason for specifying all of this, is that the ship could be moving in some direction because of external forces (e.g. black hole), and this is the velocity, but the velocity might not match the direction the ship is facing, its forward.
Example 1: black hole: velocity points to black hole, direction points in the opposite, because the ship is trying to escape it.
Example 2: space battle: a ship passes by an enemy ship, and the ship might need to change its rotation (slowly) with side thrusts to face the enemy and fire.

Now, I wrote a script which is considering only the Z axis (the forward, in fact) and then the plan was to proceed adding X and Y variables, but it quickly became clear that the script was too long, and I’m sure there’s a better way of doing this, maybe involving a better use of quaternion (tbh, I don’t know how they work).

This is the uncompleted script:

    public class EnergyBlast : MonoBehaviour
    {
        private ParticleSystem _particleSystem;
        private Rigidbody _unitRigidBody;
        private float _unitMaxSpeed;
      
        [SerializeField] private float maxStartLifetime = 0.2f;
        private void Start()
        {
            _particleSystem = ComponentUtils.RequireComponent<ParticleSystem>(this);
            _unitRigidBody = ComponentUtils.RequireComponent<Rigidbody>(transform.root);
            _unitMaxSpeed = ComponentUtils.RequireComponent<ShipAttributes>(transform.root).maxSpeed;
        }

        private void Update()
        {
            //Vector3 directionRotation = transform.root.transform.rotation.eulerAngles; // Equals _unitRigidBody.rotation.eulerAngles
            Vector3 unitRotation = _unitRigidBody.rotation.eulerAngles;
          
          
            // Ensure the direction vector is normalized
            Vector3 normalizedVelocity = _unitRigidBody.velocity.normalized;

            // Get the quaternion representing the rotation towards the direction vector
            Quaternion rotation = Quaternion.LookRotation(normalizedVelocity);

          
            // Convert the quaternion to Euler angles
            Vector3 velocityRotation = rotation.eulerAngles;

            var main = _particleSystem.main;
          
          
            var abs = Math.Abs(unitRotation.z - velocityRotation.z);
            if (abs < 90)
            {
                // Se abs = 0 --> max emissions
                // Se abs = 89.9 --> zero emissions
                var factor = 1 - abs/90;

                Vector3 normalizedirectionVector = transform.root.transform.rotation * Vector3.forward;
                float speedFactor = normalizedirectionVector.z - _unitRigidBody.velocity.z; // assume that max velocity z = 1, but it's not true
              
                // (Blocked)
                // normalizedirectionVector.z = 0 ; _unitRigidBody.velocity.z = 0, then
                // normalizedirectionVector.z = 1 ; _unitRigidBody.velocity.z = 0, then
                // normalizedirectionVector.z = 0 ; _unitRigidBody.velocity.z = 1, then
                // normalizedirectionVector.z = 1 ; _unitRigidBody.velocity.z = 1, then
              
                main.startLifetime = speedFactor * maxStartLifetime * factor;
            }
            else
            {
                main.startLifetime = 0;
            }

          
          
        }
      
    }

I hope that someone can help me, at least pointing me in the right direction.

As you can see, I was trying not only to consider the modules of the velocity and the direction, but also the angle.

Oh man what you ask is relatively simple but it’s a very mentally demanding thing to explain with words. It requires pictures and explanations about vector decomposition, and also I can’t do much from the top of my head, and I don’t really have much time at this moment.

If I understand you correctly, you want to model Newtonian physics, where the ship is free to rotate (via attitude thrusters), and you want to maintain the given thrust vector while the body is free to rotate, change this thrust vector, or have different configuration of thrusters, and is affected by external forces (impacts or continuous forces like gravity).

First of all let’s distinguish different components of the system. You have thrusting force and thrusting vector (producing acceleration or change in velocity), you have attitude thrusters which affect rotation, and you have external forces lined up at different angles. Some of them like impacts have a point of contact, some of them affect the body as a whole, like gravitation.

The biggest problem about these things begin when you realize you need torque r x F. For example torque happens if you want to treat your attitude thrusters as smaller variants of the main thrusters, or if your ship is hit by projectiles. This is because an off-center force will generate torque, which you want to neutralize unless you like your crew being splattered on the inside walls. What does off-center means? It means your ship has a special point called the center of mass, which is typically inside the hull (but only if your ship’s shape is convex). A coaxial thruster run directly through this center of mass, while attitude and off-center thrusters do not, but the most fun in such games comes from thruster arrays which have to be balanced.

You need to be able to decompose each thrusting force in such a way to easily tell how much translates into linear inertia, and how much into torque.

Of course if you use physics engine, these things should be automatic for you. However…

This part tells me that you want to decompose the thrust vectors so that you can apply forces to various thrusters as individual components of that thrust, such that ||thrust||^2 = Sum(||thrusters||^2)

Let’s think about a simple model, let’s imagine you have only a back thruster (capable of 10 m/s), and you want to thrust forward 5m/s. It’s easy to see that you need to fire your back thruster at 50%.

Now imagine what happens when you ship rotates by 30 degrees to the left. The back thruster vector is split at this angle, one component contributes to your thrust vector, the other turns into torque which is a waste (edit: or if it’s a coaxial thruster, no torque is generated but it sends you in the wrong direction, and so you want to fire off different thrusters to compensate). So you want to cancel this torque with another thruster which is detected as producing the best outcome, plus you also want to engage another thruster to help with the main thrust with the minimum of added torque. Typically, when designing a ship, you want all primary thrusters to not be off-center, and all your attitude thrusters to be extremely off-center.

These things get complicated very quickly if you’re approaching this naively. (edit: It’s in fact a very fun problem typically solved with explicit and implicit methods.)

Of course, maybe you want this only as a visual information, and so you only want to update thruster visual to reflect the motion of the ship, without actually generating forces, without considering torque etc. So to find out how much the thruster (visually) contributes to your current flight direction, you still want to decompose the force vector in such a way to have one of its perpendicular components perfectly aligned with the direction vector.

This is done by means of trigonometry and it’s easier to compute with a dot product, because it gives you a cosine between the angles. This in turn practically makes the thrust vector project onto the velocity line, which you can find with Vector3.Project.

A dot product will compare two vectors and spit out a value between -1 and 1, which is great for you, because 1 means they point at the same direction, and -1 that they’re directly opposite to each other. 0 means that they’re perfectly perpendicular to each other. So you can very early (and cheaply) tell if a main thruster contributes to your direction vector and how much. This is the length of the projected vector basically.

Also notice that thrusters never generate thrust backwards, so what you can do

// these vectors are not unit; the result is a % of velocity's magnitude
// (i.e. 1 would mean all 100% of thruster is contributing, 0.5 would mean only 50% of thruster is contributing)
float getThrusterContribution(Vector3 thruster, Vector3 velocity) {
  return MathF.Max(0f, Vector3.Dot(thruster, velocity));
}

The problem here is when velocity differs from your thruster’s capacity (for example if your velocity is twice that of your thruster’s output, you’d get a maximum of 50% even though they are perfectly aligned), so if you want only a geometric solution you need to pass unit vectors (thruster.normalized and velocity.normalized).

This is enough to attenuate the VFX of the thruster. But if you need to find the wasted thrust component you can do

float decomposeThruster(Vector3 thruster, Vector3 velocity, out Vector3 waste) {
  var proj = Vector3.zero;
  var tm = thruster.magnitude;
  var dp = Vector3.Dot(thruster / tm, velocity.normalized); // thruster / tm equals thruster.normalized
  if(dp > 0f) proj = Vector3.Project(thruster, velocity);
  waste = thruster - proj;
  return proj.magnitude / tm; // returns contribution % as before
}

I’m doing this in the browser, so this could be completely wrong, but I can explain the reasoning behind it.

The dot will take the geometric relationship between two directions. If these directions are perpendicular or opposite, there is no valid contribution, otherwise the vector projection is computed (considering true vector magnitudes).
This projection is collinear with the original velocity, so we have locked the two vectors on their common plane.

Finally, because we know the common plane, and sum of vector components is a parallelogram on this plane, we can find the other component by subtracting projection from the original thrust vector.


I strongly advise to you to get better with quaternions (and linear algebra in general). If you don’t get them, maybe work on something non-physics based until you’re ready, because it is really important that you stay on top of this and not feel lost at every step. Here’s a crash course of mine.

Check out my Spaceship3D demo in my Proximity Buttons project. It is a first person thrust-and-turn spaceship controller that operates purely with a single massive thruster and individual axis torque controls.

The torques are produced as the output of PD filters for each axis (roll, pitch and yaw) and with a little bit of damping it is pretty controllably fun.

The code is (largely) in one class broken into multiple partial class files. Just run from the Spaceship3D scene.

9870384--1422675--Screenshot-MyGame-133618969160312690.png

proximity_buttons is presently hosted at these locations:

https://bitbucket.org/kurtdekker/proximity_buttons

Hi, thanks for the explanation! My navigation system is quite complex, but not as complex as what you describe. It’s not a simulator. Still, the way I made my implementation is an interesting little story. I’ll tell it to you if you want.

That’s exactly what I need. It’s in fact the way you said: it’s easy but at the same time complex. The navigation system will be completely separated from the particle system. As you assumed, I just want to convey information, make it look decent. That’s it:)

Damn, I always forget about the dot product! I used it in this project in the past, but when I stop using it, I forget about its existence. Super useful tool. Thanks for the reminder!

This is a very useful implementation tip I was handling this case in a different way in my code, but yours is much more elegant.

You lost me here. Are you talking about paying attention to the fact that max module of the velocity doesn’t necesserely match with the maxStartLifetime of the rear thruster? If yes, then yes, I haven’t forgot about it. I definetly need to make a commensurate.

Why you’re doing this isn’t obvious for me. I need to think about this a bit.

As always, I forget how thungs when I don’t use them. I remember that the properties of quaternions are very interesting and useful. Defintely will do!

Oh man, thanks for taking the time. That’s just great. I wish I had a person like you working alongside me :wink:

Could you elabote how this is useful in my case? Since I want to change the startLifeTime of a ParticleSystem

@orionsyndrome This code does the job:

public class EnergyBlast : MonoBehaviour
{
    private ParticleSystem _particleSystem;
    private Rigidbody _unitRigidBody;
    private float _unitMaxSpeed;
    
    [SerializeField] private float maxStartLifetime = 0.2f;
    private void Start()
    {
        _particleSystem = ComponentUtils.RequireComponent<ParticleSystem>(this);
        _unitRigidBody = ComponentUtils.RequireComponent<Rigidbody>(transform.root);
        _unitMaxSpeed = ComponentUtils.RequireComponent<ShipAttributes>(transform.root).maxSpeed;
    }

    private void Update()
    {

        Vector3 facingDirection = transform.root.transform.eulerAngles;
        Vector3 forward = Quaternion.Euler(facingDirection) * Vector3.forward;


        var main = _particleSystem.main;
        main.startLifetime = getThrusterContribution(forward, _unitRigidBody.velocity) * maxStartLifetime;
        
    }
    
    float getThrusterContribution(Vector3 facingDirection, Vector3 velocity) {
        Debug.Log($"facingDirection.normalized = {facingDirection.normalized}");
        Debug.Log($"Prodotto scalare = {Vector3.Dot(facingDirection.normalized, velocity.normalized)}");
        Debug.Log(
            $"MathF.Max(...) = {MathF.Max(0f, Vector3.Dot(facingDirection.normalized, velocity.normalized))}");
        return MathF.Max(0f, Vector3.Dot(facingDirection.normalized, velocity.normalized));
    }
    
}

There’s also some flickering in the thrust, but by smoothing the transition with Lerp, things shoud be fine:smile:

1 Like

Great, I’m glad that I covered the simple VFX case as well.

Don’t do this btw

Vector3 facingDirection = transform.root.transform.eulerAngles;
Vector3 forward = Quaternion.Euler(facingDirection) * Vector3.forward;

That’s a pretty roundabout (and fairly expensive) way to just do

var forward = transform.root.transform.forward;

It’s already there.
Another way to obtain it would be

var xf = transform.root.transform;
var localForward = xf.localRotation * Vector3.forward;
var globalForward = xf.rotation * Vector3.forward;

This is the likely cause of that flicker btw due to numerical imprecisions that accumulate throughout the conversion. Euler angles are sloppy and unstable. In fact, don’t ever use Euler angles if you can help it, you only need them if you need actual angles in degrees (for display or for a more readable explicit rotation in code).

I don’t even use the Euler method, I nearly always use AngleAxis, but via shortcuts like so

var arbitraryRot = rotationY(TAU / 5f); // 72 degrees (TAU = 2PI)
 or
var arbitraryRot = rotationY(72f * TAU / 180f); // 72f * Mathf.Deg2Rad

The difference in the amount of instructions emitted through this vs the Euler method is immense.
Also you get to combine rotations in order of your preference

var arbitraryRot = rotationY(TAU / 5f) * rotationX(TAU / 4f); // first X then Y
/// <summary> Returns a quaternion rotation around X axis. </summary>
[MethodImpl(INLINE)] static public Quaternion rotationX(float rad)
  => angleAxis(rad, v3_x);

/// <summary> Returns a quaternion rotation around Y axis. </summary>
[MethodImpl(INLINE)] static public Quaternion rotationY(float rad)
  => angleAxis(rad, v3_y);

/// <summary> Returns a quaternion rotation around Z axis. </summary>
[MethodImpl(INLINE)] static public Quaternion rotationZ(float rad)
  => angleAxis(rad, v3_z);

/// <summary> Constructs a quaternion based on angle (in radians) and axis of rotation (unit vector). </summary>
static public Quaternion angleAxis(float rad, Vector3 axis) {
  const float TOO_SMALL = 3E-8f;
  if(sqrMag(axis) >= TOO_SMALL) {
    rad *= half;
    return v4_xyz(mul(sin(rad), axis), w: cos(rad)).AsQuat();
  }
  return quat_id;
}

Here’s a rewrite if you can’t cope with my style

// axis is supposed to be unit vector
// positive angles rotate counter-clockwise (left handed system)
static public Quaternion angleAxis(float rad, Vector3 axis) {
  const float TOO_SMALL = 3E-8f;
  if(axis.sqrMagnitude >= TOO_SMALL) {
    rad *= half;
    var s = MathF.Sin(rad);
    return new Quaternion(
      s * axis.x,
      s * axis.y,
      s * axis.z,
      MathF.Cos(rad);
    );
  }
  return Quaternion.identity;
}

This is what quaternion is, in its essence, it bakes in trigonometric relations with some extradimensional fluff (stored in w). It’s a compact piece of a linear matrix (that is yet to unfurl into an actual 3x3 rotation matrix) that you can apply to vectors to reliably preserve the idea of an arbitrary rotation in 3D. You can think of x, y, z, as if they’re Euler sliders, but with sines of half angles in radians. However they are not independent, when one of these change, a quantity in w changes as well, so there is this conservation because rotation quaternion are always unit (x^2 + y^2 + z^2 + w^2 = 1). If you pay close attention, w is practically a register for “no-rotation”, hence identity = 0,0,0,1 (i.e. “do nothing”)

Edit:
Also btw, try to cache frequently accessed transforms in a variable, don’t rely on constantly traversing the tree with transform.root.transform etc. Ideally you hook your object directly in the inspector, then locally cache its transform in Awake. Assuming that the root object is always the valid target is a bad practice prone to accidents.

What I posted is obviously lacking a particle system, but since I’ve made about 500 spaceship control games in my life, every one that has a PDFilter has always felt better than the ones that didn’t, and I thought you might care to know about the process.

If not, well, move on!

Btw, you can slightly improve this expression (if you don’t already have normalized versions lying around), if you know that

dot(a, b) = length(a) * length(b) * dot(normalize(a), normalize(b));

This is why you need to normalize the vectors, to get rid of this anomaly, however normalization is a bit expensive because it’s computed as sqrt(a.x * a.x + a.y * a.y + a.z * a.z) and you have to do it twice, just to take the dot (which is cheap).

So if you know the above, you can do

dot(n(a), n(b)) = dot(a, b) / length(a) * length(b)
dot(n(a), n(b)) = dot(a, b) / sqrt(sqrlength(a) * sqrlength(b))

right? This is just one square root, instead of two.

Thus we can implement a more optimal variant of this

float dotNormalized(Vector3 a, Vector3 b)  => Vector3.Dot(a, b) / MathF.Sqrt(a.sqrMagnitude * b.sqrMagnitude);

// finally
return MathF.Max(0f, dotNormalized(facingDirection, velocity));

All of this is funny to me, because sqrMagnitude is already a kind of dot product as well. It’s a self dot-product. Look

float sqrt(float n) => MathF.Sqrt(n);
float sqrMag(Vector3 v) => dot(v);
float mag(Vector3 v) => sqrt(sqrMag(v));
Vector3 normalize(Vector3 v) => v / mag(v);
float dot(Vector3 a, Vector3 b) => a.x * b.x + a.y * b.y + a.z * b.z;
float dot(Vector3 v) => v.x * v.x + v.y * v.y + v.z * v.z;
float dotNormalized(Vector3 a, Vector3 b) => dot(a, b) / sqrt(dot(a) * dot(b));

So 10 multiplications, 1 division, 1 square root

whereas the previous solution

dot(a / sqrt(dot(a)), b / sqrt(dot(b)))

had 9 multiplications in total, but 6 divisions and 2 square roots

This matters because this kind of code will always lie in your hot paths. It’s not something you do once and forget about it, this is something you do in every frame, and likely many times over. Frankly divisions aren’t that much more expensive than multiplication (which is brutally optimal), but every single shaved off operation counts.