#!/usr/bin/perl
# 
# Use trial-and-error to find the elevation angle theta of an artillery piece
#    which will cause the shell to travel a maximum distance
#
# MWR 4/26/2018

$debug = 0;

# convert degrees to radians
$DEGTORAD = (3.14159/180.0);

# gravitational acceleration at Earth's surface (m/s)
$g = 9.80;
# height of gun above the target (m)
$H = 1300.0;
# muzzle velocity of shell (m/s)
$v = 800.0;




$start_theta_deg = 43.0;
$end_theta_deg = 47.0;
$d_theta_deg = 0.01;

for ($theta_deg = $start_theta_deg; $theta_deg <= $end_theta_deg;
         $theta_deg += $d_theta_deg) {

  $theta_rad = $theta_deg * $DEGTORAD;

  # first, we figure out how much time the shell must spend in air
  #    it involves solving a quadratic equation
  $rad = sqrt(  ($v*sin($theta_rad))**2 + 2*$g*$H );
  $top_plus = 0.0 - $v*sin($theta_rad) + $rad;
  $top_minus = 0.0 - $v*sin($theta_rad) - $rad;
  $bot = 0.0 - $g;

  $plus_time = $top_plus / $bot;
  $minus_time = $top_minus / $bot;

  if (($minus_time > 0) && ($plus_time < 0)) {
     # this is what we expect -- carry on with the plus_time
     $t_air = $minus_time;
  }
  else {
     # this is what unexpected -- print results and quit
     printf " uh-oh   plus_time %8.2f  minus_time %8.2f ?? \n",
              $plus_time, $minus_time;
     exit(1);
  }

  # compute horizontal distance travelled
  $dist = $v*cos($theta_rad)*$t_air;


  printf " theta_deg %8.2f  t_air %9.3f   dist %9.3f   \n", 
          $theta_deg, $t_air, $dist;
}


exit 0;









# we end up with this equation:
#
#                      sin(theta)                 L^2        1
#    0  =  H  +  L * [ --------- ]  -  0.5 * g * ----- * ----------------
#                      cos(theta)                 v^2     [cos(theta)]^2
#
# and we want to find the value of angle "theta" which satisfies 
#    this equation.
#

$start_theta_deg = 0.0;
$end_theta_deg = 89.0;
$d_theta_deg = 1.0;

for ($theta_deg = $start_theta_deg; $theta_deg <= $end_theta_deg;
         $theta_deg += $d_theta_deg) {

  $theta_rad = $theta_deg * $DEGTORAD;

  $x1 = $H + $L*(sin($theta_rad)/cos($theta_rad));
  $x2 = 0.5*$g*($L*$L)/($v*$v) * (1.0 / (cos($theta_rad)*cos($theta_rad)));
  $sum = $x1 - $x2;

  printf " theta_deg %8.2f  x1 %9.3f   x2 %9.3f   sum %9.3f \n", 
          $theta_deg, $x1, $x2, $sum;
}


exit 0;