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

$debug = 0;
$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( 1.38197  * ($k/$m) );
$omega_3 = sqrt( 2.61803  * ($k/$m) );
$omega_4 = sqrt( 3.61803  * ($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 numbers come up again and again in calculations
$alpha = 0.138196;
$beta = 0.361804;
$gamma = 0.223608; 

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

# this allows us to pick the mode of interest
$s1 = 0.0;
$s2 = 1.0;
$s3 = 0.0;
$s4 = 0.0;

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

  # position of each block
  $x1 = $L*(   $s1* $alpha*cos($omega_1*$t)   - $s2*  $beta*cos($omega_2*$t)  
             + $s3*  $beta*cos($omega_3*$t)   - $s4* $alpha*cos($omega_4*$t)  );
  $x2 = $L*(   $s1* $gamma*cos($omega_1*$t)   - $s2* $gamma*cos($omega_2*$t)   
             - $s3* $gamma*cos($omega_3*$t)   + $s4* $gamma*cos($omega_4*$t)  );
  $x3 = $L*(   $s1* $gamma*cos($omega_1*$t)   + $s2* $gamma*cos($omega_2*$t)   
             - $s3* $gamma*cos($omega_3*$t)   - $s4* $gamma*cos($omega_4*$t)  );
  $x4 = $L*(   $s1* $alpha*cos($omega_1*$t)   + $s2*  $beta*cos($omega_2*$t)   
             + $s3*  $beta*cos($omega_3*$t)   + $s4* $alpha*cos($omega_4*$t)  );



  if (1 == 1) {
    printf "t %9.3f  x1 %9.3f  x2 %9.3f  x3 %9.3f   x4 %9.3f \n",
         $t, $x1, $x2, $x3, $x4;
  }
  else {
    printf "t %9.3f  1 %9.3f \n", $t, $x1;
    printf "t %9.3f  2 %9.3f \n", $t, $x2;
    printf "t %9.3f  3 %9.3f \n", $t, $x3;
    printf "t %9.3f  4 %9.3f \n", $t, $x4;
    printf "\n\n";
  }


}


exit 0;
