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,2015,2016, 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 <cassert>
  42 #include <cctype>
  43 #include <cstring>
  44
  45 #include <string>
  46
  47 #include "gromacs/fileio/gmxfio.h"
  48 #include "gromacs/fileio/oenv.h"
  49 #include "gromacs/math/vec.h"
  50 #include "gromacs/utility/binaryinformation.h"
  51 #include "gromacs/utility/coolstuff.h"
  52 #include "gromacs/utility/cstringutil.h"
  53 #include "gromacs/utility/exceptions.h"
  54 #include "gromacs/utility/fatalerror.h"
  55 #include "gromacs/utility/futil.h"
  56 #include "gromacs/utility/smalloc.h"
  57 #include "gromacs/utility/sysinfo.h"
  58
  59 gmx_bool output_env_get_print_xvgr_codes(const gmx_output_env_t *oenv)
  60 {
  61     int xvg_format;
  62
  63     xvg_format = output_env_get_xvg_format(oenv);
  64
  65     return (xvg_format == exvgXMGRACE || xvg_format == exvgXMGR);
  66 }
  67
  68 static char *xvgrstr(const char *gmx, const gmx_output_env_t *oenv,
  69                      char *buf, int buflen)
  70 {
  71     /* Supported greek letter names and corresponding xmgrace/xmgr symbols */
  72     const char *sym[]  = { "beta", "chi", "delta", "eta", "lambda", "mu", "omega", "phi", "psi", "rho", "theta", NULL };
  73     const char  symc[] = { 'b',    'c',   'd',     'h',   'l',      'm',  'w',     'f',   'y',   'r',   'q',     '\0' };
  74     int         xvgf;
  75     gmx_bool    bXVGR;
  76     int         g, b, i;
  77     char        c;
  78
  79     xvgf  = output_env_get_xvg_format(oenv);
  80     bXVGR = (xvgf == exvgXMGRACE || xvgf == exvgXMGR);
  81
  82     g = 0;
  83     b = 0;
  84     while (gmx[g] != '\0')
  85     {
  86         /* Check with the largest string we have ("lambda"), add one for \0 */
  87         if (b + 6 + 1 >= buflen)
  88         {
  89             gmx_fatal(FARGS, "Output buffer length in xvgstr (%d) too small to process xvg input string '%s'", buflen, gmx);
  90         }
  91         if (gmx[g] == '\\')
  92         {
  93             g++;
  94             if (gmx[g] == 's')
  95             {
  96                 /* Subscript */
  97                 if (bXVGR)
  98                 {
  99                     buf[b++] = '\\';
 100                     buf[b++] = 's';
 101                 }
 102                 else
 103                 {
 104                     buf[b++] = '_';
 105                 }
 106                 g++;
 107             }
 108             else if (gmx[g] == 'S')
 109             {
 110                 /* Superscript */
 111                 if (bXVGR)
 112                 {
 113                     buf[b++] = '\\';
 114                     buf[b++] = 'S';
 115                 }
 116                 else
 117                 {
 118                     buf[b++] = '^';
 119                 }
 120                 g++;
 121             }
 122             else if (gmx[g] == 'N')
 123             {
 124                 /* End sub/superscript */
 125                 if (bXVGR)
 126                 {
 127                     buf[b++] = '\\';
 128                     buf[b++] = 'N';
 129                 }
 130                 else
 131                 {
 132                     if (gmx[g+1] != ' ')
 133                     {
 134                         buf[b++] = ' ';
 135                     }
 136                 }
 137                 g++;
 138             }
 139             else if (gmx[g] == '4')
 140             {
 141                 /* Backward compatibility for xmgr normal font "\4" */
 142                 switch (xvgf)
 143                 {
 144                     case exvgXMGRACE:
 145                         sprintf(buf+b, "%s", "\\f{}");
 146                         break;
 147                     case exvgXMGR:
 148                         sprintf(buf+b, "%s", "\\4");
 149                         break;
 150                     default:
 151                         buf[b] = '\0';
 152                         break;
 153                 }
 154                 g++;
 155                 b = strlen(buf);
 156             }
 157             else if (gmx[g] == '8')
 158             {
 159                 /* Backward compatibility for xmgr symbol font "\8" */
 160                 switch (xvgf)
 161                 {
 162                     case exvgXMGRACE:
 163                         sprintf(buf+b, "%s", "\\x");
 164                         break;
 165                     case exvgXMGR:
 166                         sprintf(buf+b, "%s", "\\8");
 167                         break;
 168                     default:
 169                         buf[b] = '\0';
 170                         break;
 171                 }
 172                 g++;
 173                 b = std::strlen(buf);
 174             }
 175             else
 176             {
 177                 /* Check for special symbol */
 178                 i = 0;
 179                 while (sym[i] != NULL &&
 180                        gmx_strncasecmp(sym[i], gmx+g, std::strlen(sym[i])) != 0)
 181                 {
 182                     i++;
 183                 }
 184                 if (sym[i] != NULL)
 185                 {
 186                     c = symc[i];
 187                     if (std::isupper(gmx[g]))
 188                     {
 189                         c = std::toupper(c);
 190                     }
 191                     switch (xvgf)
 192                     {
 193                         case exvgXMGRACE:
 194                             sprintf(buf+b, "%s%c%s", "\\x", c, "\\f{}");
 195                             break;
 196                         case exvgXMGR:
 197                             sprintf(buf+b, "%s%c%s", "\\8", c, "\\4");
 198                             break;
 199                         default:
 200                             std::strncat(buf+b, gmx+g, std::strlen(sym[i]));
 201                             b += std::strlen(sym[i]);
 202                             if (gmx[g+std::strlen(sym[i])] != ' ')
 203                             {
 204                                 buf[b++] = ' ';
 205                             }
 206                             buf[b] = '\0';
 207                             break;
 208                     }
 209                     g += std::strlen(sym[i]);
 210                     b  = std::strlen(buf);
 211                 }
 212                 else
 213                 {
 214                     /* Unknown escape sequence, this should not happen.
 215                      * We do not generate a fatal error, since that might
 216                      * stop otherwise functioning code from working.
 217                      * Copy the backslash to the output and continue processing.
 218                      */
 219                     buf[b++] = '\\';
 220                 }
 221             }
 222         }
 223         else
 224         {
 225             buf[b++] = gmx[g++];
 226         }
 227     }
 228
 229     buf[b++] = '\0';
 230
 231     return buf;
 232 }
 233
 234 void xvgr_header(FILE *fp, const char *title, const char *xaxis,
 235                  const char *yaxis, int exvg_graph_type,
 236                  const gmx_output_env_t *oenv)
 237 {
 238     char buf[STRLEN];
 239
 240     if (output_env_get_print_xvgr_codes(oenv))
 241     {
 242         gmx_format_current_time(buf, STRLEN);
 243         fprintf(fp, "# This file was created %s", buf);
 244         try
 245         {
 246             gmx::BinaryInformationSettings settings;
 247             settings.generatedByHeader(true);
 248             settings.linePrefix("# ");
 249             gmx::printBinaryInformation(fp, output_env_get_program_context(oenv),
 250                                         settings);
 251         }
 252         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 253         fprintf(fp, "# %s is part of G R O M A C S:\n#\n",
 254                 output_env_get_program_display_name(oenv));
 255         fprintf(fp, "# %s\n#\n", gmx::bromacs().c_str());
 256         fprintf(fp, "@    title \"%s\"\n", xvgrstr(title, oenv, buf, STRLEN));
 257         fprintf(fp, "@    xaxis  label \"%s\"\n",
 258                 xvgrstr(xaxis, oenv, buf, STRLEN));
 259         fprintf(fp, "@    yaxis  label \"%s\"\n",
 260                 xvgrstr(yaxis, oenv, buf, STRLEN));
 261         switch (exvg_graph_type)
 262         {
 263             case exvggtXNY:
 264                 if (output_env_get_xvg_format(oenv) == exvgXMGR)
 265                 {
 266                     fprintf(fp, "@TYPE nxy\n");
 267                 }
 268                 else
 269                 {
 270                     fprintf(fp, "@TYPE xy\n");
 271                 }
 272                 break;
 273             case exvggtXYDY:
 274                 fprintf(fp, "@TYPE xydy\n");
 275                 break;
 276             case exvggtXYDYDY:
 277                 fprintf(fp, "@TYPE xydydy\n");
 278                 break;
 279         }
 280     }
 281 }
 282
 283 FILE *xvgropen_type(const char *fn, const char *title, const char *xaxis,
 284                     const char *yaxis, int exvg_graph_type,
 285                     const gmx_output_env_t *oenv)
 286 {
 287     FILE  *fp;
 288
 289     fp = gmx_fio_fopen(fn, "w");
 290
 291     xvgr_header(fp, title, xaxis, yaxis, exvg_graph_type, oenv);
 292
 293     return fp;
 294 }
 295
 296 FILE *xvgropen(const char *fn, const char *title, const char *xaxis,
 297                const char *yaxis, const gmx_output_env_t *oenv)
 298 {
 299     return xvgropen_type(fn, title, xaxis, yaxis, exvggtXNY, oenv);
 300 }
 301
 302 void
 303 xvgrclose(FILE *fp)
 304 {
 305     gmx_fio_fclose(fp);
 306 }
 307
 308 void xvgr_subtitle(FILE *out, const char *subtitle, const gmx_output_env_t *oenv)
 309 {
 310     char buf[STRLEN];
 311
 312     if (output_env_get_print_xvgr_codes(oenv))
 313     {
 314         fprintf(out, "@ subtitle \"%s\"\n", xvgrstr(subtitle, oenv, buf, STRLEN));
 315     }
 316 }
 317
 318 void xvgr_view(FILE *out, real xmin, real ymin, real xmax, real ymax,
 319                const gmx_output_env_t *oenv)
 320 {
 321     if (output_env_get_print_xvgr_codes(oenv))
 322     {
 323         fprintf(out, "@ view %g, %g, %g, %g\n", xmin, ymin, xmax, ymax);
 324     }
 325 }
 326
 327 void xvgr_world(FILE *out, real xmin, real ymin, real xmax, real ymax,
 328                 const gmx_output_env_t *oenv)
 329 {
 330     if (output_env_get_print_xvgr_codes(oenv))
 331     {
 332         fprintf(out, "@ world xmin %g\n"
 333                 "@ world ymin %g\n"
 334                 "@ world xmax %g\n"
 335                 "@ world ymax %g\n", xmin, ymin, xmax, ymax);
 336     }
 337 }
 338
 339 void xvgr_legend(FILE *out, int nsets, const char** setname,
 340                  const gmx_output_env_t *oenv)
 341 {
 342     int  i;
 343     char buf[STRLEN];
 344
 345     if (output_env_get_print_xvgr_codes(oenv))
 346     {
 347         xvgr_view(out, 0.15, 0.15, 0.75, 0.85, oenv);
 348         fprintf(out, "@ legend on\n");
 349         fprintf(out, "@ legend box on\n");
 350         fprintf(out, "@ legend loctype view\n");
 351         fprintf(out, "@ legend %g, %g\n", 0.78, 0.8);
 352         fprintf(out, "@ legend length %d\n", 2);
 353         for (i = 0; (i < nsets); i++)
 354         {
 355             if (setname[i])
 356             {
 357                 if (output_env_get_xvg_format(oenv) == exvgXMGR)
 358                 {
 359                     fprintf(out, "@ legend string %d \"%s\"\n",
 360                             i, xvgrstr(setname[i], oenv, buf, STRLEN));
 361                 }
 362                 else
 363                 {
 364                     fprintf(out, "@ s%d legend \"%s\"\n",
 365                             i, xvgrstr(setname[i], oenv, buf, STRLEN));
 366                 }
 367             }
 368         }
 369     }
 370 }
 371
 372 void xvgr_new_dataset(FILE *out, int nr_first, int nsets,
 373                       const char **setname,
 374                       const gmx_output_env_t *oenv)
 375 {
 376     int  i;
 377     char buf[STRLEN];
 378
 379     if (output_env_get_print_xvgr_codes(oenv))
 380     {
 381         fprintf(out, "@\n");
 382         for (i = 0; (i < nsets); i++)
 383         {
 384             if (setname[i])
 385             {
 386                 if (output_env_get_xvg_format(oenv) == exvgXMGR)
 387                 {
 388                     fprintf(out, "@ legend string %d \"%s\"\n",
 389                             i+nr_first, xvgrstr(setname[i], oenv, buf, STRLEN));
 390                 }
 391                 else
 392                 {
 393                     fprintf(out, "@ s%d legend \"%s\"\n",
 394                             i+nr_first, xvgrstr(setname[i], oenv, buf, STRLEN));
 395                 }
 396             }
 397         }
 398     }
 399     else
 400     {
 401         fprintf(out, "\n");
 402     }
 403 }
 404
 405 void xvgr_line_props(FILE *out, int NrSet, int LineStyle, int LineColor,
 406                      const gmx_output_env_t *oenv)
 407 {
 408     if (output_env_get_print_xvgr_codes(oenv))
 409     {
 410         fprintf(out, "@    with g0\n");
 411         fprintf(out, "@    s%d linestyle %d\n", NrSet, LineStyle);
 412         fprintf(out, "@    s%d color %d\n", NrSet, LineColor);
 413     }
 414 }
 415
 416 static const char *LocTypeStr[] = { "view", "world" };
 417 static const char *BoxFillStr[] = { "none", "color", "pattern" };
 418
 419 void xvgr_box(FILE *out,
 420               int LocType,
 421               real xmin, real ymin, real xmax, real ymax,
 422               int LineStyle, int LineWidth, int LineColor,
 423               int BoxFill, int BoxColor, int BoxPattern, const gmx_output_env_t *oenv)
 424 {
 425     if (output_env_get_print_xvgr_codes(oenv))
 426     {
 427         fprintf(out, "@with box\n");
 428         fprintf(out, "@    box on\n");
 429         fprintf(out, "@    box loctype %s\n", LocTypeStr[LocType]);
 430         fprintf(out, "@    box %g, %g, %g, %g\n", xmin, ymin, xmax, ymax);
 431         fprintf(out, "@    box linestyle %d\n", LineStyle);
 432         fprintf(out, "@    box linewidth %d\n", LineWidth);
 433         fprintf(out, "@    box color %d\n", LineColor);
 434         fprintf(out, "@    box fill %s\n", BoxFillStr[BoxFill]);
 435         fprintf(out, "@    box fill color %d\n", BoxColor);
 436         fprintf(out, "@    box fill pattern %d\n", BoxPattern);
 437         fprintf(out, "@box def\n");
 438     }
 439 }
 440
 441 /* reads a line into ptr, adjusting len and renewing ptr if neccesary */
 442 static char *fgets3(FILE *fp, char **ptr, int *len, int maxlen)
 443 {
 444     int   len_remaining = *len; /* remaining amount of allocated bytes in buf */
 445     int   curp          = 0;    /* current position in buf to read into */
 446
 447     do
 448     {
 449         if (len_remaining < 2)
 450         {
 451             if (*len + STRLEN < maxlen)
 452             {
 453                 /* This line is longer than len characters, let's increase len! */
 454                 *len          += STRLEN;
 455                 len_remaining += STRLEN;
 456                 srenew(*ptr, *len);
 457             }
 458             else
 459             {
 460                 /*something is wrong, we'll just keep reading and return NULL*/
 461                 len_remaining = STRLEN;
 462                 curp          = 0;
 463             }
 464         }
 465         if (fgets(*ptr + curp, len_remaining, fp) == NULL)
 466         {
 467             /* if last line, skip */
 468             return NULL;
 469         }
 470         curp         += len_remaining-1; /* overwrite the nul char in next iteration */
 471         len_remaining = 1;
 472     }
 473     while ((std::strchr(*ptr, '\n') == NULL) && (!feof(fp)));
 474
 475     if (*len + STRLEN >= maxlen)
 476     {
 477         return NULL; /* this line was too long */
 478     }
 479
 480     if (feof(fp))
 481     {
 482         /* We reached EOF before '\n', skip this last line. */
 483         return NULL;
 484     }
 485     {
 486         /* now remove newline */
 487         int slen = std::strlen(*ptr);
 488         if ((*ptr)[slen-1] == '\n')
 489         {
 490             (*ptr)[slen-1] = '\0';
 491         }
 492     }
 493
 494     return *ptr;
 495 }
 496
 497 static int wordcount(char *ptr)
 498 {
 499     int i, n = 0, is[2];
 500     int cur = 0;
 501 #define prev (1-cur)
 502
 503     if (NULL != ptr)
 504     {
 505         for (i = 0; (ptr[i] != '\0'); i++)
 506         {
 507             is[cur] = std::isspace(ptr[i]);
 508             if ((0 == i) && !is[cur])
 509             {
 510                 n++;
 511             }
 512             else if ((i > 0)  && (!is[cur] && is[prev]))
 513             {
 514                 n++;
 515             }
 516             cur = prev;
 517         }
 518     }
 519     return n;
 520 }
 521
 522 static char *read_xvgr_string(const char *line)
 523 {
 524     const char *ptr0, *ptr1;
 525     char       *str;
 526
 527     ptr0 = std::strchr(line, '"');
 528     if (ptr0 != NULL)
 529     {
 530         ptr0++;
 531         ptr1 = std::strchr(ptr0, '"');
 532         if (ptr1 != NULL)
 533         {
 534             str            = gmx_strdup(ptr0);
 535             str[ptr1-ptr0] = '\0';
 536         }
 537         else
 538         {
 539             str = gmx_strdup("");
 540         }
 541     }
 542     else
 543     {
 544         str = gmx_strdup("");
 545     }
 546
 547     return str;
 548 }
 549
 550 int read_xvg_legend(const char *fn, double ***y, int *ny,
 551                     char **subtitle, char ***legend)
 552 {
 553     FILE    *fp;
 554     char    *ptr;
 555     char    *base = NULL;
 556     char    *fmt  = NULL;
 557     int      k, line = 0, nny, nx, maxx, rval, legend_nalloc, set, nchar;
 558     double   lf;
 559     double **yy = NULL;
 560     char    *tmpbuf;
 561     int      len = STRLEN;
 562     *ny  = 0;
 563     nny  = 0;
 564     nx   = 0;
 565     maxx = 0;
 566     fp   = gmx_fio_fopen(fn, "r");
 567
 568     snew(tmpbuf, len);
 569     if (subtitle != NULL)
 570     {
 571         *subtitle = NULL;
 572     }
 573     legend_nalloc = 0;
 574     if (legend != NULL)
 575     {
 576         *legend = NULL;
 577     }
 578
 579     while ((ptr = fgets3(fp, &tmpbuf, &len, 10*STRLEN)) != NULL && ptr[0] != '&')
 580     {
 581         line++;
 582         trim(ptr);
 583         if (ptr[0] == '@')
 584         {
 585             if (legend != NULL)
 586             {
 587                 ptr++;
 588                 trim(ptr);
 589                 set = -1;
 590                 if (std::strncmp(ptr, "subtitle", 8) == 0)
 591                 {
 592                     ptr += 8;
 593                     if (subtitle != NULL)
 594                     {
 595                         *subtitle = read_xvgr_string(ptr);
 596                     }
 597                 }
 598                 else if (std::strncmp(ptr, "legend string", 13) == 0)
 599                 {
 600                     ptr += 13;
 601                     sscanf(ptr, "%d%n", &set, &nchar);
 602                     ptr += nchar;
 603                 }
 604                 else if (ptr[0] == 's')
 605                 {
 606                     ptr++;
 607                     sscanf(ptr, "%d%n", &set, &nchar);
 608                     ptr += nchar;
 609                     trim(ptr);
 610                     if (std::strncmp(ptr, "legend", 6) == 0)
 611                     {
 612                         ptr += 6;
 613                     }
 614                     else
 615                     {
 616                         set = -1;
 617                     }
 618                 }
 619                 if (set >= 0)
 620                 {
 621                     if (set >= legend_nalloc)
 622                     {
 623                         legend_nalloc = set + 1;
 624                         srenew(*legend, legend_nalloc);
 625                         (*legend)[set] = read_xvgr_string(ptr);
 626                     }
 627                 }
 628             }
 629         }
 630         else if (ptr[0] != '#')
 631         {
 632             if (nny == 0)
 633             {
 634                 (*ny) = nny = wordcount(ptr);
 635                 /* fprintf(stderr,"There are %d columns in your file\n",nny);*/
 636                 if (nny == 0)
 637                 {
 638                     return 0;
 639                 }
 640                 snew(yy, nny);
 641                 snew(fmt, 3*nny+1);
 642                 snew(base, 3*nny+1);
 643             }
 644             /* Allocate column space */
 645             if (nx >= maxx)
 646             {
 647                 maxx += 1024;
 648                 for (k = 0; (k < nny); k++)
 649                 {
 650                     srenew(yy[k], maxx);
 651                 }
 652             }
 653             /* Initiate format string */
 654             fmt[0]  = '\0';
 655             base[0] = '\0';
 656
 657             /* fprintf(stderr,"ptr='%s'\n",ptr);*/
 658             for (k = 0; (k < nny); k++)
 659             {
 660                 std::strcpy(fmt, base);
 661                 std::strcat(fmt, "%lf");
 662                 rval = sscanf(ptr, fmt, &lf);
 663                 /* fprintf(stderr,"rval = %d\n",rval);*/
 664                 if ((rval == EOF) || (rval == 0))
 665                 {
 666                     break;
 667                 }
 668                 yy[k][nx] = lf;
 669                 srenew(fmt, 3*(nny+1)+1);
 670                 srenew(base, 3*nny+1);
 671                 std::strcat(base, "%*s");
 672             }
 673             if (k != nny)
 674             {
 675                 fprintf(stderr, "Only %d columns on line %d in file %s\n",
 676                         k, line, fn);
 677                 for (; (k < nny); k++)
 678                 {
 679                     yy[k][nx] = 0.0;
 680                 }
 681             }
 682             nx++;
 683         }
 684     }
 685     gmx_fio_fclose(fp);
 686
 687     *y = yy;
 688     sfree(tmpbuf);
 689     sfree(base);
 690     sfree(fmt);
 691
 692     if (legend_nalloc > 0)
 693     {
 694         if (*ny - 1 > legend_nalloc)
 695         {
 696             assert(legend);
 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 gmx_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 }