error_analysis.c: Implement error analysis

Currently it is hard-coded to read column five, but adding a command
line switch to change columns wouldn't be too difficult.

diff --git a/error_analysis.c b/error_analysis.c
index c2423049f8ba3e94dc3161704eb9fce828895f8c..924af912922bc41ab442c7df9359b9fde1ce57cd 100644 (file)
--- a/error_analysis.c
+++ b/error_analysis.c
@@ -1,7 +1,7 @@
 /*
-A simple "hello world" example in C.
+Calculate the maximum relative energy error in a trajectory.
 
-Copyright (C) 2012  W. Trevor King
+Copyright (C) 2013  W. Trevor King
 
 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -17,9 +17,56 @@ You should have received a copy of the GNU General Public License
 along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include <stdio.h>
+#include <stdio.h>  /* for *printf(), getc(), fscanf(), EOF, stdin, stderr */
+#include <math.h>   /* for fabs() */
 
-int main() {
-       printf("hello, world\n");
+/* Read a double column from the current line in a tab-delimited
+ * stream.
+ */
+int read_column(FILE *stream, int column, size_t line, double *value)
+{
+       int c, col = 0;
+
+       while ((c = getc(stream)) != EOF) {
+               if (c == '\t') {  /* found a new column */
+                       col++;
+                       if (col == column) {  /* read our value */
+                               c = fscanf(stream, "%lf", value);
+                               if (c != 1) {
+                                       fprintf(stderr, "could not parse column %d in line %ld\n",
+                                               column, line);
+                                       return -1;
+                               }
+                               while ((c = getc(stream)) != EOF)
+                                       if (c == '\n')
+                                               return 0;  /* advanced to next line */
+                       }
+               } else if (c == '\n') {
+                       fprintf(stderr, "only %d columns in line %ld\n", column, line);
+                       return -1;
+               }
+       }
+       return -1;  /* reached the end of the file */
+}
+
+int main(int argc, char **argv) {
+       int line, column = 5;
+       double val_i, val, error, max_error = 0;
+
+       for (line=0; read_column(stdin, column, line, &val) == 0; line++) {
+               if (line == 0) {
+                       val_i = val;
+                       if (val_i == 0) {
+                               fprintf(stderr,
+                                   "cannot calculate relative error when initial value is zero\n");
+                               return -1;
+                       }
+               } else {
+                       error = fabs((val - val_i)/(val_i));
+                       if (error > max_error)
+                               max_error = error;
+               }
+       }
+       printf("%g\n", max_error);
        return 0;
 }