#!usr/bin/perl
#
# MWR 4/11/2007

use POSIX;


##########################################################################
# first, the ideal case
#

$outfile = "./ideal_data.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");
$k = 30;
$m = 20.0;
$amp = 1.0;
$phi = 0.0;

$omega = sqrt($k/$m);

$t_start = 0;
$t_end = 100;
$dt = 0.1;
for ($t = $t_start; $t <= $t_end; $t += $dt) {

  $x = $amp * cos($omega*$t + $phi);
  printf OUTFILE ("  t %12.5e  x %12.5e \n", $t, $x);
}
close(OUTFILE);


######################################################################
# Now the damped case

$outfile = "./damped_a_data.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");
$k = 30;
$m = 20.0;
$b = 0.3;
$amp = 1.0;
$phi = 0.0;

$omega = sqrt($k/$m - (($b*$b)/(4*$m*$m)));
$tau = 2*$m/$b;

$t_start = 0;
$t_end = 100;
$dt = 0.1;
for ($t = $t_start; $t <= $t_end; $t += $dt) {

  $x = $amp * cos($omega*$t + $phi);
  $x_exp = exp(-$t/$tau);
  $x = $x * $x_exp;
  printf OUTFILE ("  t %12.5e  x %12.5e \n", $t, $x);
}
close(OUTFILE);



######################################################################
# Now a heavilty damped case

$outfile = "./damped_b_data.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");
$k = 30;
$m = 20.0;
$b = 0.9;
$amp = 1.0;
$phi = 0.0;

$omega = sqrt($k/$m - (($b*$b)/(4*$m*$m)));
$tau = 2*$m/$b;

$t_start = 0;
$t_end = 100;
$dt = 0.1;
for ($t = $t_start; $t <= $t_end; $t += $dt) {

  $x = $amp * cos($omega*$t + $phi);
  $x_exp = exp(-$t/$tau);
  $x = $x * $x_exp;
  printf OUTFILE ("  t %12.5e  x %12.5e \n", $t, $x);
}
close(OUTFILE);




######################################################################
# Now a more heavily damped case

$outfile = "./damped_c_data.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");
$k = 30;
$m = 20.0;
$b = 2.9;
$amp = 1.0;
$phi = 0.0;

$omega = sqrt($k/$m - (($b*$b)/(4*$m*$m)));
$tau = 2*$m/$b;

$t_start = 0;
$t_end = 100;
$dt = 0.1;
for ($t = $t_start; $t <= $t_end; $t += $dt) {

  $x = $amp * cos($omega*$t + $phi);
  $x_exp = exp(-$t/$tau);
  $x = $x * $x_exp;
  printf OUTFILE ("  t %12.5e  x %12.5e \n", $t, $x);
}
close(OUTFILE);

######################################################################
# Now a really heavily damped case

$outfile = "./damped_d_data.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");
$k = 30;
$m = 20.0;
$b = 15;
$amp = 1.0;
$phi = 0.0;

$omega = sqrt($k/$m - (($b*$b)/(4*$m*$m)));
$tau = 2*$m/$b;

$t_start = 0;
$t_end = 100;
$dt = 0.1;
for ($t = $t_start; $t <= $t_end; $t += $dt) {

  $x = $amp * cos($omega*$t + $phi);
  $x_exp = exp(-$t/$tau);
  $x = $x * $x_exp;
  printf OUTFILE ("  t %12.5e  x %12.5e \n", $t, $x);
}
close(OUTFILE);



######################################################################
# Now an OVER-damped oscillator
#   note the slightly different form of the solution

$outfile = "./damped_e_data.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");
$k = 30;
$m = 20.0;
$b = 70;
$amp = 1.0;
$phi = 0.0;

# normal omega would be negative ...
#  $omega = sqrt($k/$m - (($b*$b)/(4*$m*$m)));
# so we use this form instead
$omega = sqrt((($b*$b)/(4*$m*$m)) - ($k/$m));
$tau = 2*$m/$b;

$t_start = 0;
$t_end = 100;
$dt = 0.01;
for ($t = $t_start; $t <= $t_end; $t += $dt) {

  # $x = $amp * cos($omega*$t + $phi);
  # $x_exp = exp(-$t/$tau);

  $x = $amp * cosh($omega*$t);
  $x_exp = exp(-$t/$tau);
  $x = $x * $x_exp;
  printf OUTFILE ("  t %12.5e  x %12.5e \n", $t, $x);
}
close(OUTFILE);




######################################################################
# Now critically-damped oscillator
#   note the slightly different form of the solution

$outfile = "./damped_f_data.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");
$k = 30;
$m = 20.0;
$b = sqrt(4.0*$k*$m);
$amp = 1.0;
$phi = 0.0;

# normal omega would be negative ...
#  $omega = sqrt($k/$m - (($b*$b)/(4*$m*$m)));
# so we use this form instead
$omega = sqrt((($b*$b)/(4*$m*$m)) - ($k/$m));
$tau = 2*$m/$b;
$v0 = 0.0;

$t_start = 0;
$t_end = 100;
$dt = 0.01;
for ($t = $t_start; $t <= $t_end; $t += $dt) {

  # $x = $amp * cos($omega*$t + $phi);
  # $x_exp = exp(-$t/$tau);

  $x = $amp + $amp*(1.0/$tau)*$t + $v0*$t;
  $x_exp = exp(-$t/$tau);
  $x = $x * $x_exp;
  printf OUTFILE ("  t %12.5e  x %12.5e \n", $t, $x);
}
close(OUTFILE);






exit 0;