From: W. Trevor King <wking@tremily.us>
Date: Tue, 12 Feb 2013 21:54:47 +0000 (-0500)
Subject: error_analysis.c: Implement error analysis
X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=d8c1238ec3ded1886970fa6ecacc1e5d823394b5;p=assignment-template.git

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 c242304..924af91 100644
--- 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;
 }