src/gromacs/gmxlib/xvgr.cpp

   1 /*  -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROningen Mixture of Alchemy and Childrens' Stories
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <string.h>
  41 #include <ctype.h>
  42 #include <time.h>
  43
  44 #ifdef HAVE_SYS_TIME_H
  45 #include <sys/time.h>
  46 #endif
  47
  48 #include "sysstuff.h"
  49 #include "string2.h"
  50 #include "futil.h"
  51 #include "statutil.h"
  52 #include "copyrite.h"
  53 #include "smalloc.h"
  54 #include "xvgr.h"
  55 #include "viewit.h"
  56 #include "vec.h"
  57 #include "gmxfio.h"
  58
  59 gmx_bool output_env_get_print_xvgr_codes(const 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 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 = 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, strlen(sym[i])) != 0)
 181                 {
 182                     i++;
 183                 }
 184                 if (sym[i] != NULL)
 185                 {
 186                     c = symc[i];
 187                     if (isupper(gmx[g]))
 188                     {
 189                         c = 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                             strncat(buf+b, gmx+g, strlen(sym[i]));
 201                             b += strlen(sym[i]);
 202                             if (gmx[g+strlen(sym[i])] != ' ')
 203                             {
 204                                 buf[b++] = ' ';
 205                             }
 206                             buf[b] = '\0';
 207                             break;
 208                     }
 209                     g += strlen(sym[i]);
 210                     b  = 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 output_env_t oenv)
 237 {
 238     char   pukestr[100], buf[STRLEN];
 239     time_t t;
 240
 241     if (output_env_get_print_xvgr_codes(oenv))
 242     {
 243         time(&t);
 244         gmx_ctime_r(&t, buf, STRLEN);
 245         fprintf(fp, "# This file was created %s", buf);
 246         fprintf(fp, "# by the following command:\n# %s\n#\n", command_line());
 247         fprintf(fp, "# %s is part of G R O M A C S:\n#\n", ShortProgram());
 248         bromacs(pukestr, 99);
 249         fprintf(fp, "# %s\n#\n", pukestr);
 250         fprintf(fp, "@    title \"%s\"\n", xvgrstr(title, oenv, buf, STRLEN));
 251         fprintf(fp, "@    xaxis  label \"%s\"\n",
 252                 xvgrstr(xaxis, oenv, buf, STRLEN));
 253         fprintf(fp, "@    yaxis  label \"%s\"\n",
 254                 xvgrstr(yaxis, oenv, buf, STRLEN));
 255         switch (exvg_graph_type)
 256         {
 257             case exvggtXNY:
 258                 if (output_env_get_xvg_format(oenv) == exvgXMGR)
 259                 {
 260                     fprintf(fp, "@TYPE nxy\n");
 261                 }
 262                 else
 263                 {
 264                     fprintf(fp, "@TYPE xy\n");
 265                 }
 266                 break;
 267             case exvggtXYDY:
 268                 fprintf(fp, "@TYPE xydy\n");
 269                 break;
 270             case exvggtXYDYDY:
 271                 fprintf(fp, "@TYPE xydydy\n");
 272                 break;
 273         }
 274     }
 275 }
 276
 277 FILE *xvgropen_type(const char *fn, const char *title, const char *xaxis,
 278                     const char *yaxis, int exvg_graph_type,
 279                     const output_env_t oenv)
 280 {
 281     FILE  *fp;
 282
 283     fp = gmx_fio_fopen(fn, "w");
 284
 285     xvgr_header(fp, title, xaxis, yaxis, exvg_graph_type, oenv);
 286
 287     return fp;
 288 }
 289
 290 FILE *xvgropen(const char *fn, const char *title, const char *xaxis,
 291                const char *yaxis, const output_env_t oenv)
 292 {
 293     return xvgropen_type(fn, title, xaxis, yaxis, exvggtXNY, oenv);
 294 }
 295
 296 void
 297 xvgrclose(FILE *fp)
 298 {
 299     gmx_fio_fclose(fp);
 300 }
 301
 302 void xvgr_subtitle(FILE *out, const char *subtitle, const output_env_t oenv)
 303 {
 304     char buf[STRLEN];
 305
 306     if (output_env_get_print_xvgr_codes(oenv))
 307     {
 308         fprintf(out, "@ subtitle \"%s\"\n", xvgrstr(subtitle, oenv, buf, STRLEN));
 309     }
 310 }
 311
 312 void xvgr_view(FILE *out, real xmin, real ymin, real xmax, real ymax,
 313                const output_env_t oenv)
 314 {
 315     if (output_env_get_print_xvgr_codes(oenv))
 316     {
 317         fprintf(out, "@ view %g, %g, %g, %g\n", xmin, ymin, xmax, ymax);
 318     }
 319 }
 320
 321 void xvgr_world(FILE *out, real xmin, real ymin, real xmax, real ymax,
 322                 const output_env_t oenv)
 323 {
 324     if (output_env_get_print_xvgr_codes(oenv))
 325     {
 326         fprintf(out, "@ world xmin %g\n"
 327                 "@ world ymin %g\n"
 328                 "@ world xmax %g\n"
 329                 "@ world ymax %g\n", xmin, ymin, xmax, ymax);
 330     }
 331 }
 332
 333 void xvgr_legend(FILE *out, int nsets, const char** setname,
 334                  const output_env_t oenv)
 335 {
 336     int  i;
 337     char buf[STRLEN];
 338
 339     if (output_env_get_print_xvgr_codes(oenv))
 340     {
 341         xvgr_view(out, 0.15, 0.15, 0.75, 0.85, oenv);
 342         fprintf(out, "@ legend on\n");
 343         fprintf(out, "@ legend box on\n");
 344         fprintf(out, "@ legend loctype view\n");
 345         fprintf(out, "@ legend %g, %g\n", 0.78, 0.8);
 346         fprintf(out, "@ legend length %d\n", 2);
 347         for (i = 0; (i < nsets); i++)
 348         {
 349             if (setname[i])
 350             {
 351                 if (output_env_get_xvg_format(oenv) == exvgXMGR)
 352                 {
 353                     fprintf(out, "@ legend string %d \"%s\"\n",
 354                             i, xvgrstr(setname[i], oenv, buf, STRLEN));
 355                 }
 356                 else
 357                 {
 358                     fprintf(out, "@ s%d legend \"%s\"\n",
 359                             i, xvgrstr(setname[i], oenv, buf, STRLEN));
 360                 }
 361             }
 362         }
 363     }
 364 }
 365
 366 void xvgr_new_dataset(FILE *out, int nr_first, int nsets,
 367                       const char **setname,
 368                       const output_env_t oenv)
 369 {
 370     int  i;
 371     char buf[STRLEN];
 372
 373     if (output_env_get_print_xvgr_codes(oenv))
 374     {
 375         fprintf(out, "@\n");
 376         for (i = 0; (i < nsets); i++)
 377         {
 378             if (setname[i])
 379             {
 380                 if (output_env_get_xvg_format(oenv) == exvgXMGR)
 381                 {
 382                     fprintf(out, "@ legend string %d \"%s\"\n",
 383                             i+nr_first, xvgrstr(setname[i], oenv, buf, STRLEN));
 384                 }
 385                 else
 386                 {
 387                     fprintf(out, "@ s%d legend \"%s\"\n",
 388                             i+nr_first, xvgrstr(setname[i], oenv, buf, STRLEN));
 389                 }
 390             }
 391         }
 392     }
 393     else
 394     {
 395         fprintf(out, "\n");
 396     }
 397 }
 398
 399 void xvgr_line_props(FILE *out, int NrSet, int LineStyle, int LineColor,
 400                      const output_env_t oenv)
 401 {
 402     if (output_env_get_print_xvgr_codes(oenv))
 403     {
 404         fprintf(out, "@    with g0\n");
 405         fprintf(out, "@    s%d linestyle %d\n", NrSet, LineStyle);
 406         fprintf(out, "@    s%d color %d\n", NrSet, LineColor);
 407     }
 408 }
 409
 410 static const char *LocTypeStr[] = { "view", "world" };
 411 static const char *BoxFillStr[] = { "none", "color", "pattern" };
 412
 413 void xvgr_box(FILE *out,
 414               int LocType,
 415               real xmin, real ymin, real xmax, real ymax,
 416               int LineStyle, int LineWidth, int LineColor,
 417               int BoxFill, int BoxColor, int BoxPattern, const output_env_t oenv)
 418 {
 419     if (output_env_get_print_xvgr_codes(oenv))
 420     {
 421         fprintf(out, "@with box\n");
 422         fprintf(out, "@    box on\n");
 423         fprintf(out, "@    box loctype %s\n", LocTypeStr[LocType]);
 424         fprintf(out, "@    box %g, %g, %g, %g\n", xmin, ymin, xmax, ymax);
 425         fprintf(out, "@    box linestyle %d\n", LineStyle);
 426         fprintf(out, "@    box linewidth %d\n", LineWidth);
 427         fprintf(out, "@    box color %d\n", LineColor);
 428         fprintf(out, "@    box fill %s\n", BoxFillStr[BoxFill]);
 429         fprintf(out, "@    box fill color %d\n", BoxColor);
 430         fprintf(out, "@    box fill pattern %d\n", BoxPattern);
 431         fprintf(out, "@box def\n");
 432     }
 433 }
 434
 435 /* reads a line into ptr, adjusting len and renewing ptr if neccesary */
 436 static char *fgets3(FILE *fp, char **ptr, int *len, int maxlen)
 437 {
 438     int   len_remaining = *len; /* remaining amount of allocated bytes in buf */
 439     int   curp          = 0;    /* current position in buf to read into */
 440
 441     do
 442     {
 443         if (len_remaining < 2)
 444         {
 445             if (*len + STRLEN < maxlen)
 446             {
 447                 /* This line is longer than len characters, let's increase len! */
 448                 *len          += STRLEN;
 449                 len_remaining += STRLEN;
 450                 srenew(*ptr, *len);
 451             }
 452             else
 453             {
 454                 /*something is wrong, we'll just keep reading and return NULL*/
 455                 len_remaining = STRLEN;
 456                 curp          = 0;
 457             }
 458         }
 459         if (fgets(*ptr + curp, len_remaining, fp) == NULL)
 460         {
 461             /* if last line, skip */
 462             return NULL;
 463         }
 464         curp         += len_remaining-1; /* overwrite the nul char in next iteration */
 465         len_remaining = 1;
 466     }
 467     while ((strchr(*ptr, '\n') == NULL) && (!feof(fp)));
 468
 469     if (*len + STRLEN >= maxlen)
 470     {
 471         return NULL; /* this line was too long */
 472     }
 473
 474     if (feof(fp))
 475     {
 476         /* We reached EOF before '\n', skip this last line. */
 477         return NULL;
 478     }
 479     {
 480         /* now remove newline */
 481         int slen = strlen(*ptr);
 482         if ((*ptr)[slen-1] == '\n')
 483         {
 484             (*ptr)[slen-1] = '\0';
 485         }
 486     }
 487
 488     return *ptr;
 489 }
 490
 491 static int wordcount(char *ptr)
 492 {
 493     int i, n, is[2];
 494     int cur = 0;
 495 #define prev (1-cur)
 496
 497     if (strlen(ptr) == 0)
 498     {
 499         return 0;
 500     }
 501     /* fprintf(stderr,"ptr='%s'\n",ptr); */
 502     n = 1;
 503     for (i = 0; (ptr[i] != '\0'); i++)
 504     {
 505         is[cur] = isspace(ptr[i]);
 506         if ((i > 0)  && (is[cur] && !is[prev]))
 507         {
 508             n++;
 509         }
 510         cur = prev;
 511     }
 512     return n;
 513 }
 514
 515 static char *read_xvgr_string(const char *line)
 516 {
 517     const char *ptr0, *ptr1;
 518     char       *str;
 519
 520     ptr0 = strchr(line, '"');
 521     if (ptr0 != NULL)
 522     {
 523         ptr0++;
 524         ptr1 = strchr(ptr0, '"');
 525         if (ptr1 != NULL)
 526         {
 527             str            = strdup(ptr0);
 528             str[ptr1-ptr0] = '\0';
 529         }
 530         else
 531         {
 532             str = strdup("");
 533         }
 534     }
 535     else
 536     {
 537         str = strdup("");
 538     }
 539
 540     return str;
 541 }
 542
 543 int read_xvg_legend(const char *fn, double ***y, int *ny,
 544                     char **subtitle, char ***legend)
 545 {
 546     FILE    *fp;
 547     char    *ptr;
 548     char    *base = NULL;
 549     char    *fmt  = NULL;
 550     int      k, line = 0, nny, nx, maxx, rval, legend_nalloc, set, nchar;
 551     double   lf;
 552     double **yy = NULL;
 553     char    *tmpbuf;
 554     int      len = STRLEN;
 555     *ny  = 0;
 556     nny  = 0;
 557     nx   = 0;
 558     maxx = 0;
 559     fp   = gmx_fio_fopen(fn, "r");
 560
 561     snew(tmpbuf, len);
 562     if (subtitle != NULL)
 563     {
 564         *subtitle = NULL;
 565     }
 566     legend_nalloc = 0;
 567     if (legend != NULL)
 568     {
 569         *legend = NULL;
 570     }
 571
 572     while ((ptr = fgets3(fp, &tmpbuf, &len, 10*STRLEN)) != NULL && ptr[0] != '&')
 573     {
 574         line++;
 575         trim(ptr);
 576         if (ptr[0] == '@')
 577         {
 578             if (legend != NULL)
 579             {
 580                 ptr++;
 581                 trim(ptr);
 582                 set = -1;
 583                 if (strncmp(ptr, "subtitle", 8) == 0)
 584                 {
 585                     ptr += 8;
 586                     if (subtitle != NULL)
 587                     {
 588                         *subtitle = read_xvgr_string(ptr);
 589                     }
 590                 }
 591                 else if (strncmp(ptr, "legend string", 13) == 0)
 592                 {
 593                     ptr += 13;
 594                     sscanf(ptr, "%d%n", &set, &nchar);
 595                     ptr += nchar;
 596                 }
 597                 else if (ptr[0] == 's')
 598                 {
 599                     ptr++;
 600                     sscanf(ptr, "%d%n", &set, &nchar);
 601                     ptr += nchar;
 602                     trim(ptr);
 603                     if (strncmp(ptr, "legend", 6) == 0)
 604                     {
 605                         ptr += 6;
 606                     }
 607                     else
 608                     {
 609                         set = -1;
 610                     }
 611                 }
 612                 if (set >= 0)
 613                 {
 614                     if (set >= legend_nalloc)
 615                     {
 616                         legend_nalloc = set + 1;
 617                         srenew(*legend, legend_nalloc);
 618                         (*legend)[set] = read_xvgr_string(ptr);
 619                     }
 620                 }
 621             }
 622         }
 623         else if (ptr[0] != '#')
 624         {
 625             if (nny == 0)
 626             {
 627                 (*ny) = nny = wordcount(ptr);
 628                 /* fprintf(stderr,"There are %d columns in your file\n",nny);*/
 629                 if (nny == 0)
 630                 {
 631                     return 0;
 632                 }
 633                 snew(yy, nny);
 634                 snew(fmt, 3*nny+1);
 635                 snew(base, 3*nny+1);
 636             }
 637             /* Allocate column space */
 638             if (nx >= maxx)
 639             {
 640                 maxx += 1024;
 641                 for (k = 0; (k < nny); k++)
 642                 {
 643                     srenew(yy[k], maxx);
 644                 }
 645             }
 646             /* Initiate format string */
 647             fmt[0]  = '\0';
 648             base[0] = '\0';
 649
 650             /* fprintf(stderr,"ptr='%s'\n",ptr);*/
 651             for (k = 0; (k < nny); k++)
 652             {
 653                 strcpy(fmt, base);
 654                 strcat(fmt, "%lf");
 655                 rval = sscanf(ptr, fmt, &lf);
 656                 /* fprintf(stderr,"rval = %d\n",rval);*/
 657                 if ((rval == EOF) || (rval == 0))
 658                 {
 659                     break;
 660                 }
 661                 yy[k][nx] = lf;
 662                 srenew(fmt, 3*(nny+1)+1);
 663                 srenew(base, 3*nny+1);
 664                 strcat(base, "%*s");
 665             }
 666             if (k != nny)
 667             {
 668                 fprintf(stderr, "Only %d columns on line %d in file %s\n",
 669                         k, line, fn);
 670                 for (; (k < nny); k++)
 671                 {
 672                     yy[k][nx] = 0.0;
 673                 }
 674             }
 675             nx++;
 676         }
 677     }
 678     gmx_fio_fclose(fp);
 679
 680     *y = yy;
 681     sfree(tmpbuf);
 682
 683     if (legend_nalloc > 0)
 684     {
 685         if (*ny - 1 > legend_nalloc)
 686         {
 687             srenew(*legend, *ny-1);
 688             for (set = legend_nalloc; set < *ny-1; set++)
 689             {
 690                 (*legend)[set] = NULL;
 691             }
 692         }
 693     }
 694
 695     return nx;
 696 }
 697
 698 int read_xvg(const char *fn, double ***y, int *ny)
 699 {
 700     return read_xvg_legend(fn, y, ny, NULL, NULL);
 701 }
 702
 703 void write_xvg(const char *fn, const char *title, int nx, int ny, real **y,
 704                const char **leg, const output_env_t oenv)
 705 {
 706     FILE *fp;
 707     int   i, j;
 708
 709     fp = xvgropen(fn, title, "X", "Y", oenv);
 710     if (leg)
 711     {
 712         xvgr_legend(fp, ny-1, leg, oenv);
 713     }
 714     for (i = 0; (i < nx); i++)
 715     {
 716         for (j = 0; (j < ny); j++)
 717         {
 718             fprintf(fp, "  %12.5e", y[j][i]);
 719         }
 720         fprintf(fp, "\n");
 721     }
 722     xvgrclose(fp);
 723 }
 724
 725 real **read_xvg_time(const char *fn,
 726                      gmx_bool bHaveT, gmx_bool bTB, real tb, gmx_bool bTE, real te,
 727                      int nsets_in, int *nset, int *nval, real *dt, real **t)
 728 {
 729     FILE      *fp;
 730 #define MAXLINELEN 16384
 731     char       line0[MAXLINELEN];
 732     char      *line;
 733     int        t_nalloc, *val_nalloc, a, narg, n, sin, set, nchar;
 734     double     dbl;
 735     gmx_bool   bEndOfSet, bTimeInRange, bFirstLine = TRUE;
 736     real     **val;
 737
 738     t_nalloc   = 0;
 739     *t         = NULL;
 740     val        = NULL;
 741     val_nalloc = NULL;
 742     *nset      = 0;
 743     *dt        = 0;
 744     fp         = gmx_fio_fopen(fn, "r");
 745     for (sin = 0; sin < nsets_in; sin++)
 746     {
 747         if (nsets_in == 1)
 748         {
 749             narg = 0;
 750         }
 751         else
 752         {
 753             narg = bHaveT ? 2 : 1;
 754         }
 755         n         = 0;
 756         bEndOfSet = FALSE;
 757         while (!bEndOfSet && fgets(line0, MAXLINELEN, fp))
 758         {
 759             line = line0;
 760             /* Remove whitespace */
 761             while (line[0] == ' ' || line[0] == '\t')
 762             {
 763                 line++;
 764             }
 765             bEndOfSet = (line[0] == '&');
 766             if (line[0] != '#' && line[0] != '@' && line[0] != '\n' && !bEndOfSet)
 767             {
 768                 if (bFirstLine && bHaveT)
 769                 {
 770                     /* Check the first line that should contain data */
 771                     a = sscanf(line, "%lf%lf", &dbl, &dbl);
 772                     if (a == 0)
 773                     {
 774                         gmx_fatal(FARGS, "Expected a number in %s on line:\n%s", fn, line0);
 775                     }
 776                     else if (a == 1)
 777                     {
 778                         fprintf(stderr, "Found only 1 number on line, "
 779                                 "assuming no time is present.\n");
 780                         bHaveT = FALSE;
 781                         if (nsets_in > 1)
 782                         {
 783                             narg = 1;
 784                         }
 785                     }
 786                 }
 787
 788                 a            = 0;
 789                 bTimeInRange = TRUE;
 790                 while ((a < narg || (nsets_in == 1 && n == 0)) &&
 791                        sscanf(line, "%lf%n", &dbl, &nchar) == 1 && bTimeInRange)
 792                 {
 793                     /* Use set=-1 as the time "set" */
 794                     if (sin)
 795                     {
 796                         if (!bHaveT || (a > 0))
 797                         {
 798                             set = sin;
 799                         }
 800                         else
 801                         {
 802                             set = -1;
 803                         }
 804                     }
 805                     else
 806                     {
 807                         if (!bHaveT)
 808                         {
 809                             set = a;
 810                         }
 811                         else
 812                         {
 813                             set = a-1;
 814                         }
 815                     }
 816                     if (set == -1 && ((bTB && dbl < tb) || (bTE && dbl > te)))
 817                     {
 818                         bTimeInRange = FALSE;
 819                     }
 820
 821                     if (bTimeInRange)
 822                     {
 823                         if (n == 0)
 824                         {
 825                             if (nsets_in == 1)
 826                             {
 827                                 narg++;
 828                             }
 829                             if (set >= 0)
 830                             {
 831                                 *nset = set+1;
 832                                 srenew(val, *nset);
 833                                 srenew(val_nalloc, *nset);
 834                                 val_nalloc[set] = 0;
 835                                 val[set]        = NULL;
 836                             }
 837                         }
 838                         if (set == -1)
 839                         {
 840                             if (sin == 0)
 841                             {
 842                                 if (n >= t_nalloc)
 843                                 {
 844                                     t_nalloc = over_alloc_small(n);
 845                                     srenew(*t, t_nalloc);
 846                                 }
 847                                 (*t)[n] = dbl;
 848                             }
 849                             /* else we should check the time of the next sets with set 0 */
 850                         }
 851                         else
 852                         {
 853                             if (n >= val_nalloc[set])
 854                             {
 855                                 val_nalloc[set] = over_alloc_small(n);
 856                                 srenew(val[set], val_nalloc[set]);
 857                             }
 858                             val[set][n] = (real)dbl;
 859                         }
 860                     }
 861                     a++;
 862                     line += nchar;
 863                 }
 864                 if (line0[strlen(line0)-1] != '\n')
 865                 {
 866                     fprintf(stderr, "File %s does not end with a newline, ignoring the last line\n", fn);
 867                 }
 868                 else if (bTimeInRange)
 869                 {
 870                     if (a == 0)
 871                     {
 872                         fprintf(stderr, "Ignoring invalid line in %s:\n%s", fn, line0);
 873                     }
 874                     else
 875                     {
 876                         if (a != narg)
 877                         {
 878                             fprintf(stderr, "Invalid line in %s:\n%s"
 879                                     "Using zeros for the last %d sets\n",
 880                                     fn, line0, narg-a);
 881                         }
 882                         n++;
 883                     }
 884                 }
 885                 if (a > 0)
 886                 {
 887                     bFirstLine = FALSE;
 888                 }
 889             }
 890         }
 891         if (sin == 0)
 892         {
 893             *nval = n;
 894             if (!bHaveT)
 895             {
 896                 snew(*t, n);
 897                 for (a = 0; a < n; a++)
 898                 {
 899                     (*t)[a] = a;
 900                 }
 901             }
 902             if (n > 1)
 903             {
 904                 *dt = (real)((*t)[n-1]-(*t)[0])/(n-1.0);
 905             }
 906             else
 907             {
 908                 *dt = 1;
 909             }
 910         }
 911         else
 912         {
 913             if (n < *nval)
 914             {
 915                 fprintf(stderr, "Set %d is shorter (%d) than the previous set (%d)\n",
 916                         sin+1, n, *nval);
 917                 *nval = n;
 918                 fprintf(stderr, "Will use only the first %d points of every set\n",
 919                         *nval);
 920             }
 921         }
 922     }
 923     gmx_fio_fclose(fp);
 924
 925     sfree(val_nalloc);
 926
 927     return val;
 928 }