error_analysis.c

   1 /*
   2 Calculate the maximum relative energy error in a trajectory.
   3
   4 Copyright (C) 2013  W. Trevor King
   5
   6 This program is free software: you can redistribute it and/or modify
   7 it under the terms of the GNU General Public License as published by
   8 the Free Software Foundation, either version 3 of the License, or
   9 (at your option) any later version.
  10
  11 This program is distributed in the hope that it will be useful,
  12 but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 GNU General Public License for more details.
  15
  16 You should have received a copy of the GNU General Public License
  17 along with this program.  If not, see <http://www.gnu.org/licenses/>.
  18 */
  19
  20 #include <stdio.h>  /* for *printf(), getc(), fscanf(), EOF, stdin, stderr */
  21 #include <math.h>   /* for fabs() */
  22
  23 /* Read a double column from the current line in a tab-delimited
  24  * stream.
  25  */
  26 int read_column(FILE *stream, int column, size_t line, double *value)
  27 {
  28         int c, col = 0;
  29
  30         while ((c = getc(stream)) != EOF) {
  31                 if (c == '\t') {  /* found a new column */
  32                         col++;
  33                         if (col == column) {  /* read our value */
  34                                 c = fscanf(stream, "%lf", value);
  35                                 if (c != 1) {
  36                                         fprintf(stderr, "could not parse column %d in line %ld\n",
  37                                                 column, line);
  38                                         return -1;
  39                                 }
  40                                 while ((c = getc(stream)) != EOF)
  41                                         if (c == '\n')
  42                                                 return 0;  /* advanced to next line */
  43                         }
  44                 } else if (c == '\n') {
  45                         fprintf(stderr, "only %d columns in line %ld\n", column, line);
  46                         return -1;
  47                 }
  48         }
  49         return -1;  /* reached the end of the file */
  50 }
  51
  52 int main(int argc, char **argv) {
  53         int line, column = 5;
  54         double val_i, val, error, max_error = 0;
  55
  56         for (line=0; read_column(stdin, column, line, &val) == 0; line++) {
  57                 if (line == 0) {
  58                         val_i = val;
  59                         if (val_i == 0) {
  60                                 fprintf(stderr,
  61                                     "cannot calculate relative error when initial value is zero\n");
  62                                 return -1;
  63                         }
  64                 } else {
  65                         error = fabs((val - val_i)/(val_i));
  66                         if (error > max_error)
  67                                 max_error = error;
  68                 }
  69         }
  70         printf("%g\n", max_error);
  71         return 0;
  72 }