#!/usr/bin/perl
#
# Compute the motions of two identical blocks connected by 
#    three identical springs, in the presence of damping
#
# MWR 3/5/2019

$debug = 1;
$mypi = 3.14159;

# mass of each block (kg)
$m = 1.0;
# force constant of each spring (N/m)
$k = 1.0;
# resistive force coeff (N/m/s/)
$b = 0.1;
# frequencies of the normal modes
$omega_1 = sqrt( ($k/$m) - (($b*$b/(4.0*$m*$m))) );
$omega_2 = sqrt( (3.0*$k/$m) - (($b*$b/(4.0*$m*$m))) );
# periods of the normal modes
$period_1 = (2.0*$mypi)/($omega_1);
$period_2 = (2.0*$mypi)/($omega_2);
# time constant tau (s)
$tau = 2.0*$m/$b;
# initial conditions --
#     block 1 is moved +L, all other blocks at zero,
#     all blocks motionless
$L = 1.0;
$A1 = 0.5*$L;
$A2 = -0.5*$L;

printf "# mass m %9.3f \n", $m;
printf "# force constant k %9.3f \n", $k;
printf "# resistive force coeff b %9.3f \n", $b;
printf "# omega_1 %9.3f    period_1 %9.3f \n", $omega_1, $period_1;
printf "# omega_2 %9.3f    period_2 %9.3f \n", $omega_2, $period_2;
printf "# initial x1 L %9.3f \n", $L;
printf "#  \n";


$num_period = 20;
$start_time = 0; 
$end_time = $num_period * $period_1;
$num_per_period = 100;
$d_t = $period_1 / $num_per_period;

for ($t = $start_time; $t <= $end_time; $t += $d_t) {

  # position of each block
  $s1 = $A1*cos($omega_1*$t)*exp(-$t/$tau);
  $s2 = $A2*cos($omega_2*$t)*exp(-$t/$tau);

  $x1 = 0.5*($s1 - $s2);
  $x2 = 0.5*($s1 + $s2);

  printf "t %9.3f  x1 %9.3f  x2 %9.3f   \n",
         $t, $x1, $x2;


}


exit 0;
