From c1a768dd8c8e08eb284f4953f948d1f8ac323c1e Mon Sep 17 00:00:00 2001
From: caes <caes>
Date: Thu, 22 Jun 2017 07:18:00 -0400
Subject: [PATCH] added tophat stuff

---
 scripts/plot.sh       |  9 +++--
 scripts/tophat_fft.pl | 85 +++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 90 insertions(+), 4 deletions(-)
 create mode 100755 scripts/tophat_fft.pl

diff --git a/scripts/plot.sh b/scripts/plot.sh
index 4cf8767..1c16c39 100755
--- a/scripts/plot.sh
+++ b/scripts/plot.sh
@@ -27,11 +27,12 @@ case $1 in
     "lags"|"lag"|"delay"|"delays")
         gnuplot_file=timelag_atlas.gp
         gnuplot_input=$(cat scripts/templates/${gnuplot_file}|perl -pe 's|\n|␤|g')
-        scripts/propagate_tables.sh
-        for tabfile in data/tables/timelag_*${errtype}*.tab;
+        for tabfile in data/tables/lag_*.tab;
         do
-            ref_band_extracted=$(basename $tabfile|sed 's|timelag_\([^≺]*\)[_ ]≺[_ ][^≺_ ]*[_ ]{[^_ ]*}.tab|\1|')
-            echo_band=$(basename $tabfile|sed 's|timelag_[^≺]*[_ ]≺[_ ]\([^≺_ ]*\)[_ ]{[^_ ]*}.tab|\1|')
+            ref_band_extracted=$(basename $tabfile|
+                                    sed 's|lag_\([0-9]\{4\}A\)_[0-9]\{4\}A.tab|\1|')
+            echo_band=$(basename $tabfile|
+                        sed 's|lag_[0-9]\{4\}A_\([0-9]\{4\}A\).tab|\1|')
             if [[ "$echo_band" == "$ref_band" ]] ; then continue; fi
             gnuplot_input_edit=$(echo "$gnuplot_input"|
                                     sed "s|%FILE%|$tabfile|"|
diff --git a/scripts/tophat_fft.pl b/scripts/tophat_fft.pl
new file mode 100755
index 0000000..c2f2979
--- /dev/null
+++ b/scripts/tophat_fft.pl
@@ -0,0 +1,85 @@
+#!/usr/bin/env perl
+
+use feature say;
+use utf8;
+use PDL;
+use PDL::FFT;
+use PDL::IO::Misc;
+
+my $_2π = 2*3.1415926539;
+
+# space parameters
+our $xres= 0.1;
+our $xmin = 0;
+our $xmax = 2500;
+our $xbins = ($xmax-$xmin)/$xres + 1;
+
+say "Using $xbins bins";
+
+# tophat parameters
+my $μ1=7;
+my $Δ1=5; # full width
+
+my @tophat_list = ([$μ1,$Δ1],[28,8],[15,11]);
+
+our $tophat_count = 0;
+foreach (@tophat_list) {
+    $tophat_count++;
+
+    # Capture tophat into piddle
+    # Complex Number z(x) = u(x) + iv(x)
+    my $u = tophat(@$_[0],@$_[1]);
+    my $v = zeroes($u);
+
+    $x_coords = $u->xlinvals($xmin,$xmax);
+
+    # Output x-domain tophat
+    wcols $x_coords,$u,"analyses/tables/tophat${tophat_count}.tab";
+
+    # Transform z(x) to Z(x⁻¹) = U(x⁻¹) + iV(x⁻¹)
+    my $U = $u;
+    my $V = $v;
+    fft($U,$V);
+
+    # Determine frequency coordinates
+    my $num_elements = nelem($U);
+    say "Found $num_elements elements.";
+
+    $f = ($num_elements % 2 == 0) ?
+        $U->xlinvals(
+            -(${num_elements}/2-1)/${num_elements}/${xres},
+            1/2/${xres}
+        )->rotate(-(${num_elements}/2 -1))
+    :
+        $U->xlinvals(
+            -(${num_elements}/2-0.5)/${num_elements}/${xres},
+            (${num_elements}/2-0.5)/${num_elements}/${xres}
+        )->rotate(-(${num_elements}-1)/2);
+
+    # Output frequency-domain tophat
+    wcols $f, $U, $V, "fft${tophat_count}.tab";
+
+    # Calculate x offset
+    #   This currently multiplies the imaginary component by -1,
+    #   and I really need to figure out why this is necessary for
+    #   proper output.
+    my $φdiff = atan2(-1*$V,$U);
+
+    wcols $f,$φdiff,"analyses/tables/tophat_φdiff${tophat_count}.tab";
+
+    my $offset = $φdiff/($_2π*$f);
+
+    # Output frequency-domain time delay for given tophat
+    wcols $f,$offset,"analyses/tables/tophat_fft${tophat_count}.tab";
+}
+
+sub tophat {
+    (my $mean,my $width) = @_;
+    my $halfwidth = $width/2;
+    my @vals = ();
+    for (my $x_coord = $xmin; $x_coord <= $xmax; $x_coord += $xres) {
+        push @vals, ($x_coord >= ($mean - $halfwidth) &&
+                    $x_coord <= ($mean + $halfwidth)) ? 1/$width : 0;
+    }
+    return pdl(@vals);
+}
\ No newline at end of file