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

$debug = 1;

$mypi = 3.14159;

# force constant of first spring (N/m)
$k1 = 30.0;
# force constant middle spring (N/m)
$k2 = 2.0*$k1;
# force constant right-most spring (N/m)
$k3 = 3.0*$k1;
# mass of each block (kg)
$m = 2.0;
# initial displacement of right-hand block, block2 (m)
$L = 1.00;
# angular frequency of mode 1 (rad/s)
$omega_1 = sqrt(  (4.0 - sqrt(5.0))*($k1/$m) );
# angular frequency of mode 2 (rad/s)
$omega_2 = sqrt(  (4.0 + sqrt(5.0))*($k1/$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 k1 %lf k2 %lf k3 %lf  (N/m) \n", $k1, $k2, $k3;
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";


# these are the amplitudes of "x1" component of the normal modes
$s1_1 = 1.618;
$s1_2 = -0.618;
# these are the amplitudes of "x2" component of the normal modes
$s2_1 = 1.0;
$s2_2 = 1.0;
# and these are the amplitudes of x1, in terms of the normal modes
$x1_1 = 1.0 / (($s1_1 - $s1_2));
$x1_2 = -$x1_1;
# and these are the amplitudes of x2, in terms of the normal modes
$x2_1 = 1.0 / (($s1_1/abs($s1_2)) + 1.0);
$x2_2 = $x2_1 * (($s1_1/abs($s1_2)));

printf "# x1_1  %8.3f  x1_2  %8.3f \n", $x1_1, $x1_2;
printf "# x2_1  %8.3f  x2_2  %8.3f \n", $x2_1, $x2_2;

# parameters controlling the time slices in simulation
$num_period = 10.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 = ($x1_1*$L)*cos($omega_1*$t) + ($x1_2*$L)*cos($omega_2*$t);
  $x2 = ($x2_1*$L)*cos($omega_1*$t) + ($x2_2*$L)*cos($omega_2*$t);

  $v1 =     ($x1_1*$L)*(0.0 - $omega_1)*sin($omega_1*$t)
         +  ($x1_2*$L)*(0.0 - $omega_2)*sin($omega_2*$t);
  $v2 =     ($x2_1*$L)*(0.0 - $omega_1)*sin($omega_1*$t)
         +  ($x2_2*$L)*(0.0 - $omega_2)*sin($omega_2*$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;


