Skip to main content
gamedev.net
🔒 Locked

math for a sidewinder-like missile

Started by fluke 9 replies 3.7k views
Original Post
fluke
fluke
i'll start with the source. be careful not to fall out of your chair when you start laughing..

class Missile { 
public:
    float x,x2,y,y2,z,z2;
    float curAngle, destAngle;

    Missile(){
        x=0;
        y=0;
        z=0;
        curAngle = 0.0f;
    }

    void Init(float x1, float y1, float z1, float angle){
        x=x1;
        y=y1;
        z=z1;
        curAngle = angle;
    }

    void Update(float x2, float y2, float z2, float speed){        
        float destAngle = atan2((z2-z),(x2-x));

        if(curAngle>destAngle+PI) curAngle -= 1.0f * PI;
        if(curAngle<destAngle-PI) curAngle += 1.0f * PI;

        curAngle += (destAngle - curAngle) * .1f;
        
        x += cosf(curAngle) * speed;
        y = 0; //???
        z += sinf(curAngle) * speed;                
    }
};

ooookay. quite plainly, i suck at math and found this 2d implementation that has some problems. first of all, in its simplicity, it doesn't even do what i want, which is to behave like a sidewinder missile eradically painting the sky with plumes of smoke using its proportional navigation before finally plowing into its target. instead, it "attempts" to curve around and eventually hit the target, as opposed to going in a straight, boring line. however, it even fails at this in some ways. when it reaches the target, the x and z values go up and down causing a "ghosting" type effect where there are seemingly two instances of the missile co-existing a short distance from one another. finally, the last problems is that i don't know how to solve for y. anyone care to help or better yet provide some sample code on how to create a more accurate sidewinder effect? :)
Emergent
Emergent
You don't need angles here, and they are just making your life complicated. Just use the missile's position and velocity (possibly normalized) as its state, rather than the position and angle.

Now,
1 - Your current heading is a unit vector. I'll call it v. It lives on a unit sphere.
2 - If your position is the vector p and the target's is z, then your goal heading is g=normalize(z-p). It also lives on a sphere.

So your problem is actually: How do I implement a proportional controller on a sphere?

(I'm assuming that your "real" problem is 3d here.)

The axis of rotation you always want to instantaneously rotate about is,
a = normalize(v x g)

and the instantaneous rate of change in v that such a rotation would cause is the tangent vector at v to the geodesic connecting it to g on the sphere; this is,

t = a x g

so you simply want to always adjust v in the direction of t. By how much?

The angle is,

theta = asin(norm(v x g))

so your controller can be,

dv/dt = k*theta*(a x g)

where k is your controller gain (some positive constant).

Note that this is in the perfect world of differential equations where your stepsizes are infinitely small. But if you just integrate this using Euler integration, norm(v) will slowly get bigger and bigger. You can just snap it back to the sphere after each step, like

v = v + dt*k*theta*(a x g)
v = normalize(v)

which is the simplest implementation I can think of right now. :-)
Nanoha
Nanoha
I agree, vectors will definatly make your life so much easier.

I couldn't quite understand if the missile is supposed to home in or not, I guess it is otherwise you wouldn't have any problem. Definatly drop the angles, make a simple vector class:

class Vector
{
public:
float x;
float y;
float z;
};

its going to make your life ALOT easier in the long run.

Now give your missile position, velocity and acceleration vectors.

When you fire your missile, set its position to where it is, set the velocity to the direction it is facing (same as acceleration). Magnitudes of acceleration and initial velocity will depend, trial and error.

You want your missile to go in a particular direction (another vector, desired velocity). Just set your acceleration in the direction you want to go and each update do:

velocity += acceleration*timePassed
position += velocity*timePassed;

where timePassed is the time since the last update (otherwise things will speed up and slow down with your frame rate).

The problem with that is that eventually your missile will travel super fast so you need some kind of drag (like air resistance) and you can do a very simple method (though not physically correct) by using current velocity * a drag factor like:

Vector drag = -velocity*0.1;
Vector finalAcceleration = (drag+acceleration);
velocity += finalAcceleration*timePassed;
position += velocity*timePassed

Its not a perfect solution but it is pretty simple and looks good.
Interested in Fractals? Check out my App, Fractal Scout, free on the Google Play store.
fluke
fluke
Emergent, you're using a few expressions i don't understand. what is "x" in "v x g" and "a x g"?

Nanoha, how can i get the missile to curve using your approach? and to answer your question, yes, i guess it is effectively a homing missile. sidewinders (which i assume get their name from their behavior which is similar to the snake of the same name) overcompensate their turns to "predict" the position of their target. their curves are wide and windy when initially fired and the path is less windy as it closes in on the enemy.

i'll try to get both of these implementations working and any others that anyone provides.
Trihex
Trihex
The x most likely stands for the cross product of the two vectors.

You can find more information on the cross product at :

Cross Product
Eric ArmstrongTrihex Softwarewww.trihex.com
Emergent
Emergent
I used 'x' to denote cross product.

It's funny that you want overshoot. I understand why; it'd look cool. It's just that normally overshoot is a bad thing for a controller!

The controller I proposed will not overshoot. It will overshoot if you use a (badly tuned. :-)) PID controller instead of a proportional controller.

How do you build a PID controller on a sphere? Hmmm. That's kind of interesting! Let's see. The most natural representation for your error is,

e = theta*a

so I guess the way to do a PID controller on this manifold is,

dv/dt = (k1*e + k2*integral(e) + k3*(d/dt)e) x g .

where really "integral" should probably not be an actual integral, but a "leaky integral," like

dI/dt = e - k4 I.

If k4 were zero, 'I' would be the integral of e. Choosing a small positive k4 will make your integral "leaky" or a little "forgetful," and this is a good thing.

Then, the "right" way to discretize this is really with a thing called the "Bilinear Transform" (or maybe I should have just modeled the system as a discrete-time one to begin with) but that's probably way more than you're interested in. So my guess is that using a sum to approximate the integral and finite difference to approximate the derivative should work.
Emergent
Emergent
Actually you don't need the 'D' part of a PID controller here. Just a PI (with leaky I -- which is really called a "lag compensator") should do the trick.

[EDIT: Giving this more thought, I'm quite convinced that this will indeed work -- nicely.]

[Edited by - Emergent on September 3, 2009 9:00:28 AM]
Nanoha
Nanoha
Quote:
Original post by fluke
Nanoha, how can i get the missile to curve using your approach?


if you missile is flying in direction (0, 0, 1) but your target is in direction (1, 0, 0) then all you do is accelerate in that direction.

(0, 0, 1) + (1, 0, 0) = (1, 0, 1) which is sorta half way between the current and the desired.
next update more so:
(1, 0, 1) + (1, 0, 0) = (2, 0, 1)

You'll notice the velocity is getting high and still going down the z axis, thats where the drag I mentioned comes in, it'll eventually reduce that to 0 so your just flying down the x axis (which is where the target is).

The missile will automaticly curve towards its target.

The only "problem" is drawing your missile but I couldn't see from your post how you were representing orientation. You can get the orientation of your missile from its velocity (assuming it will always be facing the way it is traveling).
Interested in Fractals? Check out my App, Fractal Scout, free on the Google Play store.
Emergent
Emergent
Quote:
Original post by Nanoha
The only "problem" is drawing your missile but I couldn't see from your post how you were representing orientation. You can get the orientation of your missile from its velocity (assuming it will always be facing the way it is traveling).


Yeah, I didn't address this either...

Velocity will not uniquely give you an orientation; consider spinning your missile about an axis parallel to its velocity. If your missile is axially symmetric though, then velocity would be enough information.

One thing you could do would be to use the Frenet-Serret frame as your orientation. Again, this is not always unique. You can deal with this by, whenever a Frenet-Serret vector is not uniquely specified, using whatever the last value for that vector was. Following this convention, you now have unique orientations.

The real answer to this problem though would be to just keep track of the missile's orientation directly. Then, your controller would always just rotate your orientation about the axis I talked about.

So whereas for the lag compensator, working with just velocities you'd have,

dv/dt = (k1*e + k2*I) x g
dI/dt = e - k4 I ,

if you work with orientations you instead have,

dR/dt = hat(k1*e + k2*I) R
dI/dt = e - k4 I

where R is a rotation matrix giving the orientation of your missile and 'hat' refers to the hat operator; it takes a vector and returns a matrix.

Note that you're not going to want to implement the dR/dt ODE with Euler integration; you won't keep your rotation matrix a rotation matrix, among other things. Rather, you'll want to either

1 - use the exponential map, like

R = exp((dR/dt) deltat) R

implementing 'exp' using Rodrigues' formula (see here) and re-orthonormalize R occasionally, or,

2 - represent your rotation matrix R by a quaternion q. Use the chain rule to figure out what dq/dt is and just integrate that using Euler integration, always snapping q to the unit sphere after each step.

[Edited by - Emergent on September 3, 2009 9:59:26 PM]
fluke
fluke
ok, i've been working on this steadily for the last few days, but with no luck. i found an example of kind of what i'm after in a game programming book i have and it uses a PID controller like suggested... however, it's in 2D, so it uses angles for it's heading, which in the code below you see i've attempted to correct. then again below, it attempts a direction rotation based on an angle, which obviously is not going to work for 3d. i'd need to specify an axis for rotation, but i don't know how to determine the rotations for each axis based on the previous calculations.

//PART OF PID:void Missile::Clear(){	m_CurrentIndex		= -1;	m_PreviousIndex		= -1;	m_NumErrorsRecorded	= 0;	m_CurrentIntegral	= 0.0f;	for (int i = 0; i < NUM_ERROR_SLOTS; i++)	{		m_Error		= 0.0f;		m_Timestep	= 0.0f;	}}void Missile::Record(float error, float timestep){	m_PreviousIndex				= m_CurrentIndex;	m_CurrentIndex				= (m_CurrentIndex + 1) % NUM_ERROR_SLOTS;	if (m_NumErrorsRecorded == NUM_ERROR_SLOTS){		m_CurrentIntegral		-= m_Error[m_CurrentIndex] * m_Timestep[m_CurrentIndex];	}	m_Error[m_CurrentIndex]		= error;	m_Timestep[m_CurrentIndex]	= timestep;	m_CurrentIntegral			+= error * timestep;	m_NumErrorsRecorded			= min(m_NumErrorsRecorded + 1, NUM_ERROR_SLOTS);}float Missile::GetError(){	if (m_NumErrorsRecorded >= 1){		return m_Error[m_CurrentIndex];	}else{		return 0.0f;	}}float Missile::GetErrorIntegral(){	return m_CurrentIntegral;}float Missile::GetErrorDerivative(){	if (m_NumErrorsRecorded >= 2){		float difference	= m_Error[m_CurrentIndex] - m_Error[m_PreviousIndex];		float time_interval	= m_Timestep[m_CurrentIndex];		if (time_interval > 0.001f){			return (difference / time_interval);		}else{			return 999999.0f;		}	}else{		return 0.0f;	}}float Missile::GetOutput(){	m_P_Coefficient = 2.0f;				// Our current P coefficient	m_I_Coefficient = .5f;				// Our current I coefficient	m_D_Coefficient = 2.9f;				// Our current D coefficient	return m_P_Coefficient * GetError() + m_I_Coefficient * GetErrorIntegral() + m_D_Coefficient * GetErrorDerivative();}float Missile::CalculateIntegralDirectly(){	float integral = 0.0f;	for (int i = 0; i < m_NumErrorsRecorded; i++){		integral += (m_Error * m_Timestep);	}	return integral;}//MISSILE MOVEMENT:void Missile::Move(float timestep){        	//Vector3 vector_to_target = vTargetPos - vMissilePos;//Vector3 vector_to_target = *(m_pTarget->GetPosition()) - vMissilePos;	//float heading_error = vDirection.GetAngle() - vector_to_target.GetAngle();            	Vector3 vector_to_target = vTargetPos - vMissilePos;        		    Vector3 v1 = vDirection;    Vector3 v2 = vector_to_target;    v1.normalize();    v2.normalize();    float heading_error = RAD2DEG(asin(Vector3::dot(&v1,&v2)));		if (heading_error > 180.0f){		heading_error -= 360.0f;	}else if (heading_error < -180.0f){		heading_error += 360.0f;	}	Record(heading_error, timestep);    float	desired_acceleration			= GetMaxAcceleration();	float	desired_angular_acceleration	= GetOutput();	desired_acceleration					= Clamp(desired_acceleration,			0.0f,		1000.0f);	desired_angular_acceleration			= Clamp(desired_angular_acceleration,	-180.0f,	180.0f);	SetAcceleration(desired_acceleration);	SetAngularAcceleration(desired_angular_acceleration);    //move	m_Speed += m_Acceleration * timestep;	Vector3 velocity = vDirection;	velocity.normalize_especial(m_Speed);	m_AngularVelocity += m_AngularAcceleration * timestep;	float drag		= -m_Speed * (float)fabs(m_Speed);	drag			*= (MissileDragFactor * timestep);	m_Speed			+= drag;	float rotational_drag	= -m_AngularVelocity * (float)fabs(m_AngularVelocity);	rotational_drag			*= (MissileRotationalDragFactor * timestep);	m_AngularVelocity		+= rotational_drag;		float delta_angle = m_AngularVelocity * timestep;		vDirection.Rotate(delta_angle); //not correct for 3d	vDirection.normalize();		vMissilePos += velocity * timestep;    Vector3 p = vMissilePos;    p.normalize();    x = p.x;    y = p.y;    z = p.z;}
fluke
fluke
ok, i re-wrote the whole thing and got it working-ish. i ditched the pid and just went with a simpler implementation like Nanoha suggested and using a direction and velocity vector to calculate the new pos. anyway, thanks for the help everyone!

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.

Sign in to reply to this topic.