#!/usr/bin/perl
# 
# generate some data to illustrate behavior of 
#   amplitude for different Q values
#
# MWR 1/29/2019

$debug = 1;

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


$mypi = 3.14159;

# initial amplitude (m)
$A = 4.0;
# 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;
# phase offset (rad)
$phi = 0;
# 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;

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


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

printf OUTFILE "# data generated to illustrate Q \n";
printf OUTFILE "# amplitude A = $A (m) \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 ($t = $start_t; $t <= $end_t; $t += $d_t) {

  # amplitude
  $exp_factor = exp(0.0 - ($t/$tau));
  $pos = $A * cos($omega0*$t + $phi) * $exp_factor;

  printf OUTFILE " t %12.5f  exp_factor %12.5f  pos %12.5f \n", 
            $t, $exp_factor, $pos;
 
}

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

