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