#!/usr/bin/perl
# 
# generate some data to illustrate behavior of 
#   amplitude-as-function-of-time for a critical damping case
#
#   This version FAILS to solve the desired initial conditions
#
#     x(0) = H
#     vx(0) = 0
#
# MWR 1/30/2019

$debug = 1;

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


$mypi = 3.14159;

# this is the desired starting position at time t = 0 (m)
$start_pos = 5.0;
# mass of weight (kg)
$m = 2.0;
# force constant of spring (N/m)
$k = 32.0;
# 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 = sqrt(4.0*$k*$m);
# frequency of un-damped system (rad/s)
$omega0 = sqrt($k / $m);
# period of un-damped system
$period = (2.0*$mypi)/($omega0);
# time constant (s)
$tau = 2.0*$m/$b;
# Q factor (no units)
$Q = ($omega0*$m)/$b;
# coeffs (m)
$A = $start_pos;
$B = 0.0;

# range of time to use in calculations
$n_per_period = 500;
$n_period = 1.5;
$start_t = 0;
$end_t = $n_period*$period;
$d_t = ($period/$n_per_period);

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

printf OUTFILE "# data generated to illustrate Q \n";
printf OUTFILE "# force F = $F (N) \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 "# coeff A = $A (m) \n";
printf OUTFILE "# coeff B = $B (m) \n";
printf OUTFILE "# undamped frequency omega0 = $omega0  (rad/s) \n";
printf OUTFILE "# time constant tau = $tau  (s) \n";
printf OUTFILE "# Q = $Q  \n";
printf OUTFILE "# \n";

for ($t = $start_t; $t <= $end_t; $t += $d_t) {

  # pos
  $y = $A*exp(0.0-($t/$tau)) + $B*$t*exp(0.0-($t/$tau));

  # velocity
  $v1 = -(1.0/$tau)*$A*exp(0.0-($t/$tau));
  $v2 = $B*exp(0.0-($t/$tau));
  $v3 = -(1.0/$tau)*$B*$t*exp(0.0-($t/$tau));
  $v = $v1 + $v2 + $v3;


  printf OUTFILE " t %12.5f  pos %12.5f  v1 %12.5f  v2 %12.5f  v3 %12.5f  v %12.5f \n", 
            $t, $y, $v1, $v2, $v3, $v;
 
}


#  
#  # now print a second stanza, which shows the behavior for
#  #   a larger damping 
#  $b = 0.8;
#  
#  
#  printf OUTFILE "\n\n\n";
#  
#  printf OUTFILE "# data for an intermediate damping \n";
#  printf OUTFILE "# force F = $F (N) \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 "# phase offset  phi = $phi (rad) \n";
#  printf OUTFILE "# undamped frequency omega0 = $omega0  (rad/s) \n";
#  printf OUTFILE "# time constant tau = $tau  (s) \n";
#  printf OUTFILE "# Q = $Q  \n";
#  printf OUTFILE "# \n";
#  
#  for ($omega = $start_omega; $omega <= $end_omega; $omega += $d_omega) {
#  
#    # amplitude
#    $top = $F / $m;
#    $bot = (($omega0**2.0 - $omega**2.0)**2.0) + ((($b*$b)/($m*$m))*($omega*$omega));
#    $bot = sqrt($bot);
#  
#    $amp = $top / $bot;
#  
#    printf OUTFILE " omega %12.5f  top %12.5f  bot %12.5f  amp %12.5f \n", 
#              $omega, $top, $bot, $amp;
#   
#  }
#  
#  
#  
#  # now print a third stanza, which shows the behavior for
#  #   a very large damping 
#  $b = 2.0;
#  
#  
#  printf OUTFILE "\n\n\n";
#  
#  printf OUTFILE "# data for a large damping \n";
#  printf OUTFILE "# force F = $F (N) \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 "# phase offset  phi = $phi (rad) \n";
#  printf OUTFILE "# undamped frequency omega0 = $omega0  (rad/s) \n";
#  printf OUTFILE "# time constant tau = $tau  (s) \n";
#  printf OUTFILE "# Q = $Q  \n";
#  printf OUTFILE "# \n";
#  
#  for ($omega = $start_omega; $omega <= $end_omega; $omega += $d_omega) {
#  
#    # amplitude
#    $top = $F / $m;
#    $bot = (($omega0**2.0 - $omega**2.0)**2.0) + ((($b*$b)/($m*$m))*($omega*$omega));
#    $bot = sqrt($bot);
#  
#    $amp = $top / $bot;
#  
#    printf OUTFILE " omega %12.5f  top %12.5f  bot %12.5f  amp %12.5f \n", 
#              $omega, $top, $bot, $amp;
#   
#  }
#  
#  

# 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;

