restructured conditional prints in diff, added parabolic focus solution to quad_reg output

1 files changed, 19 insertions, 8 deletions

diff --git a/quad_reg.awk b/quad_reg.awk
index 7d9304a..4267412 100644
--- a/quad_reg.awk
+++ b/quad_reg.awk
@@ -4,7 +4,6 @@
 # quadratic regression along columns
 
 BEGIN {
-    OFMT="%.9g"
     sign = "[+-]?"
     decimal = "[0-9]+[.]?[0-9]*"
     fraction = "[.][0-9]*"
@@ -23,6 +22,8 @@ NF > 0 {
 
     ### iterate over columns
     for (y=1; y<=nf_max; y++) {
+        if ($y == header[n])
+            continue
         if ($y ~ number) {
 
             ### mean
@@ -77,7 +78,7 @@ NF > 0 {
                 }
                 else {
                     c[x,y] = 0
-                    b[x,y] = 0
+                    b[x,y] = 1
                 }
                 a[x,y] = mean[y] - b[x,y]*mean[x] - c[x,y]*mean2[x]
 
@@ -88,16 +89,24 @@ NF > 0 {
                 # correlation
                 sum_delta2[y] ? r[x,y] = sqrt(1 - sum_err2[x,y]/sum_delta2[y]) : r[x,y] = 0
 
+                D[x,y] = b[x,y]*b[x,y] - 4.0*a[x,y]*c[x,y]
+
+                # quadratic roots (x-intercepts)
+                if (c[x,y]) {
+                    rx0[x,y] = (-1.0*b[x,y] - sqrt(D[x,y])/(2.0*c[x,y]))
+                    rx1[x,y] = (-1.0*b[x,y] + sqrt(D[x,y])/(2.0*c[x,y]))
+                }
+
                 # vertex of parabola
                 if (c[x,y]) {
-                    xv[x,y] = -1.0*b[x,y]/(2.0*c[x,y])
-                    yv[x,y] = -1.0*(b[x,y]*b[x,y])/(4.0*c[x,y]) + a[x,y]
+                    xv[x,y] = -0.5*b[x,y]/c[x,y]
+                    yv[x,y] = -0.25*(b[x,y]*b[x,y])/c[x,y] + a[x,y]
                 }
 
-                # roots (x-intercept)
+                # focus of parabola
                 if (c[x,y]) {
-                    rx0[x,y] = (-1.0*b[x,y] - sqrt(b[x,y]*b[x,y] - 4.0*a[x,y]*c[x,y]))/(2.0*c[x,y])
-                    rx1[x,y] = (-1.0*b[x,y] + sqrt(b[x,y]*b[x,y] - 4.0*a[x,y]*c[x,y]))/(2.0*c[x,y])
+                    xf[x,y] = -0.5*b[x,y]/c[x,y]
+                    yf[x,y] = (4.0*a[x,y]*c[x,y] - b[x,y]*b[x,y] + 1.0)/(4.0*c[x,y])
                 }
             }
         }
@@ -112,8 +121,10 @@ END {
             if (x != y && r[x,y]) {
                 printf(OFMT OFS "(%s)" OFS " = (" OFMT ")(%s)^2" OFS " + (" OFMT ")(%s)" OFS " + (" OFMT ")",
                     10.0*log(r[x,y]*r[x,y])/log(10), header[y], c[x,y], header[x], b[x,y], header[x], a[x,y])
-                printf("\t[" OFMT "," OFMT "][" OFMT "," OFMT "][" OFMT "," OFMT "]" OFS" [" OFMT "," OFMT "]" ORS,
+                printf("  [" OFMT "," OFMT "][" OFMT "," OFMT "][" OFMT "," OFMT "]" OFS" [" OFMT "," OFMT "]",
                     rx0[x,y], 0, rx1[x,y], 0, 0, a[x,y], xv[x,y], yv[x,y])
+                printf("  [" OFMT "," OFMT "]", xf[x,y], yf[x,y])
+                printf("  [" OFMT "]" ORS, sqrt((yf[x,y] - yv[x,y])*(yf[x,y] - yv[x,y])))
             }
         }
     }