#!/usr/bin/perl
# 
# generate some data for illustrating damped harmonic motion
#
# MWR 1/24/2019

$debug = 1;

$mypi = 3.14159;

# put output into this file
$outfile = "./make_damped.out";

# amplitude (m)
$A = 4.0;
# mass (kg)
$m = 1.0;
# force constant (N/m)
$k = 10.0;
# resistive force coeff (N/m/s)
$b = 0.0;
# phase offset (rad)
$phi = 0;

open(OUTFILE, ">$outfile") || die("can't open file $outfile");


# generate undamped motion

$omega = sqrt($k/$m);
$period = (2*$mypi)/$omega;
$start_t = 0;
$end_t = 30*$period;
$dt = $period*0.05;


printf OUTFILE "# here comes undamped motion \n";
printf OUTFILE "# period = %12.5f  dt = %12.5f  \n", $period, $dt;
for ($t = $start_t; $t <= $end_t; $t += $dt) {
  $x = $A * cos($omega*$t + $phi);
  printf OUTFILE " %8.2f  %8.2f \n", $t, $x;
}
printf OUTFILE "\n\n\n";



# generate lightly damped motion

$b = 0.02*sqrt(4.0*$k*$m);
$omega = sqrt( ($k/$m) - ($b*$b)/(4.0*$m*$m) );
$tau = (2.0*$m/$b);
$period = (2*$mypi)/$omega;
$start_t = 0;
$end_t = 30*$period;
$dt = $period*0.05;


printf OUTFILE "# here comes lightly damped motion \n";
printf OUTFILE "# b = %12.5f  tau = %12.5f \n", $b, $tau;
for ($t = $start_t; $t <= $end_t; $t += $dt) {
  $x = $A * cos($omega*$t + $phi) * exp( -$t/$tau);
  printf OUTFILE " %8.2f  %8.2f \n", $t, $x;
}
printf OUTFILE "\n\n\n";



# generate critically damped motion

$b = sqrt(4.0*$k*$m);
#$omega = sqrt( ($k/$m) - ($b*$b)/(4.0*$m*$m) );
$omega = 0.0;
$tau = (2.0*$m/$b);

# $period = (2*$mypi)/$omega;
$start_t = 0;
#$end_t = 30*$period;
$end_t = 60.0;
$dt = 0.09935;


printf OUTFILE "# here comes lightly damped motion \n";
printf OUTFILE "# b = %12.5f  tau = %12.5f \n", $b, $tau;
for ($t = $start_t; $t <= $end_t; $t += $dt) {
  $x = $A * cos($omega*$t + $phi) * exp( -$t/$tau);
  printf OUTFILE " %8.2f  %8.2f \n", $t, $x;
}
printf OUTFILE "\n\n\n";




# generate over-damped motion

$b = 5.0*sqrt(4.0*$k*$m);
$gamma = sqrt( ($b*$b)/(4.0*$m*$m) - ($k/$m) );
$tau = (2.0*$m/$b);

# $period = (2*$mypi)/$omega;
$start_t = 0;
#$end_t = 30*$period;
$end_t = 60.0;
$dt = 0.09935;


printf OUTFILE "# here comes over damped motion \n";
printf OUTFILE "# b = %12.5f  tau = %12.5f  gamma = %12.5f  \n", $b, $tau, $gamma;
for ($t = $start_t; $t <= $end_t; $t += $dt) {
  $x = $A * exp($gamma*$t + $phi) * exp( -$t/$tau);
  printf OUTFILE " %8.2f  %8.2f \n", $t, $x;
}
printf OUTFILE "\n\n\n";




close(OUTFILE);

exit 0;



