#!/usr/bin/perl
#
# Compute some data we can use to illustrate the idea
#     of wave pulses.
#
# MWR 2/8/2019

$debug = 0;

$mypi = 3.14159;

# fundamental frequency (rad/s)
$omega0 = 1.0;
# amplitude of fundamental frequency (m)
$amp0 = 1.0;
# starting time for the "wide" view (s)
$start_time_wide = 0.0;
$start_time_wide = -10.0;
# ending time for the "wide" view (s)
$end_time_wide = 100.0;
# sampling of the "wide" view (s)
$d_time_wide = 0.0008;


# compute values of a sinc function,
#    which we will use to set amplitudes for off-center frequencies
#    save $num_sinc results in "sinc_amp_array[]",
#    and corresponding frequencies in "sinc_freq_array[]"
$sinc_center_freq = 1000.0;
$freq_width = 40;
$sinc_width = 5000.0;
$num_sinc = 19.0;

$sinc_start = ($sinc_center_freq - ($sinc_width));
$sinc_end   = ($sinc_center_freq + ($sinc_width));
$d_sinc     = ($sinc_end - $sinc_start) / ($num_sinc - 1.0);
$freq_start = ($sinc_center_freq - ($freq_width));
$freq_end   = ($sinc_center_freq + ($freq_width));
$d_freq     = ($freq_end - $freq_start) / ($num_sinc - 1.0);
if ($debug > 0) {
  printf " sinc_start %.3f  sinc_end %.3f  d_sinc %.3f \n",
      $sinc_start, $sinc_end, $d_sinc;
}

printf "# first, the frequencies and amplitudes we'll add \n";
for ($i = 0; $i < $num_sinc; $i++) {

  $x = (($sinc_start + $i*$d_sinc) / $sinc_center_freq) - 1.0;
  $freq = $freq_start + $i*$d_freq;

  if ($debug > 0) {
    printf " i %5d   freq %8.3f  x %12.5f  ", $i, $freq, $x;
  }

  # a normalization factor
  $norm = 3.0/($num_sinc);
  if ($x == 0) {
    $sinc_amp = $norm;
  }
  else {
    $sinc_amp = $norm*(sin($x) / ($x));
  }
  if ($debug > 0) {
    printf " %12.5f \n", $sinc_amp;
  }

  $sinc_amp_array[$i] = $sinc_amp;
  $sinc_freq_array[$i] = $freq;


  printf " i %5d   freq %8.3f  x %12.5f  amp %12.5f \n", $i, $freq, $x, $sinc_amp;
}
printf "\n\n\n";




  # compute data for a sum of frequencies
  printf "# now, values of the combined waves \n";
  for ($t = $start_time_wide; $t <= $end_time_wide; $t += $d_time_wide) {
 
     $sum = 0;
     for ($i = 0; $i < $num_sinc; $i++) {

       $this_freq = $sinc_freq_array[$i];
       $this_amp  = $sinc_amp_array[$i];

       if ($debug > 0) {
         printf "#  t %12.5f  freq %8.3f  amp %8.3f ", 
               $t, $this_freq, $this_amp;
       }

       $val = $this_amp*$amp0*cos($this_freq*$t);
       $sum += $val;
      
       if ($debug > 0) {
         printf "   val %12.5f sum %12.5f \n", $val, $sum;
       }


     }


     printf " %12.5f %12.5f \n", $t, $sum;

  }








exit 0;
