#!/usr/bin/perl
# 
# Compute positions and velocities of a pair of blocks
#    which are connected by springs of force constant
#
#             k   block   kc  block  k
#
# MWR 2/24/2019

$debug = 1;

$mypi = 3.14159;

# force constant of side springs (N/m)
$k = 1.0;
# force constant of middle "coupling" spring (N/m)
$kc = 1.00;
# mass of block (kg)
$m = 1.0;
# initial displacement of right-hand block, block2 (m)
$L = 0.10;
# angular frequency of mode 1 (rad/s)
$omega_1 = sqrt( ($k/$m) + (2.0*$kc/$m) );
# angular frequency of mode 2 (rad/s)
$omega_2 = sqrt( ($k/$m) );
# Period of mode 1 (s)
$period_1 = (2.0*$mypi/$omega_1);
# Period of mode 2 (s)
$period_2 = (2.0*$mypi/$omega_2);

printf "# mass of blocks m %lf (kg) \n", $m;
printf "# spring const k %lf (N/m)  kc %lf (N/m) \n", $k, $kc;
printf "# omega_1 %lf (rad/s)  omega_2 %lf (rad/s) \n", $omega_1, $omega_2;
printf "# period_1 %lf (s)  period_2 %lf (s) \n", $period_1, $period_2;
printf "# initial position x1 %lf  x2 %lf (m) \n", 0, $L;
printf "# initial velocity v1 %lf  v2 %lf (m) \n", 0, 0;
printf "#  \n";


# parameters controlling the time slices in simulation
$num_period = 100.0;
$num_per_period = 300.0;
$start_time = 0;
$end_time = $num_period*$period_2;
$d_time = $period_2 / $num_per_period;

for ($t = $start_time; $t <= $end_time; $t += $d_time) {
  
  $x1 = ($L/2.0)*(cos($omega_2*$t) - cos($omega_1*$t));
  $x2 = ($L/2.0)*(cos($omega_2*$t) + cos($omega_1*$t));

  $v1 = (($L/2.0)*(0.0 - $omega_2))*(sin($omega_2*$t) - sin($omega_1*$t));
  $v2 = (($L/2.0)*(0.0 - $omega_2))*(sin($omega_2*$t) + sin($omega_1*$t));

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

}




exit 0;


