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