#! /bin/bash 
# Last edited on 2012-12-28 02:25:31 by stolfilocal

name="$1"; shift

# Reformat the output table to make it easier to plot the histogram:
# Also try to fit a Gaussian {p/n = exp(A + B*f^2)} to the plot:

cat ${name}.txt \
  | gawk -v tmp=${tmp} \
      ' BEGIN { 
          split("",tfmin); split("",tfmax); split("", tfr); split("", tpw);
          s11 = 0; sq1 = 0; sqq = 0; sp1 = 0; spq = 0;
          nb = 0;
          filea = (tmp "-a.txt"); # For actual spectra
          fileb = (tmp "-b.txt"); # For fitted Gaussian model
        } 
        /[0-9]/ { 
          fmin = $1; fmax = $2; fr = $3; n = $4; p = $5; 
          pw = p/n;
          if ((fr >= 0.1) && (fr <= 0.5)) 
            { # Accumulate inner products for least squares:
              y = log(pw); u = 1; v = fr*fr; 
              suu += u*u; suv += u*v; svv += v*v; 
              syu += y*u; syv += y*v;
            }
          tfmin[nb] = fmin; tfmax[nb] = fmax; tfr[nb] = fr;  tpw[nb] = pw;
          nb++;
        }
        END {
          # Compute quadratic coefs
          det = suu*svv - suv*suv;
          A = (syu*svv - syv*suv)/det;
          B = (suu*syv - suv*syu)/det;
          printf "# A = %010.6f  B = %010.6f\n", A, B > "/dev/stderr"
          
          # Write plot files
          ofr = 0;
          for (i = 0; i < nb; i++)
            { if (tfmin[i] < ofr) { printf "freq out of order\n" > "/dev/stderr"; exit(1); }
              prt(tfmin[i], tpw[i], filea);
              prt(tfmax[i], tpw[i], filea);
              printf "\n" > filea;
              ofr = tfmax[i]; 
              xfr = tfr[i];
              if ((xfr >= 0.1) && (xfr <= 0.5))
                { xpw = exp(A + B*xfr*xfr); prt(xfr, xpw, fileb); }
            }
          close(filea);
          close(fileb);
        }
        function prt(fr, pw, file) { printf "%24.16e  %24.16e\n", fr, pw > file; }
      '
       
# Plot the histograms:
gnuplot <<EOF

set term postscript eps color lw 2 "TimesRoman" 24
set size 1,2
set nokey

set yrange [+0.00000001:]
set logscale y
set format y '%.0e'
set title "binned spectrum - power/terms"
set output "${name}-r.eps"
plot \
  "${tmp}-a.txt"  using 1:2 with lines lw 2 lc rgb '#228800', \
  "${tmp}-b.txt"  using 1:2 with points pt 7 ps 1 lc rgb '#ff2200'

EOF