added gaussian.awk based on hamming script;

added pwr(x,p) function to sterling_approx;
minor clean up all around;

1 files changed, 13 insertions, 13 deletions

diff --git a/lin_reg2.awk b/lin_reg2.awk
index 8939947..f06d890 100644
--- a/lin_reg2.awk
+++ b/lin_reg2.awk
@@ -22,27 +22,26 @@ NF > 0 {
 
     ### iterate over columns
     for (y=1; y<=nf_max; y++) {
+        if ($y == header[n])
+            continue
         if ($y ~ number) {
 
             ### mean
             count[y] += 1
             sum[y] += $y
             sum2[y] += $y*$y
-            mean[y] = sum[y]/count[y]
-
-            ### difference from the mean
-            delta[y] = $y - mean[y]
-            sum_delta[y] += delta[y]
-            sum_delta2[y] += delta[y]*delta[y]
+            delta0[y] = $y - mean[y]
+            mean[y] += delta0[y]/count[y]
+            delta1[y] = $y - mean[y]
+            sum_delta[y] += delta1[y]
+            sum_delta2[y] += delta0[y]*delta1[y]
 
-            ### sample variance
-            (count[y] > 1) ? var[y] = sum_delta2[y]/(count[y] - 1) : var[y] = ""
 
-            # x = row, y = col
+            # x = row, y = col, trendline: y = A + Bx
             for (x=1; x<=nf_max; x++) {
                 count[x,y] += 1
                 sum_xy[x,y] += $x*$y
-                sum_delta_xy[x,y] += delta[x]*delta[y]
+                sum_delta_xy[x,y] += delta0[x]*delta1[y]
 
                 # correlation
                 r_den[x,y] = sqrt(sum_delta2[x]*sum_delta2[y])
@@ -58,7 +57,7 @@ NF > 0 {
                     b[x,y] = 1
                 }
 
-                ### error estimate
+                # error estimate
                 err_den[x,y] = count[x,y]*(count[x,y] - 2)
                 if (count[x,y] > 2) {
                     err[x,y] = $y - (a[x,y] + b[x,y]*$x)
@@ -82,8 +81,9 @@ END {
         for (x=1; x<=nf_max; x++) {
             if (x != y && r[x,y]) {
                 printf(OFMT OFS "(%s)" OFS " = (" OFMT " +/- " OFMT ")(%s)" OFS " + (" OFMT " +/- " OFMT ")" ORS,
-                    10.0*log(r[x,y]*r[x,y])/log(10), header[y], b[x,y], b_err[x,y], header[x],
-                    a[x,y], a_err[x,y])
+                    10.0*log(r[x,y]*r[x,y])/log(10.0),
+                    header[y], b[x,y], b_err[x,y],
+                    header[x], a[x,y], a_err[x,y])
             }
         }
     }