#!/usr/bin/perl
#
# Create an animation showing two travelling waves:
#       y1 = A sin(kx - wt)
#       y2 = A sin(kx + wt)
#    and their superposition.  Highlight the area near x=0
#
# MWR 4/9/2019

$mypi = 3.14159;

$datadir = "./stand_data";
$gif_delay = 10;

$num_snapshots = 80;
$wavelength = 5.0;
$k = (2.0*$mypi)/($wavelength);
$xmin = -10.0;
$xmax = 10.0;
$A = 1.0;
$period = 3.0;
$omega = 2.0*$mypi/($period);
$speed = $omega / $k;
$dx = 0.03;
$t_start = 0;
$t_end = 15.0;
$dt = ($t_end - $t_start) / $num_snapshots;

# where does the right-going wave start?
$right_start = 0.0 - 1.5*$wavelength;
# where does the left-going wave start?
$left_start = 0.0 - $right_start;

# when will the right-going wave reach the position x=0?
$meet_time = ($left_start / $speed);

# this y value will place a symbol far outside our graph
$big_y = 1000;

# step one: create one datafile per snapshot
for ($i = 0; $i < $num_snapshots; $i++) {
  $datafile[$i] = sprintf "%s/invert_%02d.dat", $datadir, $i;
  open(DATAFILE, ">$datafile[$i]") || die("can't open file $datafile[$i]");

  $t = $t_start + ($i*$dt);

  # constant amplitude
  $amp = $A;

  # what is the position of the leading edge of the wave going to right?
  $right_x = $right_start + $t*$speed;
  # what is the position of the leading edge of the wave going to left?
  $left_x = $left_start - $t*$speed;

  for ($x = $xmin; $x <= $xmax; $x += $dx) {

    if ($x <= $right_x) {
      $right_y = $amp*sin($k*$x - $omega*$t);
    }
    elsif (($t < $meet_time) && ($x < 0.0)) {
      $right_y = 0.0;
    }
    else {
      $right_y = $big_y;
    }

    if ($x >= $left_x) {
      $left_y = $amp*sin($k*$x + $omega*$t);
    }
    elsif (($t < $meet_time) && ($x > 0.0)) {
      $left_y = 0.0;
    }
    else {
      $left_y = $big_y;
    }

    # this is the sum of the two waves
    if ($x <= 0.0) {
      $sum_y = $right_y + $left_y;
    }
    else {
      $sum_y = $big_y;
    }

    printf DATAFILE "%lf %lf %lf %lf \n", $x, $right_y, $left_y, $sum_y;
  }
  close(DATAFILE);

}


  


# step three: plot each datafile to create an animated GIF image 
#                  which shows
#                      a) wave travelling in positive x dir
#                      b) wave travelling in negative x dir
#                      c) sum of the two waves

  my($output_file, $term_options);

  $sum_gif_delay = 2*$gif_delay;

  $output_file = "invert_a.gif";
  $term_options = "gif animate delay $sum_gif_delay  ";

  $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 'red' lw 1   \n";
  printf CMDFILE "set style line 5 lt rgb 'blue' lw 1  \n";
  printf CMDFILE "set style line 6 lt rgb 'black' lw 2  \n";

  printf CMDFILE "set object 1 rectangle from -0.5,6.0 to 0.5,10.0 fs empty border lt -1    \n";


  printf CMDFILE "set grid \n";
  printf CMDFILE "set key top right \n";
  printf CMDFILE "set xrange [$xmin:$xmax]  \n";
  printf CMDFILE "set yrange [-2:15]  \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 < $num_snapshots; $i++) {
    printf CMDFILE "    plot '$datafile[$i]' using 1:2 with points ls 4  ps 0.2 pt 6 t '' ";
    printf CMDFILE " , ";
    printf CMDFILE "       '$datafile[$i]' using 1:(\$3+4) with points ls 5  ps 0.2 pt 6  t '' ";
    printf CMDFILE " , ";
    printf CMDFILE "       '$datafile[$i]' using 1:(\$2+8) with points ls 4 ps 0.2 pt 6  t '' ";
    printf CMDFILE " , ";
    printf CMDFILE "       '$datafile[$i]' using 1:(\$3+8) with points ls 5 ps 0.2 pt 6  t '' ";
    printf CMDFILE " , ";
    printf CMDFILE "       '$datafile[$i]' using 1:(\$4+12) with points ls 6 ps 0.05 pt 6  t '' ";
    printf CMDFILE " , ";
    printf CMDFILE " 0.0 with lines ls -1 lw 1 t ''  ";
    printf CMDFILE " , ";
    printf CMDFILE " 4.0 with lines ls -1 lw 1 t ''  ";
    printf CMDFILE " , ";
    printf CMDFILE " 8.0 with lines ls -1 lw 1 t ''  ";
    printf CMDFILE " , ";
    printf CMDFILE " 12.0 with lines ls -1 lw 1 t ''  ";
    printf CMDFILE "$cmd \n";
  }

  printf CMDFILE "$cmd \n";


  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;




exit 0;
  


# step four: plot each datafile to create an animated GIF image 
#                  which shows
#                      c) sum of the two waves only, zoomed in a bit
 
  


  my($output_file, $term_options);

  $sum_gif_delay = 2*$gif_delay;

  # this file will hold a pruned version of the original datafiles
  #    that we use for plotting this figure only
  $xmin = -4.0;
  $xmax = 4.0;

  $output_file = "travel_d.gif";
  $term_options = "gif animate delay $sum_gif_delay  ";

  $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 'red' lw 1   \n";
  printf CMDFILE "set style line 5 lt rgb 'blue' lw 1  \n";
  printf CMDFILE "set style line 6 lt rgb 'black' lw 4  \n";


  printf CMDFILE "set grid \n";
  printf CMDFILE "set key top right \n";
  printf CMDFILE "set xrange [$xmin:$xmax]  \n";
  printf CMDFILE "set yrange [0-($A*2.2):0+($A*2.2)]  \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 < $num_snapshots; $i++) {

    # we need to prune the datafile so that it runs only from
    #   the fixed positions at
    #      $stand_wave_xmin to $stand_wave_xmax
    $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[1];
      if (($x >= ($stand_wave_xmin - 0.0)) && ($x <= ($stand_wave_xmax + 0.0))) {
        printf PRUNED_FILE "%s", $line;
      }
    }
    close(PRUNED_FILE);
    close(DATAFILE);



    printf CMDFILE "    plot '$pruned_file[$i]' using 1:3 with lines ls 6 lw 2 t '' ";
    printf CMDFILE "$cmd \n";
  }

  printf CMDFILE "$cmd \n";


  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;




exit 0;