init

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)

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

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.*

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.*

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])