#!/usr/bin/perl
#
# Compute data showing reflection of waves at a boundary
#    and use the data to make an animated GIF
#
# This one shows a long wave with many cycles approaching the wall.
#
# MWR 4/6/2019

use POSIX;

$mypi = 3.14159;

$debug = 1;

# put datafiles in this subdir
$datadir = "./datadir";

# wavelength of the wave (m)
$lambda = 2.0;
# length of the pulse (m)
$L = $lambda/2.0;

# incoming wave
#    amplitude (m)
$A1 = 0.8;
#    wave number (rad/m)
$k1 = 2.0*$mypi/$lambda;
#    speed of wave (m/s)
$velocity = 1.0;
#    angular frequency (rad/s)
$omega_1 = $k1*$velocity;
# virtual wave wave
#    amplitude (m)
$A2 = 0.8;
#    angular frequency (rad/s)
$omega_2 = $omega_1;
#    wave number (rad/m)
$k2 = $k1;

# plot points at these times
$start_t = 0.0;
$end_t = 8.0;
$dt = 0.10;

# plot points at these locations
$start_x = -10.0;
$end_x = 10.0;
$dx = 0.05;


# create an array of pulse y values
$pulse_num_wiggle = 5.0;
$n_pulse = 200;
$pulse_start = 0.0;
$pulse_end = $pulse_num_wiggle*$L;
$dx = ($pulse_end - $pulse_start) / $n_pulse;
for ($i = 0; $i < $n_pulse; $i++) {
  $x = $pulse_start + $i*$dx;
  $y = $A1*sin($k1*$x);
  $pulse_array_x[$i] = $x;
  $pulse_array_y[$i] = $y;

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


# where does the positive pulse start at t=0?  (m)
$pos_start = 0.0 - (1.5*$pulse_num_wiggle*$L);
# where does the negative pulse start at t=0?  (m)
$neg_start = (0.0 - $pos_start) - $L;


# counter of frames
$frame_count = 0;

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

  # create datafile
  $datafile[$frame_count] = sprintf "%s/data_%05d.dat", 
         $datadir, $frame_count;
  open(DATAFILE, ">$datafile[$frame_count]") || 
              die("can't open file $datafile[$frame_count]");

  # where does the positive pulse begin and end?
  $pos_pulse_start = $pos_start + $t*$velocity;
  $pos_pulse_end = $pos_pulse_start + $pulse_num_wiggle*$L;

  # where does the negative pulse begin and end?
  $neg_pulse_start = $neg_start - $t*$velocity;
  $neg_pulse_end = $neg_pulse_start + $L;

  # plot points at these locations
  for ($x = $start_x; $x <= $end_x; $x += $dx) {

          if (($x >= $pos_pulse_start) && ($x <= $pos_pulse_end)) {
            $index1 = floor(0.5 + ($x - $pos_pulse_start)/$dx);
            $y1 = $pulse_array_y[$index1];
          }
          else {
            $y1 = 0;
          }

          if (($x >= $neg_pulse_start) && ($x <= $neg_pulse_end)) {
            $index2 = floor(0.5 + ($x - $neg_pulse_start)/$dx);
            $y2 = $pulse_array_y[$index2];
          }
          else {
            $y2 = 0;
          }


#          if (($y1 != 0) && ($y2 != 0)) {
#            $sum = $y1 + $y2;
#          }
#          else {
#            $sum = 0.0;
#          }
           if ($x <= 0.0) {
             $sum = $y1 + $y2;
           }
           else {
             $sum = $ymin - 100.0;
           }


#			 $sum = $y1 + $y2;

			 printf DATAFILE "t %8.3f  x %8.3f  y1 %8.3f  y2 %8.3f  sum %8.3f \n", 
					  $t, $x, $y1, $y2, $sum;

  }
  close(DATAFILE);

  $frame_count++;


}


# now plot the data

$xmin = -4;
$xmax = 4;
$ymin = -2;
$ymax = 8;

  my($output_file, $term_options);

  $output_file = "pos_pulse_long.gif";
  $term_options = "gif animate delay 50 loop 1 ";

  $cmdfile = "gnuplot.in";

  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$output_file' \n";
  printf CMDFILE "set term $term_options \n";

  # revert to some oldish default colors
  printf CMDFILE "set style line 1 lt rgb 'red' lw 3  \n";
  printf CMDFILE "set style line 2 lt rgb 'green' lw 3  \n";
  printf CMDFILE "set style line 3 lt rgb 'blue' lw 3  \n";
  printf CMDFILE "set style line 4 lt rgb 'cyan' lw 3  \n";
  printf CMDFILE "set style line 5 lt rgb 'violet' lw 3  \n";



  printf CMDFILE "set grid \n";
  printf CMDFILE "set key top right \n";
  printf CMDFILE "set xrange [$xmin:$xmax]  \n";
  printf CMDFILE "set yrange [$ymin:$ymax]  \n";
  printf CMDFILE "set xlabel 'X position (m)  '  \n";
  printf CMDFILE "set ylabel 'Y position (m) ' \n";
#  printf CMDFILE "set title 'Standing wave ' \n";



  $this_frame_count = 26;

  for ($i = 0; $i < $this_frame_count; $i++) {

    # we need to prune the datafile so that it contains only x <= 0
    $pruned_file[$i] = sprintf "$datadir/pruned_%05d.dat", $i;
    open(DATAFILE, "$datafile[$i]") || die("can't open datafile $datafile[$i]");
    open(PRUNED_FILE, ">$pruned_file[$i]") || die("can't create $pruned_file[$i]
");
    while (<DATAFILE>) {
      $line = " " . $_;
      @words = split(/\s+/, $line);
      $x = $words[4];
      if ($x <= 0) {
        printf PRUNED_FILE "%s", $line;
      }
    }
    close(PRUNED_FILE);
    close(DATAFILE);





    printf CMDFILE "    plot '$pruned_file[$i]' using 4:6 with points ls 1 ps 0.2 pt 6 t '' ";
    printf CMDFILE " \n";
  }

  printf CMDFILE "$cmd \n";


  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;

  $retval = `gnuplot < $cmdfile`;








exit 0;