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

$debug = 0;

$mypi = 3.14159;

# force constant of side springs (N/m)
$k = 6.0;
# mass of left-hand block (kg)
$m = 3.0;
# mass of right-hand block (kg)
$big_m = 3.0*$m;
# initial displacement of right-hand block, block2 (m)
$L = 0.55;
# angular frequency of mode 1 (rad/s)
$omega_1 = sqrt( ($k/$m) * ((8.0 - sqrt(28.0))/3.0) );
# angular frequency of mode 2 (rad/s)
$omega_2 = sqrt( ($k/$m) * ((8.0 + sqrt(28.0))/3.0) );
# Period of mode 1 (s)
$period_1 = (2.0*$mypi/$omega_1);
# Period of mode 2 (s)
$period_2 = (2.0*$mypi/$omega_2);
# useful quantities
$alpha = ((8.0 + sqrt(28.0))/3.0) - 2.0;
$beta = 2.0 - ((8.0 - sqrt(28.0))/3.0);

printf "# mass of blocks m %lf  m2 %lf (kg) \n", $m, $big_m;
printf "# spring const k %lf (N/m)  (N/m) \n", $k;
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 "# alpha %lf  beta %lf \n", $alpha, $beta;
printf "#  \n";



# parameters controlling the time slices in simulation
$num_period = 100.0;
$num_per_period = 50.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) {
  
  $z1 = ($L/(1.0 + ($alpha/$beta)));
  $z2 = (0.0 - $alpha)*cos($omega_2*$t);
  $z3 = ($alpha)*cos($omega_1*$t);

  $x1 = $z1*($z2 + $z3);

  if ($debug > 0) {
    printf "#Q t %9.3f  x1 %9.3f   z1 %9.3f z2 %9.3f z3 %9.3f    \n",
             $t, $x1, $z1, $z2, $z3;
  }

  $z1 = ($L/($alpha + $beta));
  $z2 = $beta*cos($omega_1*$t);
  $z3 = ($alpha)*cos($omega_2*$t);

  $x2 = $z1*($z2 + $z3);

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

}




exit 0;


