#!/usr/bin/perl
# 
# Compute positions and velocities of a pair of blocks
#    which are connected by springs of force constant
#
#             k   block   k  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 left-hand block, block1 (m)
$init_x1 = 2.0;
# initial displacement of right-hand block, block2 (m)
$init_x2 = -1.0;
# initial velocity of left-hand block, block1 (m)
$init_v1 = -1.0;
# initial velocity of right-hand block, block2 (m)
$init_v2 = 1.0;
# angular frequency of mode 1 (rad/s)
$omega_1 = sqrt( ($k/$m) );
# angular frequency of mode 2 (rad/s)
$omega_2 = sqrt( (3.0*$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", $init_x1, $init_x2;
printf "# initial velocity v1 %lf  v2 %lf (m) \n", $init_v1, $init_v2;
printf "#  \n";

# after solving for the normal modes, we find 
#   amplitude of normal mode 1 (m)
$A1 = 1.0;
#   phase of normal mode 1 (rad)
$phi_1 = 0.0;
#   amplitude of normal mode 2 (m)
$A2 = -3.22;
#   phase of normal mode 2 (rad)
$phi_2 = 0.3674;



# parameters controlling the time slices in simulation
$num_period = 10.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) {
  
  $x1 = (0.5)*( ($A1*cos($omega_1*$t + $phi_1)) - 
                ($A2*cos($omega_2*$t + $phi_2)) );
  $x2 = (0.5)*( ($A1*cos($omega_1*$t + $phi_1)) + 
                ($A2*cos($omega_2*$t + $phi_2)) );

  $v1 = (0.5)*( (0.0-$omega_1)*$A1*sin($omega_1*$t + $phi_1) -
                (0.0-$omega_2)*$A2*sin($omega_2*$t + $phi_2) );
  $v2 = (0.5)*( (0.0-$omega_1)*$A1*sin($omega_1*$t + $phi_1) +
                (0.0-$omega_2)*$A2*sin($omega_2*$t + $phi_2) );

  $s1 = $A1*cos($omega_1*$t + $phi_1);
  $s2 = $A2*cos($omega_2*$t + $phi_2);

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

}




exit 0;


