Science/psdlag-agn
1
0
Fork 0
mirror of https://asciireactor.com/otho/psdlag-agn.git synced 2024-11-22 02:15:07 +00:00
Code Issues Packages Projects Releases Wiki Activity

finalized tophat fft program

Browse Source
This commit is contained in:
caes 2016-08-13 22:45:00 -04:00
parent 7cbf4db15c
commit 6ffe2c92bc
4 changed files with 48 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
pres/img/tophat_freqdomain.pdf
View File

Binary file not shown.

10
scripts/templates/tophat_freqdomain.gp
View File

@ -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"
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"

4
scripts/templates/tophat_timedomain.gp
View File

@ -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

70
scripts/tophat_fft.pl
View File

@ -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]);
our $tophat_count = 0;
foreach my $pars (@tophat_list) {
$tophat_count++;
# Capture tophat into piddle
# Complex Number z(x) = u(x) + iv(x)
my @th_list=tophat($μ1,$Δ1);
my $u = pdl(@th_list);
my @tophat=tophat(@$pars[0],@$pars[1]);
my $u = pdl(@tophat);
my $v = zeroes($u);
$x_coords = $u->xlinvals($xmin,$xmax);
foreach (('1','2','3')) {
# Output x-domain tophat
wcols $x_coords,$u,"analyses/tables/tophat${tophat_count}.tab";
# 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;
# 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;
}