Hey!!  

New user here but been around the website for a while, usually just reading articles and tutorials. 

Anyways! I'm stuck and I need help with an exercise I'm doing. I am reading a book called Mathematics and Physics for Programmers (2nd edition), currently at the very end of chapter 7, exercise 7.2.

This is the only exercise I'm struggling to work out. I'm using the following formula to calculate the position of the ball at time t: (the length is the cannon's length, the velocity is the speed) 

[ (length + time * velocity)cosθ,  (length + time * velocity)sinθ - gravity * time ^ 2 / 2 ]

I am calculating the time t by using the following formula: (let u = speed)

t = u sin(θ) + sqrt(u^2 * sin(θ)^2 + 20 * length * sin(θ) / 10

I am currently hard coding the angle (θ) to 30° because I can't seem to be able to calculate the 'best angle' at which to aim the cannon so as to hit the point. The author recommends to use an approximation method to approach the solution iteratively:

Quote

You can start by creating a function that tests whether the current aim is high or low (by working out the ball’s y position when its x position is p1). This function should return 1 for high, −1 for low, and 0 for a hit or close to it. Notice first of all that if you aim the cannon directly at the target point, it will always fall short. If you use this angle as your baseline, then work up in increments until you find a firing angle that is aiming high. You can then use a simple binary approximation method to find the solution.

I just didn't understand how... 

Here is my attempt so far:

#define GRAVITY 9.81

struct Vector2D
{
    double m_X;
    double m_Y;
};

double aimCannon(double cannonLength, double muzzleSpeed, Vector2D aimPoint)
{
    double Time = 0.0;
    double AngleTheta = 0.523599; // Equivalent to 30 degrees in radians

    Time = (muzzleSpeed * sin(AngleTheta) + (sqrt(pow(muzzleSpeed, 2) * pow(sin(AngleTheta), 2) + 20 * cannonLength * sin(AngleTheta)))) / 10;

    Vector2D aimPoint;
    aimPoint.m_X = (cannonLength + (Time * muzzleSpeed)) * cos(AngleTheta);
    aimPoint.m_Y = (cannonLength + (Time * muzzleSpeed)) * sin(AngleTheta) - (GRAVITY * pow(Time, 2)) / 2;

    return 0;
}

int main()
{
    Vector2D tempAimCoords;

    double Speed = 20;
    double CannonLength = 2.0;

    tempAimCoords.m_X = 10.0f; tempAimCoords.m_Y = 15.0f;

    aimCannon(CannonLength, Speed, tempAimCoords);

    return 0;
}

My function doesn't return an angle yet.. 

I appreciate any guidance as I seem a bit lost and confused  

Thanks and sorry for the lengthy post. 

(1) Think of how long it will take the ball to reach the target, given the angle. You can do this looking just at the x coordinate, which requires solving just a linear equation. Now you can compute the y coordinate at that time. If you fired in the right angle, the answer will be the y coordinate of the target. So this process gives you an equation to solve, which you can now do by your favorite method.

(2) Another way to approach this problem is to compute the desired muzzle speed as a function of how long it takes you to hit the target. I know the muzzle speed is given, but bear with me. If you know the time to reach the target, you know how much the projectile is going to drop due to gravity. So you can just imagine the target is that much higher and then shoot in a straight line, ignoring gravity. That will give you the desired speed. Since you know what the muzzle speed is, you can solve for the time it takes to reach the target. Once you have that, use the same trick of moving the target up the desired time and firing straight.

(3) Yet another [perhaps less practical] way to think about it: Imagine the set of points that can be reached in time t with your bullet. At time 0 this set of points is a circle that describes all the possible positions for the end of the cannon. At time t the set will be a circle of radius (length_of_cannon + t / ball_speed). If there were no gravity, the sphere would be centered at the base of the cannon. Because of gravity, all the points are lower by gravity*t^2/2. So you can think of gravity as acting on this growing sphere normally. At some point this dropping sphere will engulf the target. Compute at what time that happens. Then computing the angle from the time is easy, using the linear equation given by the horizontal coordinate, or the trick of moving the target up, from (2).

I believe (1) is what most people would do, (3) is the first thing that comes to my mind --because I am weird--, and (2) is the way I would do it after thinking about this for the last 10 minutes or so. Notice how (2) and (3) don't really require trigonometry until the last step, when you express the direction of firing as an angle (using atan2), because that's the format that was asked. You could just as well return the answer as a unit vector, and then you wouldn't have angles anywhere. Just the way I like it.

I had a bit of time, so I programmed approach (2):

#include <iostream>
#include <vector>
#include <cmath>
#include <cstdlib>

double const gravity = 9.81;

struct Polynomial {
  std::vector<double> c; // coefficients in decreasing degree order
  
  size_t degree() const {
    return c.size() - 1;
  }
  
  double &operator[](size_t i) {
    return c[i];
  }
  
  double operator[](size_t i) const {
    return c[i];
  }
  
  Polynomial(std::vector<double> const &c) : c(c) {
  }
};

bool extract_root(Polynomial &polynomial, double &root) {
  size_t degree = polynomial.degree();
  
  for (int attempt = 0; attempt < 5; ++attempt) {
    // Start with a random guess
    root = -1.0 + 2.0 * double(std::rand()) / RAND_MAX;
    
    // Refine the guess several times
    double value, derivative;
    for (int step = 0; step < 30; ++step) {
      value = derivative = 0.0;
      for (size_t i = 0; i <= degree; ++i) {
	value = root * value + polynomial[i];
	if (i < degree)
	  derivative = root * derivative + polynomial[i] * (degree - i);
      }
      double old_root = root;
      root -= value / derivative; // Newton-Raphson iteration
      if (root == old_root)
	break;
    }
    
    if (std::abs(value) < 1e-6)
      goto DONE;
  }
  
  return false; // No root found
  
 DONE:;
  
  // Divide by (x - root)
  double carry = 0.0;
  for (size_t i = degree + 1; i --> 0;)
    carry = root * (polynomial[degree - i] += carry);
  polynomial.c.pop_back();
  
  return true; // Success!
}

std::vector<double> find_all_real_roots(Polynomial polynomial) {
  size_t degree = polynomial.degree();
  std::vector<double> roots;
  
  for (int i = 0; i < degree; ++i) {
    double root;
    if (extract_root(polynomial, root))
      roots.push_back(root);
    else
      break;
  }
  
  return roots;
}

double aimCannon(double cannonLength, double muzzleSpeed, double aimPoint_x, double aimPoint_y) {
  // I am going to use short names so the formulas below don't get too long
  double const l = cannonLength;
  double const s = muzzleSpeed;
  double const x = aimPoint_x;
  double const y = aimPoint_y;
  double const g = gravity;
  
  Polynomial p(std::vector<double>({ 0.25*g*g, 0.0, g*y - s*s, -2*l*s, x*x + y*y - l*l }));
  
  double t = 1e40;
  for (double r : find_all_real_roots(p)) {
    if (r > 0.0 && r < t)
      t = r;
  }
  
  if (t == 1e40)
    throw "No solution found!";
    
  return std::atan2(y + 0.5*g*t*t, x);
}

int main() {
  try {
    double const cannonLength = 1.0;
    double const muzzleSpeed = 15.0;
    double const aimPoint_x = 10.0;
    double const aimPoint_y = 5.0;
    
    std::cout << "Aiming for (" << aimPoint_x << ',' << aimPoint_y << ") at speed " << muzzleSpeed << ".\n";
    
    double const angle = aimCannon(cannonLength, muzzleSpeed, aimPoint_x, aimPoint_y);
    std::cout << "The angle is " << angle << '\n';
    
    for (double t = 0.0; ; t += .1) {
      double x = std::cos(angle) * (cannonLength + muzzleSpeed * t);
      double y = std::sin(angle) * (cannonLength + muzzleSpeed * t) - 0.5 * gravity * t * t;
      std::cout << "t=" << t << " (" << x << "," << y << ")\n";
      if (x >= aimPoint_x)
	break;
    }
  }
  catch (char const *msg) {
    std::cout << "ERROR: " << msg << '\n';
  }
}