#!/usr/bin/perl
#
# Compute data showing reflection of waves at a boundary
#    and use the data to make an animated GIF
#
# This version shows a positive incoming and NEGATIVE reflected wave
#    appropriate for fixed end of string at x=0
#
# 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 = 4.0;
$dt = 0.10;

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


# create an array of pulse y values
$n_pulse = 200;
$pulse_start = 0.0;
$pulse_end = $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;
}


# where does the positive pulse start at t=0?  (m)
$pos_start = 0.0 - (2.0*$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 + $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 = 0.0 - $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
#    This shows the incident, reflected, AND summed waves
#         refl_b.gif


$xmin = -4;
$xmax = 4;
$ymin = -1;
$ymax = 6;

  my($output_file, $term_options);

  $output_file = "refl_b.gif";
  $term_options = "gif animate delay 50  ";

  $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";

  for ($i = 0; $i < $frame_count; $i++) {
    printf CMDFILE "    plot '$datafile[$i]' using 4:6 with lines ls 1  t '' ";
    printf CMDFILE " , ";
    printf CMDFILE "         '$datafile[$i]' using 4:(\$8+2) with lines ls 2  t '' ";
    printf CMDFILE " , ";
    printf CMDFILE "         '$datafile[$i]' using 4:(\$10+4) with points ls 3 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`;



# now plot the data
#    This shows the incident, reflected, but not summed wave
#         refl_b2.gif


$xmin = -4;
$xmax = 4;
$ymin = -1;
$ymax = 6;

  my($output_file, $term_options);

  $output_file = "refl_b2.gif";
  $term_options = "gif animate delay 50  ";

  $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";

  for ($i = 0; $i < $frame_count; $i++) {
    printf CMDFILE "    plot '$datafile[$i]' using 4:6 with lines ls 1  t '' ";
    printf CMDFILE " , ";
    printf CMDFILE "         '$datafile[$i]' using 4:(\$8+2) with lines ls 2  t '' ";
    printf CMDFILE " \n";
  }

  printf CMDFILE "$cmd \n";


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

  $retval = `gnuplot < $cmdfile`;










exit 0;