#!/usr/bin/perl
#
# Compute the distance a projectile should travel when fired
#   from a tilted gun, given a set of input measurements.
#
#   Also, use the uncertainties in those measurements to estimate
#   the uncertainty in the predicted distance travelled.
#
# MWR 1/27/2026

use POSIX;
srand();
$debug = 0;

$mypi = 3.1415926;
$DEGTORAD = ($mypi)/180.0;
# accel due to grav (cm/s^2)
$g = 980.0;

# height of table above floor (cm), and its uncertainty
$h1 = 81.0;
$d_h1 = 1.0;
# height of end of muzzle above table (cm), and its uncertainty
$h2 = 13.0;
$d_h2 = 0.5;
# distance from muzzle to edge of table (cm), and its uncertainty
$x1 = 7.0;
$d_x1 = 0.5;
# muzzle velocity of gun (cm/s), and its uncertainty
$v = 363.0;
$d_v = 8.0;
# angle of barrel above horizontal (degrees), and its uncertainty
$theta_deg = 23.0;
$d_theta_deg = 2.0;
$theta_rad = $theta_deg*$DEGTORAD;
$d_theta_rad = $d_theta_deg*$DEGTORAD;

if ($debug > 0) {
  printf " v       %9.4f  d_v    %9.4f  \n", $v, $d_v;
  printf " theta   %9.4f  d_the  %9.4f  \n", $theta_deg, $d_theta_deg;
  printf " h1      %9.4f  d_h1   %9.4f  \n", $h1, $d_h1;
  printf " h2      %9.4f  d_h2   %9.4f  \n", $h2, $d_h2;
  printf " x1      %9.4f  d_x1   %9.4f  \n", $x1, $d_x1;
  printf "\n";
}



if (0 == 1) {
  for ($i = 0; $i < 1000; $i++) {
    $x = generate_val($v, $d_v);
    printf " %8.4f \n", $x;
  }
  exit 0;
}




($x2, $d_x2) = prop_uncerts($v, $d_v, $theta_rad, $d_theta_rad,
                            $h1, $d_h1, $h2, $d_h2, $x1, $d_x1);
printf " x2        %9.4f  d_x2   %9.4f  \n", $x2, $d_x2;

$other_x2 = compute_x2($v, $theta_rad, $h1, $h2, $x1);
printf " other_x2  %9.4f  \n", $x2;


#########################
#   Enter a loop in which we generate "artificial" values
#     for all the input value,s based on their uncertainties.
#     Then, we use those artificial values, we compute
#     the displacement x2.
#
#     Keep track of all the values of x2, so we can compare
#     their distribution to the estimated range given by
#     the propagation of errors.
#
#
$num_loop = 10000;
for ($i = 0; $i < $num_loop; $i++) {

   $this_v = generate_val($v, $d_v);
   $this_theta_rad = generate_val($theta_rad, $d_theta_rad);
   $this_h1 = generate_val($h1, $d_h1);
   $this_h2 = generate_val($h2, $d_h2);
   $this_x1 = generate_val($x1, $d_x1);

   $this_x2 = compute_x2($this_v, $this_theta_rad, $this_h1, $this_h2, $this_x1);

   printf "# iter %5d  x2  %8.4f \n", $i, $this_x2;

}



exit 0;



######################################################################
# PROCEDURE: generate_val
#
# DESCRIPTION: Given a value "x" and its uncertainty, "dx",
#              generate a number from a gaussian distribution
#              with mean "x" and standard deviation "dx".
#
# USAGE:     generate_val    x  dx
#
#        where
#                   x          is the mean of the distribution
#
#                   dx         is the stdev in the distribution
#
# RETURNS:    this_x
#
#                   this_x     value from random distribution
#
sub generate_val {

  my($x, $dx, $this_x);
  
  $x = $_[0];
  $dx = $_[1];

  $this_x = $x + $dx*gaussian_rand();

  return($this_x);
}



######################################################################
# PROCEDURE: prop_uncerts
#
# DESCRIPTION: Given a set of measurements and uncertainties,
#              use propagation of errors to compute the 
#              final displacement x2 and its uncertainty.
#
# USAGE:     prop_uncerts  v dv  theta dtheta  h1 dh1  h2 dh2  x1 dx1 
#
#        where
#                   v   dv           are muzzle velocity and uncert (cm/s)
#
#                   theta dtheta         elevation angle and uncert (deg)
#
#                   h1  dh1              table height and uncert (cm)
#
#                   h2  dh2              gun height and uncert (cm)
#
#                   x1  dx1              gun-to-edge dist and uncert (cm)
#
# RETURNS:
#             (x2  dx2)
# 
#        where
#                   x2  dx2          are dist on floor table to shot (cm)
#                                           and its uncertainty
#
sub prop_uncerts {

  my($v, $d_v, $theta_rad, $d_theta_rad, $h1, $d_h1, $h2, $d_h2);
  my($x1, $d_x1);
  

  $v = $_[0];
  $d_v = $_[1];
  $theta_rad = $_[2];
  $d_theta_rad = $_[3];
  $h1 = $_[4];
  $d_h1 = $_[5];
  $h2 = $_[6];
  $d_h2 = $_[7];
  $x1 = $_[8];
  $d_x1 = $_[9];



  my($d_cos, $d_sin);
  
  # compute uncertainties (approx) in cos(theta) and sin(theta)
  $d_cos = 0.5*(cos($theta_rad + $d_theta_rad) - cos($theta_rad - $d_theta_rad));
  $d_cos = abs($d_cos);
  $d_sin = 0.5*(sin($theta_rad + $d_theta_rad) - sin($theta_rad - $d_theta_rad));
  $d_sin = abs($d_sin);
  
  if ($debug > 0) {
    printf " cos     %9.4f  d_cos  %9.4f  \n", cos($theta_rad), $d_cos;
    printf " sin     %9.4f  d_sin  %9.4f  \n", sin($theta_rad), $d_sin;
  }
  
  # compute components of initial velocity
  $vx = $v * cos($theta_rad);
  $d_vx = $vx * (  ($d_v/$v) + ($d_cos/cos($theta_rad))  );
  $vy = $v * sin($theta_rad);
  $d_vy = $vy * (  ($d_v/$v) + ($d_sin/sin($theta_rad))  );
  
  if ($debug > 0) {
    printf " vx      %9.4f  d_vx   %9.4f  \n", $vx, $d_vx;
    printf " vy      %9.4f  d_vy   %9.4f  \n", $vy, $d_vy;
  }
  
  # compute total height and its uncert
  my($H, $d_H);
  $H = $h1 + $h2;
  $d_H = $d_h1 + $d_h2;
  if ($debug > 0) {
    printf " H       %9.4f  d_H    %9.4f  \n", $H, $d_H;
  }
  
  # value of 2gH and its uncert
  my($gH2, $d_gH2);
  $gH2 = 2.0*$g*$H;
  $d_gH2 = 2.0*$g*$d_H;
  if ($debug > 0) {
    printf " gH2     %9.4f  d_gh2  %9.4f  \n", $gH2, $d_gH2;
  }
  
  
  # value of vy^2 and its uncert
  my($vysq, $d_vysq);
  $vysq = $vy*$vy;
  $d_vysq = 2.0*$vy*$d_vy;
  if ($debug > 0) {
    printf " vysq   %9.4f  d_vysq %9.4f  \n", $vysq, $d_vysq;
  }
  
  # uncert in the square root term: sqrt(vy^2 - 2gH)
  #   we define this quantity, the uncertainty in the sqrt, as "d_sqrt"
  my($top, $bot, $d_sqrt);
  $top = 0.5*($d_vysq + $d_gH2);
  $bot = sqrt($vy*$vy + 2.0*$g*$H);
  $d_sqrt = $top / $bot;
  if ($debug > 0) {
    printf " sqrt    %9.4f  d_sqrt %9.4f  \n", $bot, $d_sqrt;
  }
  
  # time in air (s), and its uncertainty
  my($tair, $d_tair);
  $tair = (  ($vy + sqrt($vy*$vy + 2.0*$g*$H)) / $g );
  $d_tair = ($d_vy + $d_sqrt) / $g;
  if ($debug > 0) {
    printf " tair    %9.4f  d_tair %9.4f  \n", $tair, $d_tair;
  }
  
  # total distance travelled (cm), and its uncert
  my($xtot, $d_xtot);
  $xtot = $tair * $vx;
  $d_xtot = $xtot * ( ($d_tair/$tair) + ($d_vx/$vx) );
  if ($debug > 0) {
    printf " xtot    %9.4f  d_xtot %9.4f  \n", $xtot, $d_xtot;
  }
  
  # distance travelled away from base of table (cm) and its uncert
  my($x2, $d_x2);
  $x2 = $xtot - $x1;
  $d_x2 = $d_xtot + $d_x1;
  if ($debug > 0) {
    printf " x2      %9.4f  d_x2   %9.4f  \n", $x2, $d_x2;
  }
  
  
  return($x2, $d_x2);
}


######################################################################
# PROCEDURE: compute_x2
#
# DESCRIPTION: Given a set of measurements (no uncertainties),
#              compute the final displacement x2.
#
# USAGE:     compute_x2  v theta  h1  h2  x1  
#
#        where
#                   v                are muzzle velocity (cm/s)
#
#                   theta                elevation angle (deg)
#
#                   h1                   table height (cm)
#
#                   h2                   gun height (cm)
#
#                   x1                   gun-to-edge dist (cm)
#
# RETURNS:
#             x2 
# 
#        where
#                   x2               is dist on floor table to shot (cm)
#
sub compute_x2 {

  my($v, $theta_rad, $h1, $h2, $x1);
  

  $v = $_[0];
  $theta_rad = $_[1];
  $h1 = $_[2];
  $h2 = $_[3];
  $x1 = $_[4];


  # compute components of initial velocity
  $vx = $v * cos($theta_rad);
  $vy = $v * sin($theta_rad);
  
  if ($debug > 1) {
    printf " vx      %9.4f  d_vx   %9.4f  \n", $vx, $d_vx;
    printf " vy      %9.4f  d_vy   %9.4f  \n", $vy, $d_vy;
  }
  
  # compute total height and its uncert
  my($H, $d_H);
  $H = $h1 + $h2;
  if ($debug > 1) {
    printf " H       %9.4f  d_H    %9.4f  \n", $H, $d_H;
  }
  
  # value of 2gH and its uncert
  my($gH2, $d_gH2);
  $gH2 = 2.0*$g*$H;
  if ($debug > 1) {
    printf " gH2     %9.4f  d_gh2  %9.4f  \n", $gH2, $d_gH2;
  }
  
  
  # value of vy^2 and its uncert
  my($vysq, $d_vysq);
  $vysq = $vy*$vy;
  if ($debug > 1) {
    printf " vysq   %9.4f  d_vysq %9.4f  \n", $vysq, $d_vysq;
  }
  
  # time in air (s), and its uncertainty
  my($tair, $d_tair);
  $tair = (  ($vy + sqrt($vy*$vy + 2.0*$g*$H)) / $g );
  if ($debug > 1) {
    printf " tair    %9.4f  d_tair %9.4f  \n", $tair, $d_tair;
  }
  
  # total distance travelled (cm), and its uncert
  my($xtot, $d_xtot);
  $xtot = $tair * $vx;
  if ($debug > 1) {
    printf " xtot    %9.4f  d_xtot %9.4f  \n", $xtot, $d_xtot;
  }
  
  # distance travelled away from base of table (cm) and its uncert
  my($x2, $d_x2);
  $x2 = $xtot - $x1;
  if ($debug > 1) {
    printf " x2      %9.4f  d_x2   %9.4f  \n", $x2, $d_x2;
  }
  
  
  return($x2);
}










######################################################################
# PROCEDURE: gaussian_rand
#
# DESCRIPTION: Returns a random number drawn from a gaussian (normal)
#              distribution with mean 0.0 and standard deviation 1.0.
#
#              Stolen from 
#
#                 http://www.unix.com.ua/orelly/perl/cookbook/ch02_11.htm
#
#
sub gaussian_rand {
    my ($u1, $u2);  # uniformly distributed random numbers
    my $w;          # variance, then a weight
    my ($g1, $g2);  # gaussian-distributed numbers

    do {
        $u1 = 2 * rand() - 1;
        $u2 = 2 * rand() - 1;
        $w = $u1*$u1 + $u2*$u2;
    } while ( $w >= 1 );

    $w = sqrt( (-2 * log($w))  / $w );
    $g2 = $u1 * $w;
    $g1 = $u2 * $w;
    # could return both if wanted
    #   return wantarray ? ($g1, $g2) : $g1;
    # but I need just one
    return($g1);
}