#!/usr/bin/perl
#
# Compute the motions of 4 identical blocks connected by 5 identical springs
#    Mode 2 only
#
# MWR 2/28/2019

$debug = 1;
$mypi = 3.14159;


# mass of each block (kg)
$m = 1.0;
# force constant of each spring (N/m)
$k = 1.0;
# frequencies of the normal modes
$omega_1 = sqrt( 0.38197  * ($k/$m) );
$omega_2 = sqrt( 3.61803  * ($k/$m) );
$omega_3 = sqrt( 2.61803  * ($k/$m) );
$omega_4 = sqrt( 1.38197  * ($k/$m) );
# periods of the normal modes
$period_1 = (2.0*$mypi)/($omega_1);
$period_2 = (2.0*$mypi)/($omega_2);
$period_3 = (2.0*$mypi)/($omega_3);
$period_4 = (2.0*$mypi)/($omega_4);
# initial conditions --
#     block 1 is moved +L, all other blocks at zero,
#     all blocks motionless
$L = 1.0;

printf "# mass m %9.3f \n", $m;
printf "# force constant k %9.3f \n", $k;
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 "# omega_3 %9.3f    period_3 %9.3f \n", $omega_3, $period_3;
printf "# omega_4 %9.3f    period_4 %9.3f \n", $omega_4, $period_4;
printf "# initial x1 L %9.3f \n", $L;
printf "#  \n";


# these two numbers come up again and again in calculations
$alpha = 0.37175;
$beta = 0.60150;

$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
  $x1 = $L*(    0*$alpha*$alpha*cos($omega_1*$t)   +   1*$beta*$beta*cos($omega_2*$t)   
             +  0*$beta*$beta*cos($omega_3*$t)     +   0*$alpha*$alpha*cos($omega_4*$t)  );
  $x2 = $L*( -  0*$alpha*$alpha*cos($omega_1*$t)   +   1*$beta*$beta*cos($omega_2*$t)   
             -  0*$beta*$beta*cos($omega_3*$t)     +   0*$alpha*$alpha*cos($omega_4*$t)  );
  $x3 = $L*(    0*$alpha*$beta*cos($omega_1*$t)    -   1*$alpha*$beta*cos($omega_2*$t)   
             -  0*$alpha*$beta*cos($omega_3*$t)    +   0*$alpha*$beta*cos($omega_4*$t)  );
  $x4 = $L*(    0*$alpha*$beta*cos($omega_1*$t)    +   1*$alpha*$beta*cos($omega_2*$t)   
             -  0*$alpha*$beta*cos($omega_3*$t)    -   0*$alpha*$beta*cos($omega_4*$t)  );


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


}


exit 0;
