UnrealScript GesBio subclass

//=============================================================================
// qGESBioRifle.
//=============================================================================
class qGESBioRifle expands GESBioRifle;

// ############ Complex Numbers for solving a quartic equation #########
struct Complex
{
    var float Re, Im;
};

static final function Complex Compl( float Re, optional float Im )
{
    local Complex A;
    A.Re = Re;
    A.Im = Im;
    return A;
}

static final preoperator Complex - ( Complex A )
{
    A.Re = - A.Re;
    A.Im = - A.Im;
    return A;
}

static final preoperator Complex / ( Complex A )
{
    local float f;
    f = A.Re*A.Re + A.Im*A.Im;
    A.Re /= f;
    A.Im /= -f;
    return A;
}

static final operator(20) Complex + ( Complex A, Complex B )
{
    A.Re += B.Re;
    A.Im += B.Im;
    return A;
}

static final operator(20) Complex - ( Complex A, Complex B )
{
    A.Re -= B.Re;
    A.Im -= B.Im;
    return A;
}

static final operator(16) Complex * ( Complex A, Complex B )
{
    local Complex C;
    C.Re = A.Re*B.Re - A.Im*B.Im;
    C.Im = A.Re*B.Im + A.Im*B.Re;
    return C;
}

static final operator(16) Complex / ( Complex A, Complex B )
{
    return A * (/B);
}

static final operator(12) Complex ** ( Complex A, float Power )
{
    local float Radius, Angle;
    local Complex B;
    // polar form: (Re,Im)= Radius*(cos(Angle),sin(Angle))
    Radius = sqrt( A.Re*A.Re + A.Im*A.Im );
    if ( A.Re > 0 )
        Angle=Atan( A.Im/A.Re );
    else if ( A.Re < 0 )
    {
        if ( A.Im >= 0 )
            Angle=Atan( A.Im/A.Re ) + pi;
        else
            Angle=Atan( A.Im/A.Re ) - pi;
    }
    else // A.Re == 0
    {
        if ( A.Im > 0 )
            Angle = pi/2.;
        else if ( A.Im < 0 )
            Angle = -pi/2.;
    }
    // calculate pow
    Radius = Radius ** Power;
    Angle *= Power;
    // transform back to cartesian form
    B.Re = Radius*cos(Angle);
    B.Im = Radius*sin(Angle);
    return B;
}

static final operator(20) Complex + ( Complex A, float f ){ A.Re += f; return A; }
static final operator(20) Complex + ( float f, Complex A ){ A.Re += f; return A; }
static final operator(20) Complex - ( Complex A, float f ){ A.Re -= f; return A; }
static final operator(20) Complex - ( float f, Complex A ){ return Compl(f)-A; }
static final operator(16) Complex * ( Complex A, float f ){ A.Re *= f; A.Im *= f; return A; }
static final operator(16) Complex * ( float f, Complex A ){ A.Re *= f; A.Im *= f; return A; }
static final operator(16) Complex / ( Complex A, float f ){ A.Re /= f; A.Im /= f; return A; }
static final operator(16) Complex / ( float f, Complex A ){ return Compl(f)/A; }
static final function Complex SquareRoot(Complex A){ return A**0.5; }
static final function Complex CubeRoot(Complex A){ return A**(1./3.); }

// returns the smallest positive solution x of e*x^4 + a*x^3 + b*x^2 + c*x + d = 0;
static function float SolveQuartic(float e, float a, float b, float c, float d)
{
    local float p, q, r, s, t;
    local Complex h, u, v, z, i, j;
    local Complex x1, x2, x3, x4;
    local Complex x,y;

    if ( e != 1 )
    {
        a /= e;
        b /= e;
        c /= e;
        d /= e;
    }
    // Now we have x^4 + a*x^3 + b*x^2 + c*x + d = 0;
    // substitute x = y - a/4, the scalar before y^3 disappears, we have
    // y^4 + p*y^2 + q*y + r = 0 with
    p = -3.*a*a/8. + b;
    q = a*a*a/8. - a*b/2. + c;
    r = -3.*a*a*a*a/256. + a*a*b/16. - a*c/4. + d;
    if ( abs(q) < 0.00001 )
    {   // (y^2 + p/2)^2 = - r + p^2/4
        // y^2 = 0.5*( -p +/- sqrt(p^2-4r) )
        i = SquareRoot( Compl(p*p-4*r) );
        j = SquareRoot( 0.5*( -p + i ) );
        x1 = -a/4. + j;
        x2 = -a/4. - j;
        j = SquareRoot( 0.5*( -p - i ) );
        x3 = -a/4. + j;
        x4 = -a/4. - j;
    }
    else
    {
    // This is equivalent to (y^2 + p/2)^2 = -q*y - r + p^2/4
    // We want squares on both sides, add a variable z, such that
    // (y^2 + p/2 + z/2)^2 = y^2*z -q*y - r + p^2/4 + p*z/2 + z*z/4 is
    // (y^2 + p/2 + z/2)^2 = z*( y - 0.5*q/z )^2 . So z must fulfill
    // 0.25*q^2/z^2 = 1/z*(-r+p^2/4) + p/2 + z/4, which is the cubic equation
    // z^3 + 2*p*z^2 + (p^2-4*r)*z - q^2 = 0
    // subsitute z = w - 2*p/3, the scalar before w^2 disappears, we have
    // w^3 + 3*s*w + 2*t = 0 with
        s = -p*p/9. - 4.*r/3.;
        t = -p*p*p/27. + 4.*p*r/3. - q*q/2.;
        if ( abs(s) < 0.00001 ) // w^3 = -2*t
        {
            // z = -2.*p/3. + CubeRoot( Compl(-2.*t) );
            if ( -2*t >= 0 )
                z = - 2.*p/3. + Compl( (-2.*t)**(1./3.) );
            else
                z = - 2.*p/3. - Compl( (+2.*t)**(1./3.) );
        }
        else
        {
    // let w = u + v, then w^3 = -2*t - 3*s*w has form
    // u^3 + v^3 + 3*u*v*(u+v) = -2*t - 3*s*(u+v)
    // assuming u^3 + v^3 = -2*t and 3*u*v = -3*s
    // then we have for f = u^3 and g = v^3:
    // f + g = -2*t   ,   f*g = -s^3
    // and (f-g)^2 = (f+g)^2 - 4*f*g = 4*t^2 + 4*s^3
    // the linear system f - g = 2*sqrt(t^2+s^3) , f + g = -2*t
    // has the solutions f = -t + sqrt(t^2+s^3), g = -t - sqrt(t^2+s^3)
    // Now we can solve u,v and finally z = u + v - 2*p/3
            h = SquareRoot( Compl(t*t+s*s*s) );
            u = CubeRoot( -t + h);
            v = -s/u;            // u*v = -s
            z = u + v - 2.*p/3.; // theoretically z should be real
        }
    // Back to (y^2 + p/2 + z/2)^2 = z*( y - 0.5*q/z )^2, take the root
    // y^2 + p/2 + z/2 = µ*sqrt(z)*( y - 0.5*q/z ) with µ=+1 or µ=-1
    // y^2 - µ*i*y = -p/2 - z/2 - µ*0.5*q/i  with i =sqrt(z)
    // ( y - µ*i/2 )^2 = z/4 - p/2 - z/2 - µ*0.5*q/i), take the root
    // y - µ*i/2 = #*0.5*sqrt(-z-2*p-µ*2*q/i ) with #=+1 or #=-1
    // x = -a/4 + 0.5*( µ*i +#*sqrt(-z-2*p-µ*2*q/i ) )
        i = SquareRoot(z);
        j = SquareRoot( -z-2*p-2*q/i );
        x1 = -a/4. + 0.5*( i + j ); // µ=+1, #=+1
        x2 = -a/4. + 0.5*( i - j ); // µ=+1, #=-1
        j = SquareRoot( -z-2*p+2*q/i );
        x3 = -a/4. + 0.5*( -i + j ); // µ=-1, #=+1
        x4 = -a/4. + 0.5*( -i - j ); // µ=-1, #=-1
    }
//  log(x1.Re@x1.Im); log(x2.Re@x2.Im); log(x3.Re@x3.Im); log(x4.Re@x4.Im);
//  x = x1; y = x*x*x*x + a*x*x*x + b*x*x + c*x + d; log(y.Re@y.Im);
//  x = x2; y = x*x*x*x + a*x*x*x + b*x*x + c*x + d; log(y.Re@y.Im);
//  x = x3; y = x*x*x*x + a*x*x*x + b*x*x + c*x + d; log(y.Re@y.Im);
//  x = x4; y = x*x*x*x + a*x*x*x + b*x*x + c*x + d; log(y.Re@y.Im);
    // choose the smallest positive real root
    e = 1000000;
    if ( abs(x1.Im) < 0.00005 && x1.Re >= 0 &&  x1.Re < e )
        e =  x1.Re;
    if ( abs(x2.Im) < 0.00005 && x2.Re >= 0 &&  x2.Re < e )
        e =  x2.Re;
    if ( abs(x3.Im) < 0.00005 && x3.Re >= 0 &&  x3.Re < e )
        e =  x3.Re;
    if ( abs(x4.Im) < 0.00005 && x4.Re >= 0 &&  x4.Re < e )
        e =  x4.Re;
    if ( e == 1000000 )
        return 0;
    else
        return e;
}

// returns the smallest positive solution x of a*x² + b*x + c = 0
static function float SolveQuadratic(float a, float b, float c)
{
    local float d;

    if ( a != 1 )
    {
        b /= a;
        c /= a;
    }
    a = b*b - 4*c;
    if ( a < 0 )
        return 0;
    if ( c < 0 )
        return 0.5*( -b + sqrt(a) );
    else if ( b < 0 )
        return 0.5*( -b - sqrt(a) );
    else
        return 0;
}

// #################### AIM FUNCTION ################################

function Projectile ProjectileFire(class<projectile> ProjClass, float ProjSpeed, bool bWarn)
{
    local Vector Start, X,Y,Z;
    local float qTime;
    local Pawn qOwner;
    local Actor qTarget;

    qOwner = Pawn(Owner);
    if ( Owner.Target == none && qOwner != none )
    {
        if ( qOwner.enemy != none && !qOwner.enemy.bDeleteme && qOwner.enemy != self )
            Owner.Target = qOwner.enemy;
    }
    qTarget = Owner.Target;
    if ( qTarget == none || qOwner == none )
        return Super.ProjectileFire(ProjClass,ProjSpeed, bWarn);

    GetAxes(qOwner.ViewRotation,X,Y,Z);
    Start = Owner.Location + CalcDrawOffset() + FireOffset.X * X + FireOffset.Y * Y + FireOffset.Z * Z;
    // call own function instead of owner.AdjustToss()
    AdjustedAim = SmartAim(ProjSpeed, -0.5*Region.Zone.ZoneGravity,
                qTarget.Velocity + qTarget.Base.Velocity - vect(0,0,120),
                qTarget.Location - Start, qTime); // <----
    // PlayerPawn(qTarget).ClientMessage(qTarget.Base);
    if ( qTime <= 0 )
        return none; // target too far away or too far above
    Owner.MakeNoise(Pawn(Owner).SoundDampening);
    if ( (bWarn || FRand() < 0.4) && Pawn(Target) != None )
        Pawn(Target).WarnTarget(qOwner, ProjSpeed, Normal(qTarget.Location - Start) );
    return Spawn(ProjClass,,, Start, AdjustedAim);
}

// solves equation t*v*d = 0.5*a*t² + b*t + c, where
// d: unknown, normed (d dot d = 1) shooting direction, which will be returned as a rotator
// t: unknown scalar standing for the time
// v: scalar standing for the projetile speed
// a: vector standing for the gravity
// b: vector standing for the initial speeds, which are independent of shooting direction
// c: vector standing for the distance to your target
function rotator SmartAim(float v, vector a, vector b, vector c, optional out float t)
{
    // solution can be obtained via comparing the lenghts of both sides
    // t²*v² = 0.25*a²*t^4 + (a dot b)*t^3 + (a dot c + b²)*t² + 2*(b dot c)*t + c²
    // which is a quartic equation in t
    // if we know t, we can easily calculate d
    if ( a != vect(0,0,0) && b != vect(0,0,0) )
        t = SolveQuartic(a dot a*0.25, a dot b, a dot c + b dot b - v*v , b dot c*2., c dot c);
    else if ( a != vect(0,0,0) )
        t = sqrt( SolveQuadratic(a dot a*0.25, a dot c - v*v, c dot c) );
    else if ( b != vect(0,0,0) )
        t = SolveQuadratic(b dot b - v*v, b dot c*2., c dot c);
    else
        t = vsize(c)/v;
    return rotator(0.5*t*t*a + t*b + c);
}

// return delta to combat style
function float SuggestAttackStyle()
{
    return +0.3;
}