src/gromacs/gmxlib/xvgr.c

   1 /*  -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROningen Mixture of Alchemy and Childrens' Stories
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <string.h>
  41 #include <ctype.h>
  42 #include <time.h>
  43
  44 #ifdef HAVE_SYS_TIME_H
  45 #include <sys/time.h>
  46 #endif
  47
  48 #include "sysstuff.h"
  49 #include "string2.h"
  50 #include "futil.h"
  51 #include "statutil.h"
  52 #include "copyrite.h"
  53 #include "smalloc.h"
  54 #include "xvgr.h"
  55 #include "viewit.h"
  56 #include "vec.h"
  57 #include "gmxfio.h"
  58
  59 /* Portable version of ctime_r implemented in src/gmxlib/string2.c, but we do not want it declared in public installed headers */
  60 char *
  61 gmx_ctime_r(const time_t *clock, char *buf, int n);
  62
  63
  64 gmx_bool output_env_get_print_xvgr_codes(const output_env_t oenv)
  65 {
  66     int xvg_format;
  67
  68     xvg_format = output_env_get_xvg_format(oenv);
  69
  70     return (xvg_format == exvgXMGRACE || xvg_format == exvgXMGR);
  71 }
  72
  73 static char *xvgrstr(const char *gmx, const output_env_t oenv,
  74                      char *buf, int buflen)
  75 {
  76     /* Supported greek letter names and corresponding xmgrace/xmgr symbols */
  77     const char *sym[]  = { "beta", "chi", "delta", "eta", "lambda", "mu", "omega", "phi", "psi", "rho", "theta", NULL };
  78     const char  symc[] = { 'b',    'c',   'd',     'h',   'l',      'm',  'w',     'f',   'y',   'r',   'q',     '\0' };
  79     int         xvgf;
  80     gmx_bool    bXVGR;
  81     int         g, b, i;
  82     char        c;
  83
  84     xvgf  = output_env_get_xvg_format(oenv);
  85     bXVGR = (xvgf == exvgXMGRACE || xvgf == exvgXMGR);
  86
  87     g = 0;
  88     b = 0;
  89     while (gmx[g] != '\0')
  90     {
  91         /* Check with the largest string we have ("lambda"), add one for \0 */
  92         if (b + 6 + 1 >= buflen)
  93         {
  94             gmx_fatal(FARGS, "Output buffer length in xvgstr (%d) too small to process xvg input string '%s'", buflen, gmx);
  95         }
  96         if (gmx[g] == '\\')
  97         {
  98             g++;
  99             if (gmx[g] == 's')
 100             {
 101                 /* Subscript */
 102                 if (bXVGR)
 103                 {
 104                     buf[b++] = '\\';
 105                     buf[b++] = 's';
 106                 }
 107                 else
 108                 {
 109                     buf[b++] = '_';
 110                 }
 111                 g++;
 112             }
 113             else if (gmx[g] == 'S')
 114             {
 115                 /* Superscript */
 116                 if (bXVGR)
 117                 {
 118                     buf[b++] = '\\';
 119                     buf[b++] = 'S';
 120                 }
 121                 else
 122                 {
 123                     buf[b++] = '^';
 124                 }
 125                 g++;
 126             }
 127             else if (gmx[g] == 'N')
 128             {
 129                 /* End sub/superscript */
 130                 if (bXVGR)
 131                 {
 132                     buf[b++] = '\\';
 133                     buf[b++] = 'N';
 134                 }
 135                 else
 136                 {
 137                     if (gmx[g+1] != ' ')
 138                     {
 139                         buf[b++] = ' ';
 140                     }
 141                 }
 142                 g++;
 143             }
 144             else if (gmx[g] == '4')
 145             {
 146                 /* Backward compatibility for xmgr normal font "\4" */
 147                 switch (xvgf)
 148                 {
 149                     case exvgXMGRACE:
 150                         sprintf(buf+b, "%s", "\\f{}");
 151                         break;
 152                     case exvgXMGR:
 153                         sprintf(buf+b, "%s", "\\4");
 154                         break;
 155                     default:
 156                         buf[b] = '\0';
 157                         break;
 158                 }
 159                 g++;
 160                 b = strlen(buf);
 161             }
 162             else if (gmx[g] == '8')
 163             {
 164                 /* Backward compatibility for xmgr symbol font "\8" */
 165                 switch (xvgf)
 166                 {
 167                     case exvgXMGRACE:
 168                         sprintf(buf+b, "%s", "\\x");
 169                         break;
 170                     case exvgXMGR:
 171                         sprintf(buf+b, "%s", "\\8");
 172                         break;
 173                     default:
 174                         buf[b] = '\0';
 175                         break;
 176                 }
 177                 g++;
 178                 b = strlen(buf);
 179             }
 180             else
 181             {
 182                 /* Check for special symbol */
 183                 i = 0;
 184                 while (sym[i] != NULL &&
 185                        gmx_strncasecmp(sym[i], gmx+g, strlen(sym[i])) != 0)
 186                 {
 187                     i++;
 188                 }
 189                 if (sym[i] != NULL)
 190                 {
 191                     c = symc[i];
 192                     if (isupper(gmx[g]))
 193                     {
 194                         c = toupper(c);
 195                     }
 196                     switch (xvgf)
 197                     {
 198                         case exvgXMGRACE:
 199                             sprintf(buf+b, "%s%c%s", "\\x", c, "\\f{}");
 200                             break;
 201                         case exvgXMGR:
 202                             sprintf(buf+b, "%s%c%s", "\\8", c, "\\4");
 203                             break;
 204                         default:
 205                             strncat(buf+b, gmx+g, strlen(sym[i]));
 206                             b += strlen(sym[i]);
 207                             if (gmx[g+strlen(sym[i])] != ' ')
 208                             {
 209                                 buf[b++] = ' ';
 210                             }
 211                             buf[b] = '\0';
 212                             break;
 213                     }
 214                     g += strlen(sym[i]);
 215                     b  = strlen(buf);
 216                 }
 217                 else
 218                 {
 219                     /* Unknown escape sequence, this should not happen.
 220                      * We do not generate a fatal error, since that might
 221                      * stop otherwise functioning code from working.
 222                      * Copy the backslash to the output and continue processing.
 223                      */
 224                     buf[b++] = '\\';
 225                 }
 226             }
 227         }
 228         else
 229         {
 230             buf[b++] = gmx[g++];
 231         }
 232     }
 233
 234     buf[b++] = '\0';
 235
 236     return buf;
 237 }
 238
 239 void xvgr_header(FILE *fp, const char *title, const char *xaxis,
 240                  const char *yaxis, int exvg_graph_type,
 241                  const output_env_t oenv)
 242 {
 243     char   pukestr[100], buf[STRLEN];
 244     time_t t;
 245
 246     if (output_env_get_print_xvgr_codes(oenv))
 247     {
 248         time(&t);
 249         gmx_ctime_r(&t, buf, STRLEN);
 250         fprintf(fp, "# This file was created %s", buf);
 251         fprintf(fp, "# by the following command:\n# %s\n#\n", command_line());
 252         fprintf(fp, "# %s is part of G R O M A C S:\n#\n", ShortProgram());
 253         bromacs(pukestr, 99);
 254         fprintf(fp, "# %s\n#\n", pukestr);
 255         fprintf(fp, "@    title \"%s\"\n", xvgrstr(title, oenv, buf, STRLEN));
 256         fprintf(fp, "@    xaxis  label \"%s\"\n",
 257                 xvgrstr(xaxis, oenv, buf, STRLEN));
 258         fprintf(fp, "@    yaxis  label \"%s\"\n",
 259                 xvgrstr(yaxis, oenv, buf, STRLEN));
 260         switch (exvg_graph_type)
 261         {
 262             case exvggtXNY:
 263                 if (output_env_get_xvg_format(oenv) == exvgXMGR)
 264                 {
 265                     fprintf(fp, "@TYPE nxy\n");
 266                 }
 267                 else
 268                 {
 269                     fprintf(fp, "@TYPE xy\n");
 270                 }
 271                 break;
 272             case exvggtXYDY:
 273                 fprintf(fp, "@TYPE xydy\n");
 274                 break;
 275             case exvggtXYDYDY:
 276                 fprintf(fp, "@TYPE xydydy\n");
 277                 break;
 278         }
 279     }
 280 }
 281
 282 FILE *xvgropen_type(const char *fn, const char *title, const char *xaxis,
 283                     const char *yaxis, int exvg_graph_type,
 284                     const output_env_t oenv)
 285 {
 286     FILE  *fp;
 287     time_t t;
 288
 289     fp = gmx_fio_fopen(fn, "w");
 290
 291     xvgr_header(fp, title, xaxis, yaxis, exvg_graph_type, oenv);
 292
 293     return fp;
 294 }
 295
 296 FILE *xvgropen(const char *fn, const char *title, const char *xaxis,
 297                const char *yaxis, const output_env_t oenv)
 298 {
 299     return xvgropen_type(fn, title, xaxis, yaxis, exvggtXNY, oenv);
 300 }
 301
 302 void
 303 xvgrclose(FILE *fp)
 304 {
 305     gmx_fio_fclose(fp);
 306 }
 307
 308 void xvgr_subtitle(FILE *out, const char *subtitle, const output_env_t oenv)
 309 {
 310     char buf[STRLEN];
 311
 312     if (output_env_get_print_xvgr_codes(oenv))
 313     {
 314         fprintf(out, "@ subtitle \"%s\"\n", xvgrstr(subtitle, oenv, buf, STRLEN));
 315     }
 316 }
 317
 318 void xvgr_view(FILE *out, real xmin, real ymin, real xmax, real ymax,
 319                const output_env_t oenv)
 320 {
 321     if (output_env_get_print_xvgr_codes(oenv))
 322     {
 323         fprintf(out, "@ view %g, %g, %g, %g\n", xmin, ymin, xmax, ymax);
 324     }
 325 }
 326
 327 void xvgr_world(FILE *out, real xmin, real ymin, real xmax, real ymax,
 328                 const output_env_t oenv)
 329 {
 330     if (output_env_get_print_xvgr_codes(oenv))
 331     {
 332         fprintf(out, "@ world xmin %g\n"
 333                 "@ world ymin %g\n"
 334                 "@ world xmax %g\n"
 335                 "@ world ymax %g\n", xmin, ymin, xmax, ymax);
 336     }
 337 }
 338
 339 void xvgr_legend(FILE *out, int nsets, const char** setname,
 340                  const output_env_t oenv)
 341 {
 342     int  i;
 343     char buf[STRLEN];
 344
 345     if (output_env_get_print_xvgr_codes(oenv))
 346     {
 347         xvgr_view(out, 0.15, 0.15, 0.75, 0.85, oenv);
 348         fprintf(out, "@ legend on\n");
 349         fprintf(out, "@ legend box on\n");
 350         fprintf(out, "@ legend loctype view\n");
 351         fprintf(out, "@ legend %g, %g\n", 0.78, 0.8);
 352         fprintf(out, "@ legend length %d\n", 2);
 353         for (i = 0; (i < nsets); i++)
 354         {
 355             if (setname[i])
 356             {
 357                 if (output_env_get_xvg_format(oenv) == exvgXMGR)
 358                 {
 359                     fprintf(out, "@ legend string %d \"%s\"\n",
 360                             i, xvgrstr(setname[i], oenv, buf, STRLEN));
 361                 }
 362                 else
 363                 {
 364                     fprintf(out, "@ s%d legend \"%s\"\n",
 365                             i, xvgrstr(setname[i], oenv, buf, STRLEN));
 366                 }
 367             }
 368         }
 369     }
 370 }
 371
 372 void xvgr_new_dataset(FILE *out, int nr_first, int nsets,
 373                       const char **setname,
 374                       const output_env_t oenv)
 375 {
 376     int  i;
 377     char buf[STRLEN];
 378
 379     if (output_env_get_print_xvgr_codes(oenv))
 380     {
 381         fprintf(out, "@\n");
 382         for (i = 0; (i < nsets); i++)
 383         {
 384             if (setname[i])
 385             {
 386                 if (output_env_get_xvg_format(oenv) == exvgXMGR)
 387                 {
 388                     fprintf(out, "@ legend string %d \"%s\"\n",
 389                             i+nr_first, xvgrstr(setname[i], oenv, buf, STRLEN));
 390                 }
 391                 else
 392                 {
 393                     fprintf(out, "@ s%d legend \"%s\"\n",
 394                             i+nr_first, xvgrstr(setname[i], oenv, buf, STRLEN));
 395                 }
 396             }
 397         }
 398     }
 399     else
 400     {
 401         fprintf(out, "\n");
 402     }
 403 }
 404
 405 void xvgr_line_props(FILE *out, int NrSet, int LineStyle, int LineColor,
 406                      const output_env_t oenv)
 407 {
 408     if (output_env_get_print_xvgr_codes(oenv))
 409     {
 410         fprintf(out, "@    with g0\n");
 411         fprintf(out, "@    s%d linestyle %d\n", NrSet, LineStyle);
 412         fprintf(out, "@    s%d color %d\n", NrSet, LineColor);
 413     }
 414 }
 415
 416 static const char *LocTypeStr[] = { "view", "world" };
 417 static const char *BoxFillStr[] = { "none", "color", "pattern" };
 418
 419 void xvgr_box(FILE *out,
 420               int LocType,
 421               real xmin, real ymin, real xmax, real ymax,
 422               int LineStyle, int LineWidth, int LineColor,
 423               int BoxFill, int BoxColor, int BoxPattern, const output_env_t oenv)
 424 {
 425     if (output_env_get_print_xvgr_codes(oenv))
 426     {
 427         fprintf(out, "@with box\n");
 428         fprintf(out, "@    box on\n");
 429         fprintf(out, "@    box loctype %s\n", LocTypeStr[LocType]);
 430         fprintf(out, "@    box %g, %g, %g, %g\n", xmin, ymin, xmax, ymax);
 431         fprintf(out, "@    box linestyle %d\n", LineStyle);
 432         fprintf(out, "@    box linewidth %d\n", LineWidth);
 433         fprintf(out, "@    box color %d\n", LineColor);
 434         fprintf(out, "@    box fill %s\n", BoxFillStr[BoxFill]);
 435         fprintf(out, "@    box fill color %d\n", BoxColor);
 436         fprintf(out, "@    box fill pattern %d\n", BoxPattern);
 437         fprintf(out, "@box def\n");
 438     }
 439 }
 440
 441 /* reads a line into ptr, adjusting len and renewing ptr if neccesary */
 442 static char *fgets3(FILE *fp, char **ptr, int *len, int maxlen)
 443 {
 444     char *p;
 445     int   len_remaining = *len; /* remaining amount of allocated bytes in buf */
 446     int   curp          = 0;    /* current position in buf to read into */
 447
 448     do
 449     {
 450         if (len_remaining < 2)
 451         {
 452             if (*len + STRLEN < maxlen)
 453             {
 454                 /* This line is longer than len characters, let's increase len! */
 455                 *len          += STRLEN;
 456                 len_remaining += STRLEN;
 457                 srenew(*ptr, *len);
 458             }
 459             else
 460             {
 461                 /*something is wrong, we'll just keep reading and return NULL*/
 462                 len_remaining = STRLEN;
 463                 curp          = 0;
 464             }
 465         }
 466         if (fgets(*ptr + curp, len_remaining, fp) == NULL)
 467         {
 468             /* if last line, skip */
 469             return NULL;
 470         }
 471         curp         += len_remaining-1; /* overwrite the nul char in next iteration */
 472         len_remaining = 1;
 473     }
 474     while ((strchr(*ptr, '\n') == NULL) && (!feof(fp)));
 475
 476     if (*len + STRLEN >= maxlen)
 477     {
 478         return NULL; /* this line was too long */
 479     }
 480
 481     if (feof(fp))
 482     {
 483         /* We reached EOF before '\n', skip this last line. */
 484         return NULL;
 485     }
 486     {
 487         /* now remove newline */
 488         int slen = strlen(*ptr);
 489         if ((*ptr)[slen-1] == '\n')
 490         {
 491             (*ptr)[slen-1] = '\0';
 492         }
 493     }
 494
 495     return *ptr;
 496 }
 497
 498 static int wordcount(char *ptr)
 499 {
 500     int i, n, is[2];
 501     int cur = 0;
 502 #define prev (1-cur)
 503
 504     if (strlen(ptr) == 0)
 505     {
 506         return 0;
 507     }
 508     /* fprintf(stderr,"ptr='%s'\n",ptr); */
 509     n = 1;
 510     for (i = 0; (ptr[i] != '\0'); i++)
 511     {
 512         is[cur] = isspace(ptr[i]);
 513         if ((i > 0)  && (is[cur] && !is[prev]))
 514         {
 515             n++;
 516         }
 517         cur = prev;
 518     }
 519     return n;
 520 }
 521
 522 static char *read_xvgr_string(const char *line)
 523 {
 524     const char *ptr0, *ptr1;
 525     char       *str;
 526
 527     ptr0 = strchr(line, '"');
 528     if (ptr0 != NULL)
 529     {
 530         ptr0++;
 531         ptr1 = strchr(ptr0, '"');
 532         if (ptr1 != NULL)
 533         {
 534             str            = strdup(ptr0);
 535             str[ptr1-ptr0] = '\0';
 536         }
 537         else
 538         {
 539             str = strdup("");
 540         }
 541     }
 542     else
 543     {
 544         str = strdup("");
 545     }
 546
 547     return str;
 548 }
 549
 550 int read_xvg_legend(const char *fn, double ***y, int *ny,
 551                     char **subtitle, char ***legend)
 552 {
 553     FILE    *fp;
 554     char    *ptr, *ptr0, *ptr1;
 555     char    *base = NULL;
 556     char    *fmt  = NULL;
 557     int      k, line = 0, nny, nx, maxx, rval, legend_nalloc, set, nchar;
 558     double   lf;
 559     double **yy = NULL;
 560     char    *tmpbuf;
 561     int      len = STRLEN;
 562     *ny  = 0;
 563     nny  = 0;
 564     nx   = 0;
 565     maxx = 0;
 566     fp   = gmx_fio_fopen(fn, "r");
 567
 568     snew(tmpbuf, len);
 569     if (subtitle != NULL)
 570     {
 571         *subtitle = NULL;
 572     }
 573     legend_nalloc = 0;
 574     if (legend != NULL)
 575     {
 576         *legend = NULL;
 577     }
 578
 579     while ((ptr = fgets3(fp, &tmpbuf, &len, 10*STRLEN)) != NULL && ptr[0] != '&')
 580     {
 581         line++;
 582         trim(ptr);
 583         if (ptr[0] == '@')
 584         {
 585             if (legend != NULL)
 586             {
 587                 ptr++;
 588                 trim(ptr);
 589                 set = -1;
 590                 if (strncmp(ptr, "subtitle", 8) == 0)
 591                 {
 592                     ptr += 8;
 593                     if (subtitle != NULL)
 594                     {
 595                         *subtitle = read_xvgr_string(ptr);
 596                     }
 597                 }
 598                 else if (strncmp(ptr, "legend string", 13) == 0)
 599                 {
 600                     ptr += 13;
 601                     sscanf(ptr, "%d%n", &set, &nchar);
 602                     ptr += nchar;
 603                 }
 604                 else if (ptr[0] == 's')
 605                 {
 606                     ptr++;
 607                     sscanf(ptr, "%d%n", &set, &nchar);
 608                     ptr += nchar;
 609                     trim(ptr);
 610                     if (strncmp(ptr, "legend", 6) == 0)
 611                     {
 612                         ptr += 6;
 613                     }
 614                     else
 615                     {
 616                         set = -1;
 617                     }
 618                 }
 619                 if (set >= 0)
 620                 {
 621                     if (set >= legend_nalloc)
 622                     {
 623                         legend_nalloc = set + 1;
 624                         srenew(*legend, legend_nalloc);
 625                         (*legend)[set] = read_xvgr_string(ptr);
 626                     }
 627                 }
 628             }
 629         }
 630         else if (ptr[0] != '#')
 631         {
 632             if (nny == 0)
 633             {
 634                 (*ny) = nny = wordcount(ptr);
 635                 /* fprintf(stderr,"There are %d columns in your file\n",nny);*/
 636                 if (nny == 0)
 637                 {
 638                     return 0;
 639                 }
 640                 snew(yy, nny);
 641                 snew(fmt, 3*nny+1);
 642                 snew(base, 3*nny+1);
 643             }
 644             /* Allocate column space */
 645             if (nx >= maxx)
 646             {
 647                 maxx += 1024;
 648                 for (k = 0; (k < nny); k++)
 649                 {
 650                     srenew(yy[k], maxx);
 651                 }
 652             }
 653             /* Initiate format string */
 654             fmt[0]  = '\0';
 655             base[0] = '\0';
 656
 657             /* fprintf(stderr,"ptr='%s'\n",ptr);*/
 658             for (k = 0; (k < nny); k++)
 659             {
 660                 strcpy(fmt, base);
 661                 strcat(fmt, "%lf");
 662                 rval = sscanf(ptr, fmt, &lf);
 663                 /* fprintf(stderr,"rval = %d\n",rval);*/
 664                 if ((rval == EOF) || (rval == 0))
 665                 {
 666                     break;
 667                 }
 668                 yy[k][nx] = lf;
 669                 srenew(fmt, 3*(nny+1)+1);
 670                 srenew(base, 3*nny+1);
 671                 strcat(base, "%*s");
 672             }
 673             if (k != nny)
 674             {
 675                 fprintf(stderr, "Only %d columns on line %d in file %s\n",
 676                         k, line, fn);
 677                 for (; (k < nny); k++)
 678                 {
 679                     yy[k][nx] = 0.0;
 680                 }
 681             }
 682             nx++;
 683         }
 684     }
 685     gmx_fio_fclose(fp);
 686
 687     *y = yy;
 688     sfree(tmpbuf);
 689
 690     if (legend_nalloc > 0)
 691     {
 692         if (*ny - 1 > legend_nalloc)
 693         {
 694             srenew(*legend, *ny-1);
 695             for (set = legend_nalloc; set < *ny-1; set++)
 696             {
 697                 (*legend)[set] = NULL;
 698             }
 699         }
 700     }
 701
 702     return nx;
 703 }
 704
 705 int read_xvg(const char *fn, double ***y, int *ny)
 706 {
 707     return read_xvg_legend(fn, y, ny, NULL, NULL);
 708 }
 709
 710 void write_xvg(const char *fn, const char *title, int nx, int ny, real **y,
 711                const char **leg, const output_env_t oenv)
 712 {
 713     FILE *fp;
 714     int   i, j;
 715
 716     fp = xvgropen(fn, title, "X", "Y", oenv);
 717     if (leg)
 718     {
 719         xvgr_legend(fp, ny-1, leg, oenv);
 720     }
 721     for (i = 0; (i < nx); i++)
 722     {
 723         for (j = 0; (j < ny); j++)
 724         {
 725             fprintf(fp, "  %12.5e", y[j][i]);
 726         }
 727         fprintf(fp, "\n");
 728     }
 729     xvgrclose(fp);
 730 }
 731
 732 real **read_xvg_time(const char *fn,
 733                      gmx_bool bHaveT, gmx_bool bTB, real tb, gmx_bool bTE, real te,
 734                      int nsets_in, int *nset, int *nval, real *dt, real **t)
 735 {
 736     FILE      *fp;
 737 #define MAXLINELEN 16384
 738     char       line0[MAXLINELEN];
 739     char      *line;
 740     int        t_nalloc, *val_nalloc, a, narg, n, sin, set, nchar;
 741     double     dbl, tend = 0;
 742     gmx_bool   bEndOfSet, bTimeInRange, bFirstLine = TRUE;
 743     real     **val;
 744
 745     t_nalloc   = 0;
 746     *t         = NULL;
 747     val        = NULL;
 748     val_nalloc = NULL;
 749     *nset      = 0;
 750     *dt        = 0;
 751     fp         = gmx_fio_fopen(fn, "r");
 752     for (sin = 0; sin < nsets_in; sin++)
 753     {
 754         if (nsets_in == 1)
 755         {
 756             narg = 0;
 757         }
 758         else
 759         {
 760             narg = bHaveT ? 2 : 1;
 761         }
 762         n         = 0;
 763         bEndOfSet = FALSE;
 764         while (!bEndOfSet && fgets(line0, MAXLINELEN, fp))
 765         {
 766             line = line0;
 767             /* Remove whitespace */
 768             while (line[0] == ' ' || line[0] == '\t')
 769             {
 770                 line++;
 771             }
 772             bEndOfSet = (line[0] == '&');
 773             if (line[0] != '#' && line[0] != '@' && line[0] != '\n' && !bEndOfSet)
 774             {
 775                 if (bFirstLine && bHaveT)
 776                 {
 777                     /* Check the first line that should contain data */
 778                     a = sscanf(line, "%lf%lf", &dbl, &dbl);
 779                     if (a == 0)
 780                     {
 781                         gmx_fatal(FARGS, "Expected a number in %s on line:\n%s", fn, line0);
 782                     }
 783                     else if (a == 1)
 784                     {
 785                         fprintf(stderr, "Found only 1 number on line, "
 786                                 "assuming no time is present.\n");
 787                         bHaveT = FALSE;
 788                         if (nsets_in > 1)
 789                         {
 790                             narg = 1;
 791                         }
 792                     }
 793                 }
 794
 795                 a            = 0;
 796                 bTimeInRange = TRUE;
 797                 while ((a < narg || (nsets_in == 1 && n == 0)) &&
 798                        sscanf(line, "%lf%n", &dbl, &nchar) == 1 && bTimeInRange)
 799                 {
 800                     /* Use set=-1 as the time "set" */
 801                     if (sin)
 802                     {
 803                         if (!bHaveT || (a > 0))
 804                         {
 805                             set = sin;
 806                         }
 807                         else
 808                         {
 809                             set = -1;
 810                         }
 811                     }
 812                     else
 813                     {
 814                         if (!bHaveT)
 815                         {
 816                             set = a;
 817                         }
 818                         else
 819                         {
 820                             set = a-1;
 821                         }
 822                     }
 823                     if (set == -1 && ((bTB && dbl < tb) || (bTE && dbl > te)))
 824                     {
 825                         bTimeInRange = FALSE;
 826                     }
 827
 828                     if (bTimeInRange)
 829                     {
 830                         if (n == 0)
 831                         {
 832                             if (nsets_in == 1)
 833                             {
 834                                 narg++;
 835                             }
 836                             if (set >= 0)
 837                             {
 838                                 *nset = set+1;
 839                                 srenew(val, *nset);
 840                                 srenew(val_nalloc, *nset);
 841                                 val_nalloc[set] = 0;
 842                                 val[set]        = NULL;
 843                             }
 844                         }
 845                         if (set == -1)
 846                         {
 847                             if (sin == 0)
 848                             {
 849                                 if (n >= t_nalloc)
 850                                 {
 851                                     t_nalloc = over_alloc_small(n);
 852                                     srenew(*t, t_nalloc);
 853                                 }
 854                                 (*t)[n] = dbl;
 855                             }
 856                             /* else we should check the time of the next sets with set 0 */
 857                         }
 858                         else
 859                         {
 860                             if (n >= val_nalloc[set])
 861                             {
 862                                 val_nalloc[set] = over_alloc_small(n);
 863                                 srenew(val[set], val_nalloc[set]);
 864                             }
 865                             val[set][n] = (real)dbl;
 866                         }
 867                     }
 868                     a++;
 869                     line += nchar;
 870                 }
 871                 if (line0[strlen(line0)-1] != '\n')
 872                 {
 873                     fprintf(stderr, "File %s does not end with a newline, ignoring the last line\n", fn);
 874                 }
 875                 else if (bTimeInRange)
 876                 {
 877                     if (a == 0)
 878                     {
 879                         fprintf(stderr, "Ignoring invalid line in %s:\n%s", fn, line0);
 880                     }
 881                     else
 882                     {
 883                         if (a != narg)
 884                         {
 885                             fprintf(stderr, "Invalid line in %s:\n%s"
 886                                     "Using zeros for the last %d sets\n",
 887                                     fn, line0, narg-a);
 888                         }
 889                         n++;
 890                     }
 891                 }
 892                 if (a > 0)
 893                 {
 894                     bFirstLine = FALSE;
 895                 }
 896             }
 897         }
 898         if (sin == 0)
 899         {
 900             *nval = n;
 901             if (!bHaveT)
 902             {
 903                 snew(*t, n);
 904                 for (a = 0; a < n; a++)
 905                 {
 906                     (*t)[a] = a;
 907                 }
 908             }
 909             if (n > 1)
 910             {
 911                 *dt = (real)((*t)[n-1]-(*t)[0])/(n-1.0);
 912             }
 913             else
 914             {
 915                 *dt = 1;
 916             }
 917         }
 918         else
 919         {
 920             if (n < *nval)
 921             {
 922                 fprintf(stderr, "Set %d is shorter (%d) than the previous set (%d)\n",
 923                         sin+1, n, *nval);
 924                 *nval = n;
 925                 fprintf(stderr, "Will use only the first %d points of every set\n",
 926                         *nval);
 927             }
 928         }
 929     }
 930     gmx_fio_fclose(fp);
 931
 932     sfree(val_nalloc);
 933
 934     return val;
 935 }