#!/usr/bin/perl
# 
# generate some data to illustrate behavior of 
#   amplitude and phase angle for different
#   driving frequencies.
#
# MWR 1/26/2019

$debug = 1;

# avoid dividing by zero
$TINY = 1.0E-10;
$BIG = 1.0/$TINY;


$mypi = 3.14159;

# mass of weight (kg)
$m = 0.10;
# force constant of spring (N/m)
$k = 7.5;
# air resistance coefficient (N/m/s)
#    this "small" b yields no visible shift between resonant and natural 
#       (max value in amp is about 12 m)
$b = 0.02;
#    this "intermediate" b yields a very small visible shift
#       (max value in amp is about 2.3 m)
#$b = 0.1;
#    this "large" b yields a small visible shift
#       (max value in amp is about 0.6 m)
#$b = 0.4;
#    this "very large" b yields a large visible shift
#       (max value in amp is about 0.33 m)
$b = 0.7;
# frequency of un-damped system (rad/s)
$omega0 = sqrt($k / $m);
# magnitude of driving force (N)
$Fmax = 2.0;

# range of frequencies to use in calculations
$n_omega = 500;
$start_omega = 0;
$end_omega = $omega0*2.5;
$d_omega = ($end_omega - $start_omega)/$n_omega;


$outfile = "./make_data.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");

printf OUTFILE "# data generated to illustrate forced harmonic motion \n";
printf OUTFILE "# mass m = $m (kg) \n";
printf OUTFILE "# spring constant  k = $k (N/m) \n";
printf OUTFILE "# air resistance coeff b = $b  (N/m/s) \n";
printf OUTFILE "# undamped frequency omega0 = $omega0  (rad/s) \n";
printf OUTFILE "# driving force Fmax = $Fmax  (N) \n";
printf OUTFILE "# \n";

for ($omega = $start_omega; $omega <= $end_omega; $omega += $d_omega) {
  # amplitude
  $bot = sqrt( ($omega0**2.0 - $omega**2.0)**2.0 + (($b**2.0)/$m**2.0)*($omega**2.0) ); 
  $amp = ($Fmax / $m) / $bot;

  # phase angle (rad)
  $bot = ( $omega0**2.0 - $omega**2.0 );
  $top = 0.0 - ( ($b/$m) * $omega );
  if (abs($bot) <= $tiny) {
    if ($bot < 0) {
      $phase = -$BIG;
    } else {
      $phase = $BIG;
    }
  } else {
    $phase = atan2($top, $bot);
  }

  printf OUTFILE " omega %12.5f  amp %12.5f  phase %12.5f \n", $omega, $amp, $phase;
 
}

# now print a second stanza which simply contains 
#   the frequency of the undamped system
printf OUTFILE "\n\n\n";
printf OUTFILE "# this is the frequency of the undamped system \n";
printf OUTFILE " omega0 %12.5f  small  -1000 \n", $omega0;
printf OUTFILE " omega0 %12.5f  zero       0 \n", $omega0;
printf OUTFILE " omega0 %12.5f  big     1000 \n", $omega0;
printf OUTFILE "\n\n\n";



close(OUTFILE);


exit 0;

