#include <iostream.h>
#include <math.h>
#include <iomanip.h>



double deflection_ratio (double Iv, double Iw, double alpha) {

  double cos_alpha = cos(alpha), sin_alpha = sin(alpha);

  return ((Iv - Iw) * cos_alpha * sin_alpha) /
         ((Iv * cos_alpha * cos_alpha) + (Iw * sin_alpha * sin_alpha));

}


double d_ratio_d_alpha (double Iv, double Iw, double alpha) {

  double cos_alpha = cos(alpha), sin_alpha = sin(alpha),
         cos2_alpha = cos_alpha * cos_alpha,
         sin2_alpha = sin_alpha * sin_alpha,
         temp =(Iv * cos2_alpha) + (Iw * sin2_alpha);

  return ( (Iv - Iw) *
             ( (temp * (cos2_alpha - sin2_alpha)) -
               (2 * cos2_alpha * sin2_alpha * (Iw - Iv)) ) ) / (temp * temp);

}


double function_to_zero (double Iv, double Iw, double alpha) {

  double cos_alpha = cos(alpha), sin_alpha = sin(alpha),
         cos2_alpha = cos_alpha * cos_alpha,
         sin2_alpha = sin_alpha * sin_alpha,
         temp =(Iv * cos2_alpha) + (Iw * sin2_alpha);

  return (temp * (cos2_alpha - sin2_alpha)) -
         (2 * cos2_alpha * sin2_alpha * (Iw - Iv));

}

// does bisection search using the whole derivative.  not used
double min_or_max_angle_using_derivative (double Iv, double Iw) {

  double low = 0, high = M_PI / 2, guess,
         derivative_at_zero = d_ratio_d_alpha (Iv, Iw, 0),
         derivative_at_guess;

  for (;;) {

    guess = (low + high) / 2;

    if (high - low < 1e-6) {
      return guess;
    }

    derivative_at_guess = d_ratio_d_alpha (Iv, Iw, guess);

    if (derivative_at_guess == 0) { // unlikely
      return guess;
    }

    if (derivative_at_zero * derivative_at_guess > 0) {
      // sign at guess = sign at zero, we're too low
      low = guess;
    } else {
      high = guess;
    }

  }

}


double min_or_max_angle (double Iw, double Iv) {

  double low = 0, high = M_PI / 2, guess,
         function_at_guess;
  const double epsilon = 0.01 * (M_PI / 180.0); // need epsolin in radians

  for (;;) {

    guess = (low + high) / 2;

    if (high - low <= 2 * epsilon) {
      return guess;
    }

    function_at_guess = function_to_zero (Iv, Iw, guess);

    if (function_at_guess == 0) { // unlikely
      return guess;
    }

    if (function_at_guess > 0) {
      // guess is too low
      low = guess;
    } else {
      high = guess;
    }

  }

}


void do_computations (double Iv, double Iw) {

  double angle;

  cout << setprecision(5) << setiosflags(ios::showpoint | ios::fixed);

  cout << "\n    Angle         Ratio       Derivative\n"
       << "--------------------------------------------\n";

  for (angle = 0; angle <= 90; angle += 90.0 / 8) {
    cout << setw(10) << angle
         << setw(14) << deflection_ratio (Iv, Iw, (M_PI / 180.0) * angle)
         << setw(14)
         << (M_PI / 180) * d_ratio_d_alpha (Iv, Iw, (M_PI / 180.0) * angle)
         << endl;
  }

  angle = (180 / M_PI) * min_or_max_angle (Iw, Iv);
  cout << "\nThe magnitude of the ratio is greatest at "
       << angle << " degrees.\n";
  cout << "At this point the ratio is "
       << deflection_ratio (Iv, Iw, (M_PI / 180) * angle) << ".\n";

  return;

}


int main () {

  double Iw, Iv;

  for (;;) {

    cout << "\nEnter Iv and Iw (any negative or zero stops program): ";
    cin >> Iv >> Iw;

    if ((Iw <= 0) || (Iv <= 0)) {
      return 0;
    }

    if (Iw == Iv) {
      cout << "\nThere are no sideways deflections in this case.\n";
    } else {
      do_computations (Iv, Iw);
    }

  }

}