[SOLVED] Help implementing complex numbers

I've finished coding the equations for quartic and cubic polynomials, but I'm running into an issue during execution of the code since the equation calculates a complex value. I have a class that I found online for using complex numbers, but I'm not exactly sure how to implement it in the script to actually return the roots in terms of (i) for those that are imaginary. The code for the cubic equation is where I've narrowed down my results to (by use of adding to CustomData and throwing exceptions to let me know where it got to). This is that code:

Code:

public double[] GetAllCubicRoots(double A, double B, double C, out double[] root)
{
    root = new double[3];

    // Depress the cubic to the form t^3 + Pt + Q = 0
    // This is done by substituting (t - A/3) for x
    double P = (((3 * B) - (A * A)) / 3);
    double Q = (((2 * A * A * A) - (9 * A * B) + (27 * C)) / 27);

    // If P and Q both equal zero, then t = 0, and x = -A/3
    // I don't think this applies to this situation, but I'll leave it just in case...
    if (P == 0 && Q == 0)
    { root[0] = -A / 3; root[1] = 0; root[2] = double.NaN;
        Me.CustomData = root.ToString();
        throw new Exception("Check Root. All Cubic Roots P And Q = 0.");
        return root;
    }

    // If P OR Q equal zero, the cubic is immediately solvable
    if (P == 0)
    {
        // The equation becomes t^3 + Q = 0
        // t = cuberoot(-Q) - (A/3) AND t = cuberoot(-Q) * the Cube Roots of Unity - (A/3)
        double cbrtNQ = CubeRoot(-Q);
        double sqrt3iOn2 = Math.Sqrt((3 * i.Imaginary) / 2);

        double root1 = cbrtNQ - (A / 3);
        double root2 = cbrtNQ * (-0.5 + sqrt3iOn2) - (A / 3);
        double root3 = cbrtNQ * (-0.5 - sqrt3iOn2) - (A / 3);

        if (root1 <= 0 && root2 <= 0 && root3 <= 0)
        {
            root[0] = root[1] = root[2] = double.NaN;
            Me.CustomData = root.ToString();
            throw new Exception("Check Root. All Cubic Roots P = 0. All roots <= 0.");

            return root;
        }
        else
        {
            root[0] = root1;
            root[1] = root2;
            root[2] = root3;
            Me.CustomData = root.ToString();
            throw new Exception("Check Root. All Cubic Roots P = 0. Some roots > 0.");

            return root;
        }
    }
    else if (Q == 0)
    {
        // The equation becomes t^3 + Pt = 0
        // Factor out a t and solve for t(t^2 + P) = 0. Roots are (-A/3) and the roots of the equation t = +- sqrt(-P) - (A/3)
        double sqrtNP = Math.Sqrt(-P);

        double root1 = -A / 3;
        double root2 = sqrtNP - (A / 3);
        double root3 = -sqrtNP - (A / 3);

        if (root1 <= 0 && root2 <= 0 && root3 <= 0)
        {
            root[0] = root[1] = root[2] = double.NaN;
            Me.CustomData = root.ToString();
            throw new Exception("Check Root. All Cubic Roots Q = 0. All roots <= 0.");

            return root;
        }
        else
        {
            root[0] = root1;
            root[1] = root2;
            root[2] = root3;
            Me.CustomData = root.ToString();
            throw new Exception("Check Root. All Cubic Roots Q = 0. Some roots > 0.");

            return root;
        }
    }

    // The discriminant (D) disctates the number and type of roots for the equation
    double discriminant = Math.Pow(Q / 2, 2) + Math.Pow(P / 3, 3);
    double sqrtD = Math.Sqrt(discriminant);

    // If D > 0, one real and two imaginary roots. Skip the imaginary roots.
    if (discriminant > 0)
    {
        double root1 = CubeRoot((-Q / 2) + sqrtD) - CubeRoot((Q / 2) + sqrtD) - (A / 3);
        Me.CustomData = root1.ToString() + "\n";
        if (root1 > 0) { root[0] = root1; root[1] = root[2] = double.NaN; }
        else { root[0] = root[1] = root[2] = double.NaN; }
        Me.CustomData += root[0].ToString() + "\n" + root[1].ToString() + "\n" + root[2].ToString();
        throw new Exception("Check Root. All Cubic Roots D > 0.");

        return root;
    }
    // If D = 0, All real, but duplicate roots
    else if (discriminant == 0)
    {
        double root1 = 2 * CubeRoot((-Q / 2)) - (A / 3);
        double root2 = CubeRoot((Q / 2)) - (A / 3);

        if (root1 <= 0 && root2 <= 0)
        {
            root[0] = root[1] = root[2] = double.NaN;
            Me.CustomData = root.ToString();
            throw new Exception("Check Root. All Cubic Roots D = 0. All roots <= 0.");

            return root;
        }
        else
        {
            root[0] = root1;
            root[1] = root2;
            Me.CustomData = root.ToString();
            throw new Exception("Check Root. All Cubic Roots D = 0. Some roots > 0.");

            return root;
        }
    }
    // If D < 0, All real and distinct roots, but u and v will be imaginary
    else
    {
        // Uses the formula: u/v = cbrt(R) * ( +/- cos(theta/3) + (i * sin(theta/3)) to find real roots from imaginary numbers
        double R = Math.Sqrt(Math.Pow((-P / 3), 3));
        double cbrtR = CubeRoot(R);
        double phi = Math.Acos(-Q / (2 * R));
        double twoPI = MathHelper.TwoPi;
        double fourPI = MathHelper.FourPi;

    Me.CustomData = "A = " + A + "\n";
    Me.CustomData += "B = " + B + "\n";
    Me.CustomData += "C = " + C + "\n\n";
    Me.CustomData += "P = " + P + "\n";
    Me.CustomData += "Q = " + Q + "\n";
    Me.CustomData += "R = " + R + "\n";
    Me.CustomData += "cbrtR = " + cbrtR + "\n";
    Me.CustomData += "phi = " + phi + "\n\n";

        double root1 = 2 * cbrtR - (A / 3);
        double root2 = 2 * cbrtR * Math.Cos((phi + twoPI) / 3) - (A / 3);
        double root3 = 2 * cbrtR * Math.Cos((phi + fourPI) / 3) - (A / 3);
    Me.CustomData += root1.ToString() + "\n" + root2.ToString() + "\n" + root3.ToString() + "\n\n";

    if (root1 <= 0 && root2 <= 0 && root3 <= 0)
        {
            root[0] = root[1] = root[2] = double.NaN;
        Me.CustomData += root[0].ToString() + "\n" + root[1].ToString() + "\n" + root[2].ToString() + "\n";
        throw new Exception("Check Root. All Cubic Roots D < 0. All roots <= 0.");

            return root;
        }
        else
        {
            root[0] = root1;
            root[1] = root2;
            root[2] = root3;
            Me.CustomData += root[0].ToString() + "\n" + root[1].ToString() + "\n" + root[2].ToString();
            throw new Exception("Check Root. All Cubic Roots D < 0. Some roots > 0.");

            return root;
        }
    }
}

The place where I have Me.CustomData is where it is ending up, so I SHOULD have 4 real roots for the quartic - three from the calculation, and one default root (x = -B/4). The output from the CustomData is:

Code:

A = -29.8642505691916
B = 31235.4973599879
C = 7.86469020078423

P = 30938.2062059681
Q = 308976.465179612

R = NaN
cbrtR = Nan
phi = Nan

Root1 = NaN
Root2 = NaN
Root3 = Nan

Obviously, I am taking the square root of (-P/3)^3, which will be the square root of a negative value, thus a complex number. I've tried changing the type to:

Code:

//Old type
double R = Math.Sqrt(Math.Pow((-P / 3), 3));

// Changed to
ComplexNumber R = new ComplexNumber(1, Math.Sqrt(Math.Pow((-P / 3), 3)));

But then I get errors all over the place as it won't allow multiplication, addition, and subtraction with doubles and complex numbers. Am I going to have to change EVERY value to a complex number? And am I doing it properly? Here is the class for the complex numbers:

Code:

public class ComplexNumber
{
    // Members of the class
    public double Real { get; set; }
    public double Imaginary { get; set; }

    public double Radius()
    {
        return Math.Sqrt(Real * Real + Imaginary * Imaginary);
    }

    public double Angle()
    {
        return Math.Atan2(Imaginary, Real);
    }

    // Functions
    public string ToString(bool asPolar)
    {
        string str = "";

        if (asPolar)
        {
            // Polar
            if (this.Angle() > 0)
            {
                // Positive angle
                str = this.Radius() + "(cos" + this.Angle() + " + isin" + this.Angle() + ")";
            }
            else
            {
                str = this.Radius() + "(cos" + this.Angle() + " - isin" + -this.Angle() + ")";
            }
        }
        else
        {
            // Rect
            if (this.Imaginary > 0)
            {
                // Positive Imaginary
                str = this.Real + " + i" + this.Imaginary;
            }
            else
            {
                str = this.Real + " - i" + -this.Imaginary;
            }
        }
        return str;
    }

    ComplexNumber Logarithm()
    {
        return new ComplexNumber()
        {
            Real = Math.Log(this.Radius()),
            Imaginary = this.Angle()
        };
    }

    ComplexNumber Conjugate()
    {
        return new ComplexNumber(this.Real, -this.Imaginary);
    }

    ComplexNumber Reciprocal()
    {
        return new ComplexNumber()
        {
            Real = (this.Real) / (this.Real * this.Real + this.Imaginary * this.Imaginary),
            Imaginary = -(this.Imaginary) / (this.Real * this.Real + this.Imaginary * this.Imaginary)
        };
    }

    // Operator overloading
    public static ComplexNumber operator +(ComplexNumber a, ComplexNumber b)
    {
        return new ComplexNumber(a.Real + b.Real, a.Imaginary + b.Imaginary);
    }

    // Negate
    public static ComplexNumber operator -(ComplexNumber a)
    {
        return new ComplexNumber(-a.Real, -a.Imaginary);
    }

    // Subtraction
    public static ComplexNumber operator -(ComplexNumber a, ComplexNumber b)
    {
        return new ComplexNumber(a.Real - b.Real, a.Imaginary - b.Imaginary);
    }

    public static ComplexNumber operator *(ComplexNumber a, ComplexNumber b)
    {
        return new ComplexNumber()
        {
            Real = (a.Real * b.Real) - (a.Imaginary * b.Imaginary),
            Imaginary = (a.Real * b.Imaginary) + (a.Imaginary * b.Real)
        };
    }

    public static ComplexNumber operator /(ComplexNumber a, ComplexNumber b)
    {
        return new ComplexNumber()
        {
            Real = (a.Real * b.Real) + (a.Imaginary * b.Imaginary) /
                ((b.Real * b.Real) + (b.Imaginary * b.Imaginary)),
            Imaginary = (a.Real * b.Real) - (a.Imaginary * b.Imaginary) /
                ((b.Real * b.Real) + (b.Imaginary * b.Imaginary))
        };
    }

    public static bool operator ==(ComplexNumber a, ComplexNumber b)
    {
        // Check whether their values are same
        return (b.Real == a.Real && b.Imaginary == a.Imaginary);
    }

    public static bool operator !=(ComplexNumber a, ComplexNumber b)
    {
        return !(b.Real == a.Real && b.Imaginary == a.Imaginary);
    }

    // Constructors
    public ComplexNumber()
    {
        // Default
        Real = 0;
        Imaginary = 0;
    }

    public ComplexNumber(double firstVal, double secondVal)
    {
        Real = firstVal;
        Imaginary = secondVal;
    }

    public override string ToString()
    {
        return base.ToString();
    }

    public override bool Equals(object obj)
    {
        return base.Equals(obj);
    }

    public override int GetHashCode()
    {
        return base.GetHashCode();
    }
}

Please help! And thank you

 

Why on earth do you need complex numbers?!

 

To calculate the intercept point of an object with another object using non-constant speed, naturally

Ends up being a quartic equation, the roots of which rely on complex numbers to find the real ones in most cases (accurate real ones, that is).

That said, I did find this snippet online last night that implements them, so I can base mine off of it - but their calculations end up being + or - several seconds of the time I actually need to intersect the target (I either fly 50 meters in front of him, or 50 meters behind him). This could be due to the remote and how slow it is on rotation and initial take off (seems like it wobbles quite a bit when settling in on the waypoint before actually moving toward it). Or it could be due to the direction in which the target is moving (ie am I chasing him down, are we flying at orthogonal directions, or are we playing chicken) combined with the distance to target (might be off by some percent, which is scaled with the distance/time needed).

Code:

/// <summary>
/// Solves a quartic equation.
/// </summary>
/// <param name="A">4th power coefficient</param>
/// <param name="B">3rd power coefficient</param>
/// <param name="C">2nd power coefficient</param>
/// <param name="D">1st power coefficient</param>
/// <param name="E">Constant</param>
/// <returns>Any real roots that exist</returns>
public double[] SolveQuartic(double A, double B, double C, double D, double E, out double[] root)
{
    // We need to make sure A isn't too small, for two reasons.  First, we need to avoid division by zero.  Second, small
    // values of A translate into large errors in the root calculation.  You may wish to play around with this value depending
    // upon the range of your inputs.
    if (Math.Abs(A) < 0.1)
    {
        if (A == 0) { A = 0.1f; }
        else { A = Math.Sign(A) * 0.1f; }
    }

    // We'll use Ferrari's method to solve these roots.
    double a = -(3 * B * B) / (8 * A * A) + C / A;
    double b = (B * B * B) / (8 * A * A * A) - (B * C) / (2 * A * A) + D / A;
    double c = -(3 * B * B * B * B) / (256 * A * A * A * A) + (C * B * B) / (16 * A * A * A) - (B * D) / (4 * A * A) + (E / A);

    // It's worth nothing that there are some special cases that arise during this calculation.  For example, 
    // if b = 0, then the quartic equation is actually a bi-quadratic equation, and has very simple roots that
    // are easy to calculate.  However, that situation almost never comes up for us, so creating a case for it
    // doesn't save any computation time.  We're just going to ignore all special cases and solve the general 
    // problem.
    double P = -(a * a) / 12 - c;
    double Q = -(a * a * a) / 108 + (a * c) / 3 - (b * b) / 8;

    // Now we have to step into Complex arithmetic, because there can be negative bases raised
    // to powers less than 1 in this equation, resulting in imaginary roots.
    Complex R = (-Q / 2) + Complex.Sqrt((Q * Q) / 4 + (P * P * P) / 27);
    Complex U = Complex.Pow(R, 1.0f / 3.0f);

    // This conditional statement is necessary to avoid division-by-zero which will occur if |U| = 0.
    // You can't divide by zero, even in the complex plane.
    Complex y;
    if (Complex.Abs(U) < 0.00001f)
    {
        y = -(5f / 6f) * a + U - Complex.Pow(Q, 1.0f / 3.0f);
    }
    else
    {
        y = -(5f / 6f) * a + U - P / (3 * U);
    }

    Complex W = Complex.Sqrt(a + 2 * y);
    double X = -B / (4 * A); // B and A are always real, so X is always real.
    Complex Y = 3 * a + 2 * y;
    Complex Z = (2 * b) / W;

    // The four possible roots of the quartic equation are t1, t2, t3, and t4...
    Complex t1 = (X + (W + Complex.Sqrt(-(Y + Z))) / 2);
    Complex t2 = (X + (W - Complex.Sqrt(-(Y + Z))) / 2);
    Complex t3 = (X + (-W + Complex.Sqrt(-(Y - Z))) / 2);
    Complex t4 = (X + (-W - Complex.Sqrt(-(Y - Z))) / 2);

    // Only add the real roots.
    List<double> result = new List<double>();
    if (!Complex.IsNaN(t1) && Complex.IsReal(t1)) { result.Add(t1.Re); }
    if (!Complex.IsNaN(t2) && Complex.IsReal(t2)) { result.Add(t2.Re); }
    if (!Complex.IsNaN(t3) && Complex.IsReal(t3)) { result.Add(t3.Re); }
    if (!Complex.IsNaN(t4) && Complex.IsReal(t4)) { result.Add(t4.Re); }
    root = result.ToArray();
    return root;
}

Honestly, however, I'm going to look at the solution when B == 0, since that coefficient pertains to my current velocity * my constant acceleration. Since the autopilot will slow down to a near-halt if the angle between my forward vector and the intercept vector is greater than ~0.5 degrees, it will give an incorrect answer, which I'm currently negating by just capping the coefficient to 2 * acceleration if my speed is greater than 0. This could also be a reason for the inaccuracy of the intercept vector.

I've spent the last two weeks learning and writing the code for several different variants of the quartic formula, only to realize that the problem wasn't in the formula, but rather in the coefficients I was passing to the method... but I'm very good at factoring 4th degree polynomials by hand, now
--- Automerge ---
Pertaining to my input coefficients, if you're willing and able, perhaps you'd take a look at how I'm computing the A, B, C, D, and E and let me know if you get the same values that I do? I won't post it unless you or someone else that is good with multivariable factoring accepts, as the process is a bit lengthy to write out.

 

The issue is that you are plugging a negative number into the Math.Sqrt function which will always return NaN.
Rather than this:

Code:

ComplexNumber R = new ComplexNumber(1, Math.Sqrt(Math.Pow((-P / 3), 3)));

You should do:

Code:

ComplexNumber R = new ComplexNumber(1, Math.Sqrt(Math.Pow(Math.Sign((-P / 3), 3))));

Also it is worth noting that doing A*A*A is faster than Math.Pow(A, 3). Since you are doing so many computations, this will decrease your performance overhead.

Alternatively, you could make a static class for your square roots, and (with some adjustment) I believe it will work for you:

Code:

public static class ComplexMath
{
    public static ComplexNumber Sqrt(double num)
    {
        if (num < 0)
            return new ComplexNumber(0, Math.Sqrt(-num));
        else
            return new ComplexNumber(Math.Sqrt(num), 0);
    }
}

 

Thanks for replying! But I am/was under the impression that any calculation performed using the type Complex will store the value regardless of the output. ie if the calculation is taking the square root of a negative number (which returns NaN as any other number type), the Complex takes the sqrt of the absolute value of the number, and then creates two different numbers, number.Real and number.Imaginary. The real part is the absolute value, and the imaginary part is the absolute value with an appended "i". Later in the code, if you end up multiplying two imaginary numbers together, since i squared = -1, it just multiplies the two real parts together and then multiplies that by -1.

I also think I may know why I was getting a bunch of errors after changing the type to Complex - in the complex class above, it only has operators for complex with complex (complex * complex, + complex, etc) where my code was mostly doubles. I will need to change all the types of my numbers to complex for those that are used with the complex numbers created (R, cbrtR, phi, and the Roots).

As I understand the Math.Sign, it returns either 1 or -1, depending on the outcome of the calculation. In your example, would it not return NaN still? As the sign of (-p/3)^3 is -1, it would still take the square root of -1, and return NaN.

Lastly, thank you for the class snippet and the bit about overhead! I created a class for taking cuberoots already, but didn't think anything about the overhead, just about how many times I had to write out math.pow(x, (double)1/3). I will also go through and adjust the squaring of numbers rather than using pow for those

 

I don't see Complex.Sqrt() defined anywhere so I can't check on that.

Yeah I messed that up, that should be Math.Abs() (Was typing on mobile lol)

 

No worries

Yeah, the complex class above was one I pulled off of a random site, but I realized it was more of a generic "How to" demo class than what I needed. Here is the class I am using currently. After looking at it, it seems that it does some funky old, dead dude calculation and returns a real number.

Code:

/// <summary>
/// A simple complex number class.
/// Source for Quartic solver + Complex class from:
/// https://code.google.com/p/xna-circle-collision-detection/
/// </summary>
public class Complex
{
    public double Re;
    public double Im;

    public Complex(double re, double im)
    {
        Re = re;
        Im = im;
    }

    public static Complex operator +(Complex a, Complex b)
    {
        Complex result = new Complex(
            a.Re + b.Re,
            a.Im + b.Im);
        return result;
    }

    public static Complex operator -(Complex a, Complex b)
    {
        Complex result = new Complex(
            a.Re - b.Re,
            a.Im - b.Im);
        return result;
    }

    public static Complex operator *(Complex a, Complex b)
    {
        Complex result = new Complex(
            a.Re * b.Re - a.Im * b.Im,
            a.Im * b.Re + b.Im * a.Re);
        return result;
    }

    public static Complex operator /(Complex a, Complex b)
    {
        Complex result = new Complex(
            (a.Re * b.Re + a.Im * b.Im) / (b.Re * b.Re + b.Im * b.Im),
            (a.Im * b.Re - a.Re * b.Im) / (b.Re * b.Re + b.Im * b.Im));
        return result;
    }

    public static Complex operator *(Complex a, double b)
    {
        Complex B = new Complex(b, 0);
        return a * B;
    }

    public static Complex operator /(Complex a, double b)
    {
        Complex B = new Complex(b, 0);
        return a / B;
    }

    public static Complex operator +(Complex a, double b)
    {
        Complex B = new Complex(b, 0);
        return a + B;
    }

    public static Complex operator -(Complex a, double b)
    {
        Complex B = new Complex(b, 0);
        return a - B;
    }


    public static Complex operator *(double a, Complex b)
    {
        Complex A = new Complex(a, 0);
        return A * b;
    }

    public static Complex operator /(double a, Complex b)
    {
        Complex A = new Complex(a, 0);
        return A / b;
    }

    public static Complex operator +(double a, Complex b)
    {
        Complex A = new Complex(a, 0);
        return A + b;
    }

    public static Complex operator -(double a, Complex b)
    {
        Complex A = new Complex(a, 0);
        return A - b;
    }

    public static Complex Sqrt(Complex a)
    {
        return Pow(a, 0.5f);
    }

    public static Complex Sqrt(double a)
    {
        return Pow(new Complex(a, 0), 0.5f);
    }

    public static Complex Pow(Complex a, double b)
    {
        Complex result;

        if (IsNaN(a) == false)
        {
            // De Moivre's Theorem:  (r*(cos(x) + i*sin(x))^n = r^n*(cos(x*n) + i*sin(x*n))

            // First convert the complex number to polar form...
            double r = (double)Math.Sqrt(a.Re * a.Re + a.Im * a.Im);
            double x = (double)Math.Atan2(a.Im, a.Re);

            // And then rebuild the components according to De Moivre... 
            double cos = (double)(Math.Pow(r, b) * Math.Cos(x * b));
            double sin = (double)(Math.Pow(r, b) * Math.Sin(x * b));
            result = new Complex(cos, sin);
        }
        else
        {
            result = new Complex(a.Re, a.Im);
        }
        return result;
    }

    public static Complex Pow(double a, double b)
    {
        return Pow(new Complex(a, 0), b);
    }

    public static Complex operator -(Complex a)
    {
        return new Complex(-a.Re, -a.Im);
    }

    public static double Abs(Complex a)
    {
        return (double)Math.Sqrt(a.Im * a.Im + a.Re * a.Re);
    }

    public static bool IsNaN(Complex a)
    {
        return double.IsNaN(a.Re) | double.IsNaN(a.Im);
    }

    public static bool IsReal(Complex a)
    {
        if (Math.Abs(a.Im) < 0.0000001)
            return true;
        return false;
    }
}

 

In the Complex Pow method:

Code:

else //implied if a is NaN
{
	result = new Complex(a.Re, a.Im);
}

Does not actually "fix" if you get a NaN for an entry which will end up propogating NaNs all over the place.

Also if you are just square rooting negative numbers, you can vastly simplify the square root like I showed above to use less trig(which isn't the best for performance).

Or heck... You could just get rid of all that complex number class mumbo jumbo if you wanted and just use a regular method like so:

Code:

//Returns true if imaginary
public bool SafeSqrt(double num, out double root)
{
	if (num < 0)
	{
		root = Math.Sqrt(-num);
		return true;
	}
	else
	{
		root = Math.Sqrt(num);
		return false;
	}
}

//Example use
void Main()
{
	double myRoot;
	if(SafeSqrt(-4, out myRoot)) //if imaginary
	{
    		Echo($"Imaginary: {myRoot}");
	}
	else
	{
    		Echo($"Real: {myRoot}");
	}	

}

 

You are absolutely right, though I haven't seen this case yet - I'm not exactly sure why they did it that way.

While this seems intuitive, it is actually not effective in solving the quartic. Here's why:

Remember, the square root of a negative number is an imaginary number, and by mathematical law, the square of any imaginary number is a real, negative number. The quad of an imaginary number is a POSITIVE number (think -1 * -1 = 1). Take this snippet, for example:

Code:

Else Statement values:
Unity = (0, 0.866025403784439)
sqrtD1_D0 = (308374000954.736, 0)

q_1 = (5290.07497054835, 0) // First cube root of unity = q * 1
q_2 = (-2645.03748527418, 4581.33931241909) // Second root of unity = q * (-1/2 + 1/2(sqrt(-3))
q_3 = (-2645.03748527418, -4581.33931241909) // Third root of unity = q * (-1/2 - 1/2(sqrt(-3))

s_1 = (3.81376950505458, 0)
s_2 = (10.3193545057878, 9.3541465614563)
s_3 = (10.3193545057878, -9.3541465614563)
S value used = (3.81376950505458, 0).

sqrt value = (2.29102946274903E-15, 37.4165862458253)
Root1 = (-3.81376950505458, 18.7082931229127)
Root2 = (-3.81376950505458, -18.7082931229127)
Root3 = (3.81376950505458, 18.7082931229127)
Root4 = (3.81376950505458, -18.7082931229127)
All roots were NaN...

This is the output from the quartic solver code. The values of q are the three cube roots of q. Note that q_1 is a real, positive number. The values of s are the calculations using the three q values (S cannot = 0, else you end up dividing by zero for the sqrt value). In this example, I forced it to use ONLY real, positive numbers for the calculation of S. This results in a real, positive number for S. However, take a loot at the resulting roots of the equation. They are ALL imaginary numbers (the second part denotes the imaginary part of the number).

Now take a look at this snippet:

Code:

Else Statement values:
Unity = (0, 0.866025403784439)
sqrtD1_D0 = (471358940742.067, 0)

q_1 = (6182.42082866452, 0)
q_2 = (-3091.21041433226, 5354.13349450952)
q_3 = (-3091.21041433226, -5354.13349450952)

s_1 = (4.08582196986483, 0)
s_2 = (10.9402847467361, 10.1655081565592)
s_3 = (10.9402847467361, -10.1655081565592)
S value used = (10.9402847467361, 10.1655081565592).

sqrt value = (13.7089255537425, -20.3310163131185)
Root1 = (-4.08582196986486, -20.3310163131185)
Root2 = (-17.7947475236073, 5.32907051820075E-15)
Root3 = (17.7947475236073, -5.32907051820075E-15)
Root4 = (4.08582196986486, 20.3310163131185)

Target Vector: {X:-0.892936441443875 Y:0.402966664118418 Z:0.200704706349462}
Forward Vector: {X:-0.898184788190721 Y:0.399337651327554 Z:0.183830156652771}
Dot Prod: 0.999837267347561
Angle Diff: 1.03366752350678
Root: 17.7947475236073
Time Before Delta: 17.7947475236073
Delta: 0.0469848874321262
Total Time: 19.8417324110395

This time, I force it to use s_2 (which is imaginary), and as you will see, I end up with one real root (even though it isn't purely real, there are 14 zeros after the decimal for imaginary part.. close enough ). This results in a time to target of ~18 seconds, which got me within 30 meters of the target this time (and now that I think of it, the fact that isn't purely real may be why I'm not hitting the target... just gave me some inspiration to tweak some calculations!).

The point being, I NEED the complex numbers in order to find real roots. Though, I did find a blog post where a developer talks about how he has written the code to find ballistic trajectory for every game he's worked on, and he shows all his math. I noticed in his quartic calculations that he splits the coordinates into three separate calculations, for x, y and z, first, and then combines them into one giant equation. I plan to do this today, and see if it gets me any closer. More to follow

 

Lol, that's all you had to say

That is zero. That is a floating point error likely due to trig and the power function. Use an epsilon of like 1e-9 to check if a number has significant value. As a learning excersise, I decided to dig up my old calculus notes and make a ComplexNumber struct, I went ahead and implemented epsilon in my Equal, IsReal, IsImaginary, and IsComplex functions. I also made some cheaper implementations of integer powers and Squaring complex numbers:


Now lets attack the issue at hand... Interception.
I've done a bunch of work with guided missiles and intercept algorithms, I'll tell you this: There are much simpler algorithms for calculating intercepts that work sufficiently well in this game. A simple quadratic can yield very good results for target interception:

Another issue I see is inaccuracy caused when the target is moving inertially.

Code:

if (Math.Abs(A) < 0.1)
{
	if (A == 0) { A = 0.1f; }
	else { A = Math.Sign(A) * 0.1f; }
}

This will undoubtedly cause some issues with regards to the accuracy of the computation (unless you switch to cubic here).

One of the main issues you'll encounter is SE's speed cap. That will throw so many complex computations out the window unless you code lots of safeguards into your solvers to compensate for this speed cap. Otherwise, you are predicting a target position using its current acceleration without considering that the target will stop accelerating once at the speed cap. That will vastly change the prediction. Also, some strange stuff happens once you hit the speed cap with respect to the lateral velocity of a target. That solution is not analytical, you have to solve the speed concatenation iteratively. Now that would require hard coding a speed limit since we can't fetch that from the game itself... and hard coding a speed limit makes the code less modular... etc etc.

 

Sorry lol I should have known you'd be familiar, but figured I'd explain just in case Maybe someone else will benefit!

That makes sense that it is just an error. The code I have does include the epsilon when it looks for which root to use. I was not tracking that it was a floating point error, however, so I just thought it wasn't QUITE a real solution.

So, I tried the quadratic first.. but could never get it to be accurate enough for point intercept. I feel the reason being that I'm using acceleration instead of a constant velocity, so I have to square t^2, ending up with 0.5a^2*t^4, which is what led me on my journey through Quartics. Below is how I set up my equation.

Code:

/* Variables:
*
* O = My Position
* A = Target Position
* B = Intercept Point
* C = Projected point of V on D
* D = Distance (A - B)
* V = Target's velocity vector
* vNorm = Target's direction (toward B)
* U = My velocity vector (toward B)
* Sc = My Constant Acceleration
* V_pt = Vt component Parallel to d (V projected onto d)
* V_ot = Vt component ORTH to d
*
*
*                                       B
*                                  /   /:
*         ut + 0.5*Sc(t)^2    /       / :
*                        /         vt/  :v_o*t
*                   /               /   :
*              /                   /e   :
*         O-----------------------A.....C
*                    D             v_p*t
*
*
*  Uses Kinematic Equation for constant acceleration: s(t) = si + sv(t) + 0.5a(t)^2
* 
*  Uses Pythagoras Theorem and the Cubic and Quartic Equations
* 
*  (OC)^2 + (CB)^2 = (OB)^2, where OC, CB, and OB are vectors in the above figure
* 
*/

I've seen a few scripts for guided missiles, and I've looked at the code for a couple of them. In them, I noticed a LOT of matrix calculations and reference values - I am not that adept with c# to implement such just yet.

In your video, what is it that is tracking the ship? Clearly not a remote by the way it drifts. I was trying to attempt such using thrust override and just having the remote utilize the gyroscope to rotate toward the intercept point - but the remote is kind of all or nothing in its current state.. Perhaps I'll have to look at how to code the rotation via gyro override, see if I can't figure that out.

The code I've written includes cases that will drop all the way down to quadratic if certain coefficients are zero.

Question about this: What is the cap for Enemy AI? I'm guessing it's a percentage of the speed cap. I've noticed enemy ships flying around 40 m/s, but do they accelerate past that? Also, in your script, do you hard code the speed for the missile? Or is it calculated using s = vt + 0.5at^2? I'm not using any acceleration for the enemy, but I am for my drone, since the remote likes to stop before turning toward the intercept point.

 

Cracks knuckles
Now we are speaking my language!

Here is some example intercept geometry for my quadratic method:

Variable definitions:
Vs: Target velocity - In retrospect I can't remember why in the hell I used "s" as a subscript lol
Vp: Projectile velocity - I use my current missile speed OR if the projectile is ballistic, I use the bullet's muzzle velocity
H : Direct heading - Line from you to the target
t : Time to intercept - In seconds

Implementation:


This system uses the law of cosines to quickly solve for the intercept point. When you have a missile that can adjust its heading, this is more than sufficient for most applications.

Most of what you are seeing is how people are extracting their desired Pitch, Yaw, and (sometimes) Roll angles. Then some of those matrix things are rotating vectors to make them body relative instead of world relative in order to do proper angular control.

The "cap" is the speed cap in the game. Enemy AI is really just a script controlling the behavior of the enemies. How fast they will go is up to the script writer. The regular cargo ships just spawn adrift at constant velocity.

Nope. I use the speed as a variable in my quadratic formula. It isn't perfect by any means, but it is solid and robust enough to get good hit percentages.

Ok... I highly recommend that you don't use autopilot for interceptions. Autopilot is sluggish and slow to respond. For my drones/missiles, all of the control is done manually using gyro and thruster overrides. If you make it yourself, you can custom tailor the response to your needs

 

This is essentially what I did in the beginning, but how do you account for acceleration as a function of velocity? My thoughts were that, since the drone starts from velocity 0 (or close to it), 1/2*acceleration*t^2 would have to be in the mix. So my legs were (using your diagram):

H - same, target position (M) - my position (P): M-P
Vrt - same, target position (M) + Vr*t: M + Vrt

Vpt - different: My position (P) + my velocity (Vp)t + 1/2 acceleration * t^2

So I get P + Vp*t + 0.5*a*t^2, and when I square that leg...

(0.25*a^2)t^4 + (Vp*a)t^3 + (Pa + Vp^2)t^2 + (2PVp)t + P^2

Set all that to zero:

(0.25*a^2)t^4 + (Vp*a)t^3 + (Pa + Vp^2)t^2 + (2Pvp - 2MVr)t + (-2M^2 + 2MP) = 0

and there's my quartic...

That is a little different from what I used. I was not using the -2ab*cosTheta. Perhaps something I overlooked for sure.

How are you using just the speed variable without it being constant? Or is that simply because the missile continuously updates its trajectory?

Also, how is updating with the target position when the target changes direction? Maybe I'm not "Locking" the target somehow (or don't know how to). I assume raycast to begin with, but raycast can only do so much, and if it misses the target on the next tick, what then?

EDIT: After closer review of your code, I see what it's doing. When missile velocity is 0, it just starts heading toward the target's current position, and continues to update the trajectory as the speed increases. Just have to ensure that the missile can keep track of the target.

You aren't kidding lol currently, I have it calculate the angle between forward and intercept vectors, and divide that by it's average rotational velocity. Then add that delta to "t" for the final calculation. And I've been playing with adding a second or two to account for the "wobble" as the remote settles in on a waypoint (seems to be ~ 2 seconds give or take).

 

You don't! You simply correct by updating your intercept frequently.

Yep! This is necessary to hit a maneuvering target.

Raycast can do a ton of things! Also, I have a pretty robust lock on system, but if you lose lock... you are SOL. However, I've not had any issues keeping track of maneuvering fighter sized vessels and ruining their days.

I'm not entirely sure if you can because the direction of vector3 is not fully defined. The orientations of vector1 and vector2 are clearly defined, but vector2's magnitude is variable. This will cause the end point of vector3 to also be variable since vector addition forms a triangle. That gives you two variables you need to solve for: heading of the intercept vector and length. I've also never seen a "time" vector, interesting concept. So... try it out and report back how it works, I honestly have no clue haha

Don't do dat. You should return the direct heading vector if no conditions are satisfied so that you are at least traveling towards the target. The way you have this set up, if an intercept fails, the vessel will pilot to the world origin.

 

Can you give me some pointers? Or a snippet? I can't seem to get mine to update anything lol here's what I'm using:

Code:

public void UpdateTargetInfo()
{
    foreach (var cam in allCams)
    {
        if (cam.CanScan(nextScanPoint))
        {
            info = new MyDetectedEntityInfo();
            info = cam.Raycast(nextScanPoint, (last_TargetPosition + (last_TargetDirection * last_TargetSpeed * (run * lastRunClock))));
            if (!info.IsEmpty())
            {
                if (info.EntityId == targets[0].EntityId)
                {
                    targetPosition = info.Position;
                    targetSpeed = info.Velocity.Length();
                    targetDirection = Vector3D.Normalize(info.Velocity);
                    targetRadius = BoundingSphereD.CreateFromBoundingBox(info.BoundingBox).Radius;
                    distToTarget = Vector3D.Distance(info.Position, currentPosition);

                    targetUpdated = true;

                    break;
                }
            }
        }
    }
    return;
}

It was worth a shot, but I removed it from the post because it was WILDLY inaccurate. In theory, it should work... if you solve the quadratic for t, then you know how vector 2 varies with time, and thus its magnitude. And if you know how its magnitude varies with time, then you can find vector 3. I honestly think the hiccup is in trying to take the square root of the vector, as the components may very well be negative, and the sign method likely isn't cutting it.

Yep, I realized that very quickly during testing. Ended up switching to the Target-Me vector.

 

That snippit looks sound to me. Where are the issues occuring? Have you tried using debug draw to see a visual representation of the raycasts?

Conceptually, time should not vary for the components of the velocity vector. That means that the object isn't truly traveling in the direction of the velocity vector you are multiplying by. It means that each components of the velocity has its... own time scale... good lord my brain hurts lol

Also solving for t suggests that you know the orientation (relative magnitudes of your desired steer vector) of vector3 (which we don't). We are essentially trying to concurrently solve 2 parameters with only one equation.

 

I really can't say where the issue lies if I didn't botch the method. I get a hit on the target, compute the future location and start to travel there, but on the next tick I get the output "Unable to locate target" which lets me know that the raycast was unable to find it on the next run...

But this is a first heard about debug draw... do tell?

EDIT: Quick google search and I found the debug menu! Never used it before. What are the spherical objects ?

Also, this is using the Quartic solver - the time is off by varying amounts depending on what my position is relative to the target. But usually I end up somewhere along the line it is flying, and usually right behind where it is.

Annnd.. looks like the raycasts stop instead of all aim for the target location.. back to the script I go!

EDIT: Figured out what the spheres were - the detected entities of the raycast, duh. Also, I noticed that per your theta calculation, theta was a constant 179.something degrees. So I switched to:

Code:

double theta = 180 * Math.Acos(Vector3D.Dot(Vector3D.Normalize(mPos - tPos), Vector3D.Normalize(tVel))) / PI;

Also, I realized I had to invert the distance vector, since we want the angle from the target's point of view, as it looks at the intercept vector.

Seems that the farther the target is from me, the greater the root is, however the distance to the intercept point doesn't take that long to get there. I have to be missing a calculation somewhere.. I just got a 25 second root, but flew 7 seconds to the point where the target will be in 25...

I think I just realized what I'm not doing while writing this lol I need to add the two vectors together to make the third. Gah I feel dumb sometimes lol

EDIT 2: Well now it's way off again.. if he is flying toward me, I fly past the intercept point, but also 75-100 meters behind him. And if I am chasing him down, I fly ~ 150 meters behind him, and still past his line of travel. It's like his velocity vector is off or something.

lol I've been swimming in physics and multivariable algorithms for a month straight trying to figure this out. Mine might explode

 

Uhh... my theta is not constant.

Also, you don't need theta at all since the dot product bypasses the need for a direct computation of theta. I didn't simplify in that page of my notes, but I did in the code snippit.

Here is the page with the simplification as well as a bunch of other doodles:

D: Keep all angles in terms of radians. The degrees in my notes were for readability as degrees are more intuitive to me. Everything in SE uses radians

 

Whip you've been awesome, and I thank you for taking the time to help me. I'm so dang frustrated with this it's not funny lol. What could I be doing wrong? Would you be willing to take a look at my script and/or save to see if you can find my error(s)?

I'm giving up for the night lol I thought I was just about there using your law of cosines. Maybe tomorrow!
--- Automerge ---
Ugh! I found the errors - two of them. One was a tab-complete oversight, where I input targetDirection - currentPosition instead of targetPosition - currentPosition. The other was the dot product we were discussing. I dotted the displacement vector with the target velocity vector, when it should have been the other way around.

EDIT: Nope, finding that it doesn't matter which way the dot goes, so I guess it was just the one that was messing everything up.

So, the point intercept seems spot on, with one hitch - the drone slows down on approach. I tried to create a vector between the intercept point and my position, and then add a portion of that vector to the intercept vector so I end up flying past the point at full speed. This is what I'm using for it, maybe I'm mucking it up somehow:

Code:

if (t > 0)
{
    nextScanDist = Vector3D.Distance(currentPosition, ((tPos + tVel * (run * lastRunClock)))) * 1.2f;
    var intVector = (targetPosition + tVel * t) - currentPosition;
    interceptVector = (targetPosition + tVel * t) + (intVector * 0.2f);

}

Is there another way to push the point out further but along the same vector?

 

Normalize the vector then multiply by desired length

 

Yes, exactly what I needed, thank you.

Now I just need to fine tune the time of rotation for the remote (until I get the gyro override situation figured out) as it is still off randomly by + or - 1 or 2 seconds depending on the displacement vector and angle between forward and intercept vectors at time of detection.

 

I blame the crappy autopilot for that.

When you get around to your own gyro control, look up PID control. It will help you design a fine tuned angular controller. Here is a class I made to help with that:


You simply plug in your angular error and fine tune your gains and it will compute the PID controller return

 

I was looking this up the other day after reading a thread about it. A bit beyond my ability at the moment, but I'll get there. Thank you for the class, tho! Creating classes isn't my strongsuit, currently. I learned everything I know from guys/gals like you and a whole lot of YouTube.

But, I start fundamentals of programming and programming methodologies tomorrow! I've fallen in love with coding due to this game lol time for formal training

 

Good luck man! I also self-taught programming so I could use it in this game haha. Never knew that playing games would give me a helpful life skill

 

@Whiplash141 I'm back! I've been researching how to rotate to target via gyro override, and I'm starting to understand the whole matrix bit. In my research, I came across an awesome video on Quaternions, and I think this is what I want to use. I saw in another thread you commented on how to go about calculating the angle between the vectors, and also get the axis of rotation by taking the cross of the forward and displacement (to target) vectors. But I realized in the video that the quaternion comes with the angle and axis needed to rotate to a given direction. My question to you, should you choose to answer it:

How do I go about taking the quaternion I have to make the ship rotate about the axis defined in the quaternion?

Code:

var fwdVec = remote.WorldMatrix.Forward;
var displacementVec = target.Position - remote.GetPosition();

QuaternionD q_Rot = QuaternionD.CreateFromTwoVectors(fwdVec, displacementVec); // As I understand it, there is no need to normalize the vectors as the CreateFrom automatically returns a unit vector, correct?

// Now the Quaternion has <w, x, y, z> where w = cos(theta/2) and x,y,z define the axis of rotation

What's the next step? Can I make the AoR given by the quaternion be an axis of rotation for the ship? ie. if the AoR is not "up" or "right", how do I shift the orientation of the gyroscopes axes to match the AoR?

EDIT: As I sit here thinking about it, I'm wondering if the above becomes a "roll" and "yaw" scenario instead of a "pitch" and "yaw"?

Or is it simpler still just to find the angle using the projection method you described in the linked thread?
--- Automerge ---
Oh, and also! I figured out why I wasn't tracking the target with raycasts - I was using a combination of two overloads for the raycast method. I finally went back and read the parameters of the method more carefully, and after slapping myself in the face a few times... 100% target lock!

That said... it's REALLY easy to see how jacked using the remote control is now, since it is recalculating the intercept vector every 60th of a second.. the remote never reaches a speed above 30 m/s before it slows down to readjust heading lol. Completely useless.

One last thing I remembered. When the drone and the target are facing each other (playing chicken), the intercept vector will randomly switch 180°, which is annoying as heck. It's like the algorithm doesn't realize that the shortest path is a head-on collision course. Any idea why this might be?

 

No idea. Quaternions are a black magic. Lol.

I know it works and how it works... so for me, yes.
--- Automerge ---

Which algorithm?

 

lol but they look so pretty in theory! I'm still plugging away at figuring out how to finalize the process with them. The videos were geared toward physics engines with "Just plug this thing into your simulator and out pops an animated rotation" - but I need to know how that animation happens

This one - it's a bit long sorry!

My guess is that a sign isn't being taken into account with a calculation somewhere.

 

Issue is that I like knowing how things work before I use them. Quats are strange to me. I prefer math that is easy to visualize, perform, and explain. Also, we only are capable of rotations about principal axes in Space Engineers, so a resultant rotation that is off axis is... problematic

This may help: https://www.euclideanspace.com/maths/geometry/rotations/conversions/angleToEuler/index.htm

Or you can just project the axis angle onto the principal axes and pretend that does the same thing since it may be good enough

Lol still using that quartic?

I can't help you with that math

(Quadratics do decently well )

 

I am the same way, but then I found this video while doing some research. They are actually fairly straight forward and simple with the way he explains it. If you've a spare 40 minutes, you may find it interesting!

Yep! lol one might say that this equation started out as something to help me learn more about programming and game mechanics all at the same time. But after MONTHS of research learning about the math, I refuse to throw it all away for a simple quadratic!

Though, I do not doubt that the quadratic would do just fine - once I have the rotation stuff worked out.

Thank you for the link! Very nice info on there

 

You misspelled "black magic"
And quats are nowhere near as intuitive as nice and simple azimuth elevation angles to me lol

Fair, but understand that you may spend more months of headache getting that to work when you could be a peasant like me and just use simple approximations