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