src/gmxlib/xvgr.c

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