+vector solve_shotdirection(vector myorg, vector myvel, vector eorg, vector evel, float spd, float newton_style)
+{
+ vector ret;
+
+ // make origin and speed relative
+ eorg -= myorg;
+ if(newton_style)
+ evel -= myvel;
+
+ // now solve for ret, ret normalized:
+ // eorg + t * evel == t * ret * spd
+ // or, rather, solve for t:
+ // |eorg + t * evel| == t * spd
+ // eorg^2 + t^2 * evel^2 + 2 * t * (eorg * evel) == t^2 * spd^2
+ // t^2 * (evel^2 - spd^2) + t * (2 * (eorg * evel)) + eorg^2 == 0
+ vector solution = solve_quadratic(evel * evel - spd * spd, 2 * (eorg * evel), eorg * eorg);
+ // p = 2 * (eorg * evel) / (evel * evel - spd * spd)
+ // q = (eorg * eorg) / (evel * evel - spd * spd)
+ if(!solution_z) // no real solution
+ {
+ // happens if D < 0
+ // (eorg * evel)^2 < (evel^2 - spd^2) * eorg^2
+ // (eorg * evel)^2 / eorg^2 < evel^2 - spd^2
+ // spd^2 < ((evel^2 * eorg^2) - (eorg * evel)^2) / eorg^2
+ // spd^2 < evel^2 * (1 - cos^2 angle(evel, eorg))
+ // spd^2 < evel^2 * sin^2 angle(evel, eorg)
+ // spd < |evel| * sin angle(evel, eorg)
+ return '0 0 0';
+ }
+ else if(solution_x > 0)
+ {
+ // both solutions > 0: take the smaller one
+ // happens if p < 0 and q > 0
+ ret = normalize(eorg + solution_x * evel);
+ }
+ else if(solution_y > 0)
+ {
+ // one solution > 0: take the larger one
+ // happens if q < 0 or q == 0 and p < 0
+ ret = normalize(eorg + solution_y * evel);
+ }
+ else
+ {
+ // no solution > 0: reject
+ // happens if p > 0 and q >= 0
+ // 2 * (eorg * evel) / (evel * evel - spd * spd) > 0
+ // (eorg * eorg) / (evel * evel - spd * spd) >= 0
+ //
+ // |evel| >= spd
+ // eorg * evel > 0
+ //
+ // "Enemy is moving away from me at more than spd"
+ return '0 0 0';
+ }
+
+ // NOTE: we always got a solution if spd > |evel|
+
+ if(newton_style == 2)
+ ret = normalize(ret * spd + myvel);
+
+ return ret;
+}
+
+vector get_shotvelocity(vector myvel, vector mydir, float spd, float newton_style, float mi, float ma)
+{
+ if(!newton_style)
+ return spd * mydir;
+
+ if(newton_style == 2)
+ {
+ // true Newtonian projectiles with automatic aim adjustment
+ //
+ // solve: |outspeed * mydir - myvel| = spd
+ // outspeed^2 - 2 * outspeed * (mydir * myvel) + myvel^2 - spd^2 = 0
+ // outspeed = (mydir * myvel) +- sqrt((mydir * myvel)^2 - myvel^2 + spd^2)
+ // PLUS SIGN!
+ // not defined?
+ // then...
+ // myvel^2 - (mydir * myvel)^2 > spd^2
+ // velocity without mydir component > spd
+ // fire at smallest possible spd that works?
+ // |(mydir * myvel) * myvel - myvel| = spd
+
+ vector solution = solve_quadratic(1, -2 * (mydir * myvel), myvel * myvel - spd * spd);
+
+ float outspeed;
+ if(solution_z)
+ outspeed = solution_y; // the larger one
+ else
+ {
+ //outspeed = 0; // slowest possible shot
+ outspeed = solution_x; // the real part (that is, the average!)
+ //dprint("impossible shot, adjusting\n");
+ }
+
+ outspeed = bound(spd * mi, outspeed, spd * ma);
+ return mydir * outspeed;
+ }
+
+ // real Newtonian
+ return myvel + spd * mydir;
+}
+