From 71bd494206276ac07cb3e760a1a1897783f9255a Mon Sep 17 00:00:00 2001
From: caes <xaw5719@wmich.edu>
Date: Fri, 12 Aug 2016 05:14:03 -0400
Subject: [PATCH] updating tophat fft

---
 plot.sh               |  7 +++--
 scripts/tophat_fft.pl | 60 ++++++++++++++++++++++++++-----------------
 2 files changed, 40 insertions(+), 27 deletions(-)

diff --git a/plot.sh b/plot.sh
index d4ceb99..9b4b55c 100755
--- a/plot.sh
+++ b/plot.sh
@@ -30,7 +30,6 @@ case $1 in
         done
         echo "$gnuplot_input"|perl -pe 's|␤|\n|g' > ${gnuplot_file}
         gnuplot $gnuplot_file
-
     ;;
 
     "lags"|"lag"|"delay"|"delays")
@@ -52,10 +51,10 @@ case $1 in
     ;;
 
     "tophat")
-        gnuplot_file=tophat_w_fft.gp
+        gnuplot_file="tophat_w_fft.gp"
     ;;
-    
+
     *)
         echo "Did not understand plot type."
-    ;; 
+    ;;
 esac
\ No newline at end of file
diff --git a/scripts/tophat_fft.pl b/scripts/tophat_fft.pl
index de9faa0..a1821ea 100755
--- a/scripts/tophat_fft.pl
+++ b/scripts/tophat_fft.pl
@@ -1,7 +1,7 @@
 #!/usr/bin/env perl
 
 use feature say;
-use utf8; 
+use utf8;
 use PDL;
 use PDL::FFT;
 use PDL::IO::Misc;
@@ -9,38 +9,52 @@ use PDL::IO::Misc;
 my $twoPI = 2*3.1415926539;
 
 # space parameters
-our $Δt=.001;
+our $xres=.001;
 our $xmin = 0;
 our $xmax = 250;
-our $num_bins = ($xmax-$xmin)/$Δt;
+our $num_bins = ($xmax-$xmin)/$xres;
+
+# tophat parameters
+my $μ1=15;
+my $Δ1=5; # full width
+my $μ2=15;
+my $Δ2=5; # full width
+my $μ3=15;
+my $Δ3=5; # full width
 
 say "Creating $num_bins bins";
 
 # Complex Number z(x) = u(x) + iv(x)
-my @th_list=tophat(10,7);
+my @th_list=tophat($μ1,$Δ1);
 my $u = pdl(@th_list);
 my $v = zeroes($u);
 
-# Output time-domain tophat
-open $tophattab, '>', 'analyses/tables/tophat.tab';
-for ($i=0; $i < $num_bins; $i++) {
-    my $x_coord = ($xmax-$xmin)*($i/$num_bins) + $xmin;
-    say $tophattab "$x_coord $th_list[$i]";
+
+foreach (('1','2','3')) {
+
+    # Output time-domain tophat
+    open $tophattab, '>', "analyses/tables/tophat${$_}.tab";
+    for ($i=0; $i < $num_bins; $i++) {
+        my $x_coord = ($xmax-$xmin)*($i/$num_bins) + $xmin;
+        say $tophattab "$x_coord $th_list[$i]";
+    }
+    close $tophattab;
+
+    # Transform z(x) to Z(1/x) = U(1/x) + iV(1/x)
+    my $U = $u;
+    my $V = $v;
+    fft($U,$V);
+    # Determine frequency coordinates
+    my $num_elements = nelem($U);
+    say "Found $num_elements elements.";
+    $f = $U->xlinvals(-($num_elements/2-1)/$num_elements/$xres,1/2/$xres)->rotate(-($num_elements/2 -1));
+
+    my $φdiff = atan2($V,$U);
+    my $timelag = $φdiff/($twoPI*$f);
+
+    wcols $f,$timelag,"analyses/tables/tophat_fft${$_}.tab";
+
 }
-close $tophattab;
-
-# Transform z(x) to z(1/x) = u(1/x) + iv(1/x)
-fft($u,$v);
-# Determine frequency coordinates
-my $num_elements = nelem($u);
-say "Found $num_elements elements.";
-$f = $u->xlinvals(-($num_elements/2-1)/$num_elements/$Δt,1/2/$Δt)->rotate(-($num_elements/2 -1));
-
-my $φdiff = atan2($v,$u);
-my $timelag = $φdiff/($twoPI*$f);
-
-wcols $f,$timelag,'analyses/tables/tophat_fft.tab';
-
 
 sub tophat {
     (my $mean,my $width) = @_;