diff --git a/pres/img/tophat_freqdomain.pdf b/pres/img/tophat_freqdomain.pdf index 1760a0c..dbff5db 100644 Binary files a/pres/img/tophat_freqdomain.pdf and b/pres/img/tophat_freqdomain.pdf differ diff --git a/scripts/templates/tophat_freqdomain.gp b/scripts/templates/tophat_freqdomain.gp index b5490fc..7eae0ae 100644 --- a/scripts/templates/tophat_freqdomain.gp +++ b/scripts/templates/tophat_freqdomain.gp @@ -1,6 +1,6 @@ set terminal pdf set termopt enhanced -set output "pres/img/tophat_freqdomain.pdf" +set output "tophat_freqdomain.pdf" set xrange [0.005:0.620] set logscale x @@ -8,7 +8,7 @@ set xlabel "Temporal Frequency [days^{-1}]" font "Times,24" set xtics font 'Times,16' set mxtics -set yrange [-11:11] +set yrange [-30:30] set ylabel "Time Delay [days]" font "Times,24" set ytics font 'Times,16' set mytics @@ -17,6 +17,6 @@ set mytics unset key -plot "analyses/tables/tophat_fft1.tab" using 1:(-1*$2) with lines lc rgb "red", \ - "analyses/tables/tophat_fft2.tab" using 1:(-1*$2) with lines lc rgb "green", \ - "analyses/tables/tophat_fft3.tab" using 1:(-1*$2) with lines lc rgb "blue" \ No newline at end of file +plot "analyses/tables/tophat_fft1.tab" using 1:2 with lines lc rgb "red", \ + "analyses/tables/tophat_fft2.tab" using 1:2 with lines lc rgb "green", \ + "analyses/tables/tophat_fft3.tab" using 1:2 with lines lc rgb "blue" \ No newline at end of file diff --git a/scripts/templates/tophat_timedomain.gp b/scripts/templates/tophat_timedomain.gp index 3aefb82..d4bf48c 100644 --- a/scripts/templates/tophat_timedomain.gp +++ b/scripts/templates/tophat_timedomain.gp @@ -1,7 +1,7 @@ set terminal pdf -set output "pres/img/tophat_timedomain.pdf" +set output "tophat_timedomain.pdf" -set xrange [0:25] +set xrange [0:35] set xlabel "Time [days]" font "Times,24" set xtics font 'Times,16' set mxtics diff --git a/scripts/tophat_fft.pl b/scripts/tophat_fft.pl index b4cdf66..297bd13 100755 --- a/scripts/tophat_fft.pl +++ b/scripts/tophat_fft.pl @@ -6,55 +6,66 @@ use PDL; use PDL::FFT; use PDL::IO::Misc; -my $twoPI = 2*3.1415926539; +my $_2π = 2*3.1415926539; # space parameters -our $xres=.001; +our $xres= 0.1; our $xmin = 0; -our $xmax = 250; -our $xbins = ($xmax-$xmin)/$xres; +our $xmax = 2500; +our $xbins = ($xmax-$xmin)/$xres + 1; + +say "Using $xbins bins"; # tophat parameters -my $μ1=15; +my $μ1=7; my $Δ1=5; # full width -my $μ2=15; -my $Δ2=5; # full width -my $μ3=15; -my $Δ3=5; # full width -say "Creating $xbins bins"; +my @tophat_list = ([$μ1,$Δ1],[28,8],[15,11]); -# Complex Number z(x) = u(x) + iv(x) -my @th_list=tophat($μ1,$Δ1); -my $u = pdl(@th_list); -my $v = zeroes($u); +our $tophat_count = 0; +foreach my $pars (@tophat_list) { + $tophat_count++; + # Capture tophat into piddle + # Complex Number z(x) = u(x) + iv(x) + my @tophat=tophat(@$pars[0],@$pars[1]); + my $u = pdl(@tophat); + my $v = zeroes($u); -foreach (('1','2','3')) { + $x_coords = $u->xlinvals($xmin,$xmax); - # Output time-domain tophat - open $tophattab, '>', "analyses/tables/tophat${$_}.tab"; - for ($i=0; $i < $xbins; $i++) { - my $x_coord = ($xmax-$xmin)*($i/$xbins) + $xmin; - say $tophattab "$x_coord $th_list[$i]"; - } - close $tophattab; + # Output x-domain tophat + wcols $x_coords,$u,"analyses/tables/tophat${tophat_count}.tab"; - # Transform z(x) to Z(1/x) = U(1/x) + iV(1/x) + # 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 = $U->xlinvals(-($num_elements/2-1)/$num_elements/$xres,1/2/$xres)->rotate(-($num_elements/2 -1)); - say "$f[1]"; + if ($num_elements % 2 == 0) { + # even number of bins + $f = $U->xlinvals( + -(${num_elements}/2-1)/${num_elements}/${xres}, 1/2/${xres} + )->rotate(-(${num_elements}/2 -1)); + } + else { + #odd number of bins + $f = $U->xlinvals( + -(${num_elements}/2-0.5)/${num_elements}/${xres}, + (${num_elements}/2-0.5)/${num_elements}/${xres} + )->rotate(-(${num_elements}-1)/2); + } + + # Calculate x offset my $φdiff = atan2($V,$U); - my $timelag = $φdiff/($twoPI*$f); + my $offset = $φdiff/($_2π*$f); - wcols $f,$timelag,"analyses/tables/tophat_fft${$_}.tab"; + wcols $f,$offset,"analyses/tables/tophat_fft${tophat_count}.tab"; } @@ -62,10 +73,9 @@ sub tophat { (my $mean,my $width) = @_; my $halfwidth = $width/2; my @vals = (); - for ($i=0; $i < $xbins; $i++) { - my $x_coord = ($xmax-$xmin)*($i/$xbins) + $xmin; - if ($x_coord >= ($mean - $halfwidth ) && $x_coord <= ($mean + $halfwidth)) { push @vals, 1/$width; } - else { push @vals, 0; } + 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 @vals; } \ No newline at end of file