How to generate random velocity vectors that can only move an object forward within a valid arc?

Question

I have an object with known coordinates in in 3D but on the ground (z=0). The object has a direction vector. My goal is to move this object on the ground (so z stays 0) using its direction vector and via randomly-generated velocity vectors with one condition: I want to ensure that the generate velocity vector can only move the object within a "valid arc", defined in degrees with respect to the direction vector of the object. More specifically, the way I determine valid range is by ensuring that the new position is within 45 degrees of the old object's direction vector (-45 degrees to the left and +45 degrees to the right). Can someone write a pseudocode on how I can achieve this?

Here's my attempt to do this but this doesn't seem to be the correct way to help me achieve what I want:

object_dir = object_position # the direction could be the same as the object's coortinates
while True:
    vel_vec=[uniform(-max_vel, max_vel), uniform(-max_vel, max_vel)] # generate a random velocity vector
    new_pos = object_dir + vel_vec # compute a new position (and/or object direction vector) for the object
    if (compute_angle(new_pos, object_dir) < 45 or compute_angle(new_pos, object_dir) > 315):
        break 

Answer 1

Based on the comments, I think you just want the direction of motion in each time step to be less than $45$ degrees from the direction of motion in the previous time step. That means you can completely reverse direction if you take enough steps to do it. You can also eventually travel into all four quadrants of the plane regardless of how you started out from the origin.

Moreover, the velocity of the object is a vector in the direction of travel of the object. The only distinction I might make between a "direction vector" and velocity is that a "direction vector" might be normalized to a certain length. Often the direction vector will be defined as the unit-length vector in the same direction as the velocity.

If you can find the angle between arbitrary vectors, not just unit vectors, I don't even see a reason to include any "direction vector" in your code. Just consider velocity.

So I would ignore positions entirely when computing the angle. But you need to remember what the direction vector was in the previous time step.

---

Here's one way you could do it. In each time step, you compute a new vel_vec. You then compute the angle between vel_vec and previous_vel_vec (which was the velocity vector in the previous time step). If that angle is in the permitted range, accept vel_vec, otherwise try again.

When you finally find an acceptable vel_vec, then you use it to compute a new position. And then make the assignment previous_vel_vec = vel_vec so that you will have the correct previous_vel_vec for the next time step.

---

But I would do it quite differently. For one thing, I don't like the fact that the velocity is more likely to be at about a $45$ degree angle from an axis than parallel to an axis regardless of what the valid arc is at that time. This happens because you distributed the velocity uniformly over a square, and there is more area in the parts of a square within $\pm22.5$ degrees of a diagonal than within $\pm22.5$ degrees of a centerline parallel to a side. I would prefer to distribute the direction uniformly within the limits of the valid arc.

I would keep track of the "compass direction" of the velocity (that is, what angle it makes with the $x$ axis in the $x,y$ plane). Then, at each time step, to find the new velocity I just add a random number uniformly distributed between $-\frac\pi4$ and $\frac\pi4$ (measuring the angle in radians -- if you insist on using degrees, which you will continually have to convert to radians anyway for every trig function, then you want a uniform distribution between $-45$ and $45$). Take another random number $\rho$ uniformly distributed between $0$ and $1,$ and multiply max_vel by $\sqrt{\rho},$ because that's how we get a uniform distribution within a circular sector. (As this answer observes, we get a uniform distribution within an entire circle if we let the angle range over $360$ degrees; we are just restricting the angle in order to sample a sector).

Answer 2

Pick a random angle $\theta$ between $-45^\circ$ and $+45^\circ$ (or equivalently $-\pi/4$ and $+\pi/4$), and rotate the current direction vector by that angle around the $z$ axis.

If the surface your object is moving on is not perpendicular to the $z$ axis, rotate around the unit surface normal of that surface.

One very efficient way to do this is to use the Rodrigues' rotation formula, $$\vec{v}_\text{new} = \vec{v}_\text{old} \cos \theta + (\vec{n} \times \vec{v}) \sin\theta + \vec{n} \left(\vec{n} \cdot \vec{v}_\text{old} \right)(1 - \cos\theta)$$ where $\vec{n}$ is the surface unit normal vector ($\lVert\vec{n}\rVert = 1$), $\vec{v}_\text{old}$ is the old direction/velocity vector, and $\vec{v}_\text{new}$ is the direction/velocity vector rotated by $\theta$ around unit axis vector $\vec{n}$.

---

Here is an example implementation in Python:

import math
import random

def turn(direction, normal, max_degrees, distribution=None):
    # Cartesian coordinate components of the two vectors
    dx, dy, dz = direction
    nx, ny, nz = normal

    # Make sure the normal vector is an unit vector.
    n = math.sqrt(nx*nx + ny*ny + nz*nz)
    if n == 0.0:
        raise ValueError("Cannot turn around a zero vector!")
    if n != 1.0:
        nx = nx / n
        ny = ny / n
        nz = nz / n

    # Pick an uniform random number between -1 and +1,
    u = random.Random().uniform(-1, +1)
    # apply the distribution function (0..1 -> 0..1) if any is specified,
    if distribution is not None:
        u = math.copysign(distribution(abs(u)), u)
    # and scale to radians, since that's what math.sin and math.cos take,
    rads = u * abs(max_degrees) * math.pi / 180.0
    s = math.sin(rads)
    c = math.cos(rads)

    # Cross product between normal (nx,ny,nz) and direction (dx,dy,dz)
    cx = ny*dz - nz*dy
    cy = nz*dx - nx*dz
    cz = nx*dy - ny*dx

    # Dot product between normal and direction, scaled by (1-cos)
    d = (nx*dx + ny*dy + nz*dz) * (1 - c)

    # Apply Rodrigues' rotation formula to dx,dy,dz
    rx = dx*c + cx*s + d*nx
    ry = dy*c + cy*s + d*ny
    rz = dz*c + cz*s + d*nz

    return rx,ry,rz


if __name__ == '__main__':
    uniform = random.Random().uniform

    # Pick a random direction on the XY plane,
    xyangle = uniform(0, 2*math.pi)
    direction = (math.cos(xyangle), math.sin(xyangle), 0)

    print("Original direction: (%.6f, %.6f, %.6f)" % direction)

    # Deflected direction can vary by up to 45 degrees.
    deflected = turn(direction, [0,0,1], 45.0)
    print("Turned direction: (%.6f, %.6f, %.6f)" % deflected)

    # Calculate and report the actual deflection angle.
    theta = math.acos( (direction[0]*deflected[0] + direction[1]*deflected[1] + direction[2]*deflected[2])
                       / math.sqrt( (direction[0]*direction[0] + direction[1]*direction[1] + direction[2]*direction[2])
                                  * (deflected[0]*deflected[0] + deflected[1]*deflected[1] + deflected[2]*deflected[2]) ))
    print("Turning angle %.3f degrees" % (theta * 180.0 / math.pi))

You can run the above example to see a random direction (in the XY plane) rotated up to $45^\circ$ around the $z$ axis.

The optional distribution argument to the turn function is a single-variate function applied to an uniform random number between $0$ and $1$ to get the desired distribution (also between $0$ and $1$). An extremely common one in Python would be lambda x: x**N where N is some constant number greater than 1, say around 2.0, so that small angles are more common than large angles. (The larger N is, the more common small angles are compared to large angles. Conversely, N between $0$ and $1$ would make larger turns more common than small turns the closer N is to zero.)

Alternatively, in Python one could use random.Random().betavariate() or random.Random().gauss() instead of random.Random().uniform() to generate a deflection scale value between $-1$ and $+1$ in the desired (symmetric?) distribution.