From de231e730f9eae5b19d7739780174027f59b4040 Mon Sep 17 00:00:00 2001 From: caes Date: Thu, 22 Jun 2017 06:10:46 -0400 Subject: [PATCH] init --- scripts/extract_tables.pl | 197 ++++++++++++++++++++++++++++++++++++ scripts/plot.sh | 57 +++++++++++ scripts/plot_psdlag_grid.sh | 76 ++++++++++++++ scripts/propagate_tables.sh | 29 ++++++ scripts/psdlag_4bin.py | 104 +++++++++++++++++++ 5 files changed, 463 insertions(+) create mode 100755 scripts/extract_tables.pl create mode 100755 scripts/plot.sh create mode 100644 scripts/plot_psdlag_grid.sh create mode 100755 scripts/propagate_tables.sh create mode 100644 scripts/psdlag_4bin.py diff --git a/scripts/extract_tables.pl b/scripts/extract_tables.pl new file mode 100755 index 0000000..12d4752 --- /dev/null +++ b/scripts/extract_tables.pl @@ -0,0 +1,197 @@ +#!/usr/bin/env perl + +# This extracts the gnuplot-ready tables from a psdlag analysis +# and saves them to temporary files. propogate_tables.sh uses this +# script to create the tables it saves to analyses/tables. + +use utf8; +use Encode qw(encode decode); +use feature 'say'; +use locale; +use Switch; + +use constant PI => 4 * atan2(1, 1); + +# Enables various levels of output. +our $verbose=0; +our $debug=0; + +if ($debug) { $verbose=1; } + + +# This section locates the output data of interest in a +# psdlab output file. +if (!${^UTF8LOCALE}) { + say encode($charset,"You are not using UTF-8 encoding. :("); +} +my $charset=$ENV{LANG}; + + +# @ARGV[0] takes the echo wavelength for use in filename + +$echo_λ = @ARGV[0]; + +# This program attempts to fill function_bin with the tabulated +# PSDs and time lags. +our %function_bin = (); + +open $freq_file,'<',"freq.out" or die $!; +open $ref_psd_file,'<',"ref_psd.out" or die $!; +open $echo_psd_file,'<',"echo_psd.out" or die $!; +open $crsspctrm_file,'<',"crsspctrm.out" or die $!; +open $timelag_file,'<',"timelag.out" or die $!; + +=pod + + This section collects the various quantities. + +=cut + +@bin_bounds = split /\s/,<$freq_file>; +@ref_psd = split /\s/,<$ref_psd_file>; +@ref_psd_σ = split /\s/,<$ref_psd_file>; +@echo_psd = split /\s/,<$echo_psd_file>; +@echo_psd_σ = split /\s/,<$echo_psd_file>; +@crsspctrm_psd = split /\s/,<$crsspctrm_file>; +@crsspctrm_psd_σ = split /\s/,<$crsspctrm_file>; +@timelag = split /\s/,<$timelag_file>; +@timelag_σ = split /\s/,<$timelag_file>; + +if ($debug) { + print encode($charset,"Found bin boundaries: "); + foreach (@bin_bounds) {print encode($charset,"$_ ");} + say encode($charset," "); +} +my $count = 0; +my $upper_bound = 0; +my $lower_bound = 0; +foreach (@bin_bounds) { + $lower_bound = $upper_bound; + $upper_bound = $_; + if ($lower_bound == 0) {next} + my $μ = ($upper_bound + $lower_bound)/2; + my $Δ = ($upper_bound - $lower_bound)/2; + #say ($μ,":",$Δ); + # push(@freq_coords_mean,$μ); + # push(@freq_coords_σ,$Δ); + $function_bin{$μ} = {"Δ" => $Δ, + "ref_PSD_μ" => $ref_psd[$count], + "ref_PSD_σ" => $ref_psd_σ[$count], + "echo_PSD_μ" => $echo_psd[$count], + "echo_PSD_σ" => $echo_psd_σ[$count], + "crsspctrm_μ" => $crosssp_psd[$count], + "crsspctrm_σ" => $crosssp_psd_σ[$count], + "timelag_μ" => $timelag[$count], + "timelag_σ" => $timelag_σ[$count]}; + # $function_bin{$μ}{"φdiff_μ"} = $μ; + # $function_bin{$μ}{"φdiff_σ"} = $σ; + # $μ = $μ/(2*PI*$μ); + # $σ = $σ/(2*PI*$μ); + $count = $count + 1; +} + +close $freq_file; +close $ref_psd_file; +close $echo_psd_file; +close $crsspctrm_file; +close $timelag_file; + +$numbins = keys %function_bin; +if ($debug) { say encode($charset, + "$numbins frequency bins captured in output."); +} + +if($verbose) { + say encode($charset,"freq μ (lin) freq σ ". + "ref_PSD_μ ref_PSD_σ ". + "echo_PSD_μ echo_PSD_σ ". + "timelag_μ timelag_σ"); + + foreach (sort { $a <=> $b } keys %function_bin) { + say encode($charset, + sprintf("%f %f %f %f %f ". + "%f %f %f ", + $_, + $function_bin{$_}{"Δ"}, + $function_bin{$_}{"ref_PSD_μ"}, + $function_bin{$_}{"ref_PSD_σ"}, + $function_bin{$_}{"echo_PSD_μ"}, + $function_bin{$_}{"echo_PSD_σ"}, + $function_bin{$_}{"timelag_μ"}, + $function_bin{$_}{"timelag_σ"} + )); + } +} + + +say encode($charset,""); + +if($debug) { + while (each %function_bin ) { + print encode($charset,$_ . ": "); + while ( my ($key,$value) = each %{$function_bin{$_}} ) { + print encode($charset,$key . " => " . $value . "; "); + } + say encode($charset,""); + } +} + +close($outputfile); + +open($datafile,'>',"tmp.refPSD") or die $!; +while( each %function_bin) { + say $datafile + $_ . " " . + $function_bin{$_}{"ref_PSD_μ"} . " " . + $function_bin{$_}{"Δ"} . " " . + $function_bin{$_}{"ref_PSD_σ"}; +} +close($datafile); + +open($datafile,'>',"tmp.echoPSD") or die $!; +while( each %function_bin) { + say $datafile + $_ . " " . + $function_bin{$_}{"echo_PSD_μ"} . " " . + $function_bin{$_}{"Δ"} . " " . + $function_bin{$_}{"echo_PSD_σ"}; +} +close($datafile); + +open($datafile,'>',"tmp.crsspctrmPSD") or die $!; +while( each %function_bin) { + say $datafile + $_ . " " . + $function_bin{$_}{"crsspctrm_PSD_μ"} . " " . + $function_bin{$_}{"Δ"} . " " . + $function_bin{$_}{"crsspctrm_PSD_σ"}; +} +close($datafile); + +open($datafile,'>',"tmp.timelag") or die $!; +while( each %function_bin) { + say $datafile + $_ . " " . + $function_bin{$_}{"timelag_μ"} . " " . + $function_bin{$_}{"Δ"} . " " . + $function_bin{$_}{"timelag_σ"}; +} +close($datafile); + + +#open($datafile,'>',"tmp.echoPSD") or die $!; + + +open($datafile,'>',"cackett_psdlag_" . $echo_λ . ".tab") or die $!; +say $datafile encode($charset, + sprintf("#freqmin freqmax psd ". + "psd error lag lag error " + )); +foreach ( sort {$a <=> $b} keys %function_bin ) { + say $datafile encode($charset, + sprintf("%e %e %e %e %e %e", + ($_ - $function_bin{$_}{"Δ"}), + ($_ + $function_bin{$_}{"Δ"}) + )); +} +close($datafile) \ No newline at end of file diff --git a/scripts/plot.sh b/scripts/plot.sh new file mode 100755 index 0000000..43bf1c4 --- /dev/null +++ b/scripts/plot.sh @@ -0,0 +1,57 @@ +#!/usr/bin/env bash + +# This metascript uses the available plotting scripts to produce a list of +# document-rady plots. + +# Only analyses with this error type will be represented in the atlas +errtype=$(cat err_type) +ref_band="1367Å" + +case $1 in + "PSD"|"psd"|"PSDs"|"PSDS"|"psds") + gnuplot_file=psd_atlas.gp + scripts/propagate_tables.sh + gnuplot_input=$(cat scripts/templates/${gnuplot_file}|perl -pe 's|\n|␤|g') + for tabfile in analyses/tables/PSD_*${errtype}*.tab; + do + echo_band=$(basename $tabfile| + sed 's|PSD[_ ]\(.\{5\}\)[_ ]{[^_ ]*}.tab|\1|') + if [[ "$echo_band" == "$ref_band" ]] ; then continue; fi + gnuplot_input_edit=$(echo "$gnuplot_input"| + sed "s|%FILE%|$tabfile|"| + sed "s|%LABEL%|$echo_band|") + gnuplot_input="${gnuplot_input_edit}" + done + echo "$gnuplot_input"|perl -pe 's|␤|\n|g' > ${gnuplot_file} + gnuplot $gnuplot_file + ;; + + "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 analyses/tables/timelag_*${errtype}*.tab; + do + ref_band_extracted=$(basename $tabfile|sed 's|timelag_\([^≺]*\)[_ ]≺[_ ][^≺_ ]*[_ ]{[^_ ]*}.tab|\1|') + echo_band=$(basename $tabfile|sed 's|timelag_[^≺]*[_ ]≺[_ ]\([^≺_ ]*\)[_ ]{[^_ ]*}.tab|\1|') + if [[ "$echo_band" == "$ref_band" ]] ; then continue; fi + gnuplot_input_edit=$(echo "$gnuplot_input"| + sed "s|%FILE%|$tabfile|"| + sed "s|%LABEL%|$echo_band|") + gnuplot_input="${gnuplot_input_edit}" + done + echo "$gnuplot_input"|perl -pe 's|␤|\n|g' > ${gnuplot_file} + gnuplot $gnuplot_file + ;; + + "tophat"|"th") + mkdir -p analyses/tables/ + scripts/tophat_fft.pl + gnuplot scripts/templates/tophat_freqdomain.gp + gnuplot scripts/templates/tophat_timedomain.gp + ;; + + *) + echo "Did not understand plot type." + ;; +esac diff --git a/scripts/plot_psdlag_grid.sh b/scripts/plot_psdlag_grid.sh new file mode 100644 index 0000000..aacce58 --- /dev/null +++ b/scripts/plot_psdlag_grid.sh @@ -0,0 +1,76 @@ +#!/usr/bin/env bash + +ywindow_of () { + max=0 + min=-5 + while read line +# for line in $(cat $1) + do + read xval yval xerr yerr <<< $line + echo $xval $yval + if (bc <<< "$yval > $max"); then max=$yval; fi + if (bc <<< "$yval < $min"); then min=$yval; fi +# done + done < $1 + echo $(bc <<< "$min - 0.5") $(bc <<< "$max + 0.5") +} + +mkdir -p analyses/tables +mkdir -p analyses/plots + + +# Restore the proper naming conventions for any files stripped +# of their unicode characters. +scripts/fix_thor_names.sh + +# Create tables and plot the power spectra and time lags. +for table in analyses/* +do + # Grab and determine labels of analyses, skip if over the same band. + ref_band=$(basename $table|sed 's|\([^≺]*\)[_ ]≺[_ ][^≺_ ]*[_ ]{[^_ ]*}|\1|') + echo_band=$(basename $table|sed 's|[^≺]*[_ ]≺[_ ]\([^≺_ ]*\)[_ ]{[^_ ]*}|\1|') + err_str=$(basename $table|sed 's|[^≺]*[_ ]≺[_ ][^≺_ ]*[_ ]{[^_ ]*;\(σ∊[CLM][MFC]\)}|\1|') + + # Prepare files + echo "Plotting PSD and time lags for $echo_band, referred to ${ref_band}." + echoPSD_tabfile=analyses/tables/PSD_${echo_band}_\{${err_str}\}.tab + refPSD_tabfile=analyses/tables/PSD_${ref_band}_\{${err_str}\}.tab + timelag_tabfile=analyses/tables/timelag_${ref_band}_≺_${echo_band}_\{${err_str}\}.tab + PSD_plotfile=analyses/plots/PSD_${echo_band}_\{${err_str}\}.pdf + timelag_plotfile=analyses/plots/timelag_${ref_band}_≺_${echo_band}_\{${err_str}\}.pdf + + # Output curves to temporary files using perl script, move tables to + # permanent location. This just assumes there are no conflicts. + scripts/extract_tables.pl $table > /dev/null + mv tmp.echoPSD $echoPSD_tabfile + mv tmp.refPSD $refPSD_tabfile + mv tmp.timelag $timelag_tabfile + + + # Plot PSD and save using gnuplot +# read ymin ymax <<< $(ywindow_of $timelag_tabfile) +# echo $ymin $ymax + cat scripts/templates/psd.gp| + sed "s|%TITLE%|Power Spectrum for Lightcurves $echo_band and $ref_band|"| + sed "s|%SUBTITLE%|as reported by Fausnaugh et. al, STORM III, 2016|"| + sed "s|%FILE1%|$echoPSD_tabfile|"| + sed "s|%LABEL1%|${echo_band}|"| + #sed "s|%FILE2%|$refPSD_tabfile|"| + #sed "s|%LABEL2%|${ref_band}|"| + sed "s|%YMIN%|$ymin|"|sed "s|%YMAX|$ymax|"| + sed "s|%OUTPUTFILE%|$PSD_plotfile|" > tmp.gp + gnuplot tmp.gp + + # Plot time lags and save using gnuplot +# read ymin ymax <<< $(ywindow_of $timelag_tabfile) + cat scripts/templates/timelag.gp| + sed "s|%TITLE%|Time Lag for Lightcurve $echo_band relative to $ref_band|"| + sed "s|%SUBTITLE%|as reported by Fausnaugh et. al, STORM III, 2016|"| + sed "s|%FILE1%|$timelag_tabfile|"| + sed "s|%LABEL1%|${echo_band}|"| + sed "s|%YMIN%|$ymin|"|sed "s|%YMAX|$ymax|"| + sed "s|%OUTPUTFILE%|$timelag_plotfile|" > tmp.gp + gnuplot tmp.gp +done + +rm tmp.* diff --git a/scripts/propagate_tables.sh b/scripts/propagate_tables.sh new file mode 100755 index 0000000..8cebd43 --- /dev/null +++ b/scripts/propagate_tables.sh @@ -0,0 +1,29 @@ +#!/usr/bin/env bash + +mkdir -p data + +ref_band="1367A" + +for lightcurve in data/lc/*.lc +do + + echo_band=$(basename $lightcurve|sed 's|[^≺]*[_ ]≺[_ ]\([^≺_ ]*\)[_ ]{[^_ ]*}|\1|') + if [[ $ref_band == $echo_band ]]; then continue; fi + err_str=$(basename $lightcurve|sed 's|[^≺]*[_ ]≺[_ ][^≺_ ]*[_ ]{[^_ ]*;\(σ∊[CLM][MFC]\)}|\1|') + +echo "Tabling PSD and time lags referred to ${ref_band} for $echo_band{${err_str}}." + + # Propagate tables into analyses/tables + echoPSD_tabfile=analyses/tables/PSD_${echo_band}_\{${err_str}\}.tab + refPSD_tabfile=analyses/tables/PSD_${ref_band}_\{${err_str}\}.tab + timelag_tabfile=analyses/tables/timelag_${ref_band}_≺_${echo_band}_\{${err_str}\}.tab + + # Output curves to temporary files using perl script, move tables to + # permanent location. This just assumes there are no conflicts. + scripts/extract_tables.pl $lightcurve > /dev/null + mv tmp.echoPSD $echoPSD_tabfile + mv tmp.refPSD $refPSD_tabfile + mv tmp.timelag $timelag_tabfile +done + +rm tmp.* \ No newline at end of file diff --git a/scripts/psdlag_4bin.py b/scripts/psdlag_4bin.py new file mode 100644 index 0000000..3bf3d3d --- /dev/null +++ b/scripts/psdlag_4bin.py @@ -0,0 +1,104 @@ +import numpy as np +from scipy.stats import norm +from scipy.stats import lognorm +import sys +import getopt +sys.path.insert(1,"/usr/local/science/clag-agn/data/") +import clag +import matplotlib +# %pylab inline + + +ref_file="data/lc/1367A.lc" +echo_file="data/lc/2246A.lc" + +#ref_file = str(sys.argv[1]) +#echo_file = str(sys.argv[2]) + + +dt = 0.01 +t1, l1, l1e = np.loadtxt(ref_file).T +# errorbar(t1, l1, yerr=l1e, fmt='o') + + + +fqL = np.array([0.0049999999, 0.018619375, 0.044733049, 0.069336227, 0.10747115, 0.16658029, + 0.25819945, 0.40020915, 0.62032418]) +fqL = np.array([0.0049999999, 0.044733049, 0.10747115, + 0.25819945, 0.62032418]) +fqL = np.array([0.0049999999, 0.018619375, 0.069336227, 0.10747115, 0.62032418]) +fqL = np.logspace(np.log10(0.005),np.log10(0.6),5) +nfq = len(fqL) - 1 +fqd = 10**(np.log10( (fqL[:-1]*fqL[1:]) )/2.) + + + + + +P1 = clag.clag('psd10r', [t1], [l1], [l1e], dt, fqL) +p1 = np.ones(nfq) +p1, p1e = clag.optimize(P1, p1) + +p1, p1e = clag.errors(P1, p1, p1e) + +# xscale('log'); ylim(-4,2) +# errorbar(fqd, p1, yerr=p1e, fmt='o', ms=10, color="black") + +ref_psd = p1 +ref_psd_err = p1e + + + +t2, l2, l2e = np.loadtxt(echo_file).T +# errorbar(t1, l1, yerr=l1e, fmt='o', color="green") +# errorbar(t2, l2, yerr=l2e, fmt='o', color="black") + + +P2 = clag.clag('psd10r', [t2], [l2], [l2e], dt, fqL) +p2 = np.ones(nfq) +p2, p2e = clag.optimize(P2, p2) + +p2, p2e = clag.errors(P2, p2, p2e) + +# xscale('log'); ylim(-6,2) +# errorbar(fqd, p1, yerr=p1e, fmt='o', ms=10, color="green") +# errorbar(fqd, p2, yerr=p2e, fmt='o', ms=10, color="black") +echo_psd = p2 +echo_psd_err = p2e + + + + +Cx = clag.clag('cxd10r', [[t1,t2]], [[l1,l2]], [[l1e,l2e]], dt, fqL, p1, p2) +p = np.concatenate( ((p1+p2)*0.5-0.3,p1*0+0.1) ) # a good starting point generally +p, pe = clag.optimize(Cx, p) + +phi, phie = p[nfq:], pe[nfq:] +lag, lage = phi/(2*np.pi*fqd), phie/(2*np.pi*fqd) +cx, cxe = p[:nfq], pe[:nfq] + +cross_spectrm = cx +cross_spectrm_err = cxe + +# xscale('log'); ylim(-2,1) +# errorbar(fqd, lag, yerr=lage, fmt='o', ms=10,color="black") + + + + +s, loc, scale = lognorm.fit(lag,loc=.01) +# xscale('log'); ylim(-4,1.5) +# errorbar(fqd, lag, yerr=lage, fmt='o', ms=10,color="black") +##plot(fqd,norm.pdf(fqd,mu,sigma)) +#plot(fqd,lognorm.pdf(fqd,s,loc,scale)) +# mu,sigma + + + +#plot(ifft(lag)) + +np.savetxt("freq.out",fqL.reshape((-1,len(fqL)))) +np.savetxt("ref_psd.out",[ref_psd,ref_psd_err]) +np.savetxt("echo_psd.out",[echo_psd,echo_psd_err]) +np.savetxt("crsspctrm.out",[cross_spectrm,cross_spectrm_err]) +np.savetxt("timelag.out",[lag,lage])