#!/usr/bin/perl
# 
# generate some data to illustrate behavior of 
#   products of sines and sines
#
# MWR 2/3/2019

$debug = 1;

# avoid dividing by zero
$TINY = 1.0E-10;
$BIG = 1.0/$TINY;


$mypi = 3.14159;

# period of fundamental 
$period = 10.0;
# fundamental frequency 
$omega_1 = (2.0*$mypi)/$period;

# range of time to use in calculations
$n_per_period = 500;
$n_period = 1.0;
$start_t = 0;
$end_t = $n_period*$period;
$d_t = ($period/$n_per_period);

$outfile = "./make_sin_sin.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");

printf OUTFILE "# data generated to illustrate sines and sines \n";
printf OUTFILE "# period P = $period  (s) \n";
printf OUTFILE "# fundamental frequency omega_1 = $omega_1  (rad/s) \n";
printf OUTFILE "# \n";

for ($t = $start_t; $t <= $end_t; $t += $d_t) {

  # sin of omega_1
  $sin_1 = sin( $omega_1*$t );

  # sin of omega_2
  $sin_2 = sin( (2.0*$omega_1)*$t );

  # sin of omega_3
  $sin_3 = sin( (3.0*$omega_1)*$t );


  printf OUTFILE " t %12.5f  sin1 %12.5f  sin2 %12.5f  sin3  %12.5f   \n",
            $t, $sin_1, $sin_2, $sin_3;
 
}


#  
#  # now print a second stanza, which shows the behavior for
#  #   a larger damping 
#  $b = 0.8;
#  
#  
#  printf OUTFILE "\n\n\n";
#  
#  printf OUTFILE "# data for an intermediate damping \n";
#  printf OUTFILE "# force F = $F (N) \n";
#  printf OUTFILE "# mass m = $m (kg) \n";
#  printf OUTFILE "# spring constant  k = $k (N/m) \n";
#  printf OUTFILE "# air resistance coeff b = $b  (N/m/s) \n";
#  printf OUTFILE "# phase offset  phi = $phi (rad) \n";
#  printf OUTFILE "# undamped frequency omega0 = $omega0  (rad/s) \n";
#  printf OUTFILE "# time constant tau = $tau  (s) \n";
#  printf OUTFILE "# Q = $Q  \n";
#  printf OUTFILE "# \n";
#  
#  for ($omega = $start_omega; $omega <= $end_omega; $omega += $d_omega) {
#  
#    # amplitude
#    $top = $F / $m;
#    $bot = (($omega0**2.0 - $omega**2.0)**2.0) + ((($b*$b)/($m*$m))*($omega*$omega));
#    $bot = sqrt($bot);
#  
#    $amp = $top / $bot;
#  
#    printf OUTFILE " omega %12.5f  top %12.5f  bot %12.5f  amp %12.5f \n", 
#              $omega, $top, $bot, $amp;
#   
#  }
#  
#  
#  
#  # now print a third stanza, which shows the behavior for
#  #   a very large damping 
#  $b = 2.0;
#  
#  
#  printf OUTFILE "\n\n\n";
#  
#  printf OUTFILE "# data for a large damping \n";
#  printf OUTFILE "# force F = $F (N) \n";
#  printf OUTFILE "# mass m = $m (kg) \n";
#  printf OUTFILE "# spring constant  k = $k (N/m) \n";
#  printf OUTFILE "# air resistance coeff b = $b  (N/m/s) \n";
#  printf OUTFILE "# phase offset  phi = $phi (rad) \n";
#  printf OUTFILE "# undamped frequency omega0 = $omega0  (rad/s) \n";
#  printf OUTFILE "# time constant tau = $tau  (s) \n";
#  printf OUTFILE "# Q = $Q  \n";
#  printf OUTFILE "# \n";
#  
#  for ($omega = $start_omega; $omega <= $end_omega; $omega += $d_omega) {
#  
#    # amplitude
#    $top = $F / $m;
#    $bot = (($omega0**2.0 - $omega**2.0)**2.0) + ((($b*$b)/($m*$m))*($omega*$omega));
#    $bot = sqrt($bot);
#  
#    $amp = $top / $bot;
#  
#    printf OUTFILE " omega %12.5f  top %12.5f  bot %12.5f  amp %12.5f \n", 
#              $omega, $top, $bot, $amp;
#   
#  }
#  
#  

# now print a second stanza which simply contains 
#   the frequency of the undamped system
#printf OUTFILE "\n\n\n";
#printf OUTFILE "# this is the frequency of the undamped system \n";
#printf OUTFILE " omega0 %12.5f  small  -1000 \n", $omega0;
#printf OUTFILE " omega0 %12.5f  zero       0 \n", $omega0;
#printf OUTFILE " omega0 %12.5f  big     1000 \n", $omega0;
#printf OUTFILE "\n\n\n";



close(OUTFILE);


exit 0;

