src/gromacs/fileio/xvgr.cpp

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