src/gmxlib/xvgr.c

   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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <string.h>
  43 #include <ctype.h>
  44 #include <time.h>
  45
  46 #ifdef HAVE_SYS_TIME_H
  47 #include <sys/time.h>
  48 #endif
  49
  50 #include "sysstuff.h"
  51 #include "string2.h"
  52 #include "futil.h"
  53 #include "statutil.h"
  54 #include "copyrite.h"
  55 #include "smalloc.h"
  56 #include "xvgr.h"
  57 #include "viewit.h"
  58 #include "vec.h"
  59 #include "gmxfio.h"
  60
  61 /* Portable version of ctime_r implemented in src/gmxlib/string2.c, but we do not want it declared in public installed headers */
  62 char *
  63 gmx_ctime_r(const time_t *clock,char *buf, int n);
  64
  65
  66 gmx_bool output_env_get_print_xvgr_codes(const output_env_t oenv)
  67 {
  68     int xvg_format;
  69
  70     xvg_format = output_env_get_xvg_format(oenv);
  71
  72     return (xvg_format == exvgXMGRACE || xvg_format == exvgXMGR);
  73 }
  74
  75 static char *xvgrstr(const char *gmx,const output_env_t oenv,
  76                      char *buf,int buflen)
  77 {
  78     /* Supported greek letter names and corresponding xmgrace/xmgr symbols */
  79     const char *sym[]={ "beta", "chi", "delta", "eta", "lambda", "mu", "omega", "phi", "psi", "rho", "theta", NULL };
  80     const char symc[]={ 'b',    'c',   'd',     'h',   'l',      'm',  'w',     'f',   'y',   'r',   'q',     '\0' };
  81     int  xvgf;
  82     gmx_bool bXVGR;
  83     int  g,b,i;
  84     char c;
  85
  86     xvgf  = output_env_get_xvg_format(oenv);
  87     bXVGR = (xvgf == exvgXMGRACE || xvgf == exvgXMGR);
  88
  89     g = 0;
  90     b = 0;
  91     while (gmx[g] != '\0')
  92     {
  93         /* Check with the largest string we have ("lambda"), add one for \0 */
  94         if (b + 6 + 1 >= buflen)
  95         {
  96             gmx_fatal(FARGS,"Output buffer length in xvgstr (%d) too small to process xvg input string '%s'",buflen,gmx);
  97         }
  98         if (gmx[g] == '\\')
  99         {
 100             g++;
 101             if (gmx[g] == 's')
 102             {
 103                 /* Subscript */
 104                 if (bXVGR)
 105                 {
 106                     buf[b++] = '\\';
 107                     buf[b++] = 's';
 108                 }
 109                 else
 110                 {
 111                     buf[b++] = '_';
 112                 }
 113                 g++;
 114             }
 115             else if (gmx[g] == 'S')
 116             {
 117                 /* Superscript */
 118                 if (bXVGR)
 119                 {
 120                     buf[b++] = '\\';
 121                     buf[b++] = 'S';
 122                 }
 123                 else
 124                 {
 125                     buf[b++] = '^';
 126                 }
 127                 g++;
 128             }
 129             else if (gmx[g] == 'N')
 130             {
 131                 /* End sub/superscript */
 132                 if (bXVGR)
 133                 {
 134                     buf[b++] = '\\';
 135                     buf[b++] = 'N';
 136                 }
 137                 else
 138                 {
 139                     if (gmx[g+1] != ' ')
 140                     {
 141                         buf[b++] = ' ';
 142                     }
 143                 }
 144                 g++;
 145             }
 146             else if (gmx[g] == '4')
 147             {
 148                 /* Backward compatibility for xmgr normal font "\4" */
 149                 switch (xvgf)
 150                 {
 151                 case exvgXMGRACE:
 152                     sprintf(buf+b,"%s","\\f{}");
 153                     break;
 154                 case exvgXMGR:
 155                     sprintf(buf+b,"%s","\\4");
 156                     break;
 157                 default:
 158                     buf[b] = '\0';
 159                     break;
 160                 }
 161                 g++;
 162                 b = strlen(buf);
 163             }
 164             else if (gmx[g] == '8')
 165             {
 166                 /* Backward compatibility for xmgr symbol font "\8" */
 167                 switch (xvgf)
 168                 {
 169                 case exvgXMGRACE:
 170                     sprintf(buf+b,"%s","\\x");
 171                     break;
 172                 case exvgXMGR:
 173                     sprintf(buf+b,"%s","\\8");
 174                     break;
 175                 default:
 176                     buf[b] = '\0';
 177                     break;
 178                 }
 179                 g++;
 180                 b = strlen(buf);
 181             }
 182             else
 183             {
 184                 /* Check for special symbol */
 185                 i = 0;
 186                 while (sym[i] != NULL &&
 187                        gmx_strncasecmp(sym[i],gmx+g,strlen(sym[i])) != 0)
 188                 {
 189                     i++;
 190                 }
 191                 if (sym[i] != NULL)
 192                 {
 193                     c = symc[i];
 194                     if (isupper(gmx[g]))
 195                     {
 196                         c = toupper(c);
 197                     }
 198                     switch (xvgf)
 199                     {
 200                     case exvgXMGRACE:
 201                         sprintf(buf+b,"%s%c%s","\\x",c,"\\f{}");
 202                         break;
 203                     case exvgXMGR:
 204                         sprintf(buf+b,"%s%c%s","\\8",c,"\\4");
 205                         break;
 206                     default:
 207                         strncat(buf+b,gmx+g,strlen(sym[i]));
 208                         b += strlen(sym[i]);
 209                         if (gmx[g+strlen(sym[i])] != ' ')
 210                         {
 211                             buf[b++] = ' ';
 212                         }
 213                         buf[b] = '\0';
 214                         break;
 215                     }
 216                     g += strlen(sym[i]);
 217                     b  = strlen(buf);
 218                 }
 219                 else
 220                 {
 221                     /* Unknown escape sequence, this should not happen.
 222                      * We do not generate a fatal error, since that might
 223                      * stop otherwise functioning code from working.
 224                      * Copy the backslash to the output and continue processing.
 225                      */
 226                     buf[b++] = '\\';
 227                 }
 228             }
 229         }
 230         else
 231         {
 232             buf[b++] = gmx[g++];
 233         }
 234     }
 235
 236     buf[b++] = '\0';
 237
 238     return buf;
 239 }
 240
 241 void xvgr_header(FILE *fp,const char *title,const char *xaxis,
 242                  const char *yaxis,int exvg_graph_type,
 243                  const output_env_t oenv)
 244 {
 245     char pukestr[100],buf[STRLEN];
 246     time_t t;
 247
 248     if (output_env_get_print_xvgr_codes(oenv))
 249     {
 250         time(&t);
 251         gmx_ctime_r(&t,buf,STRLEN);
 252         fprintf(fp,"# This file was created %s",buf);
 253         fprintf(fp,"# by the following command:\n# %s\n#\n",command_line());
 254         fprintf(fp,"# %s is part of G R O M A C S:\n#\n",ShortProgram());
 255         bromacs(pukestr,99);
 256         fprintf(fp,"# %s\n#\n",pukestr);
 257         fprintf(fp,"@    title \"%s\"\n",xvgrstr(title,oenv,buf,STRLEN));
 258         fprintf(fp,"@    xaxis  label \"%s\"\n",
 259                 xvgrstr(xaxis,oenv,buf,STRLEN));
 260         fprintf(fp,"@    yaxis  label \"%s\"\n",
 261                 xvgrstr(yaxis,oenv,buf,STRLEN));
 262         switch (exvg_graph_type)
 263         {
 264         case exvggtXNY:
 265             if (output_env_get_xvg_format(oenv) == exvgXMGR)
 266             {
 267                 fprintf(fp,"@TYPE nxy\n");
 268             }
 269             else
 270             {
 271                 fprintf(fp,"@TYPE xy\n");
 272             }
 273             break;
 274         case exvggtXYDY:
 275             fprintf(fp,"@TYPE xydy\n");
 276             break;
 277         case exvggtXYDYDY:
 278             fprintf(fp,"@TYPE xydydy\n");
 279             break;
 280         }
 281     }
 282 }
 283
 284 FILE *xvgropen_type(const char *fn,const char *title,const char *xaxis,
 285                     const char *yaxis,int exvg_graph_type,
 286                     const output_env_t oenv)
 287 {
 288     FILE *fp;
 289     time_t t;
 290
 291     fp = gmx_fio_fopen(fn,"w");
 292
 293     xvgr_header(fp,title,xaxis,yaxis,exvg_graph_type,oenv);
 294
 295     return fp;
 296 }
 297
 298 FILE *xvgropen(const char *fn,const char *title,const char *xaxis,
 299                const char *yaxis,const output_env_t oenv)
 300 {
 301     return xvgropen_type(fn,title,xaxis,yaxis,exvggtXNY,oenv);
 302 }
 303
 304 void
 305 xvgrclose(FILE *fp)
 306 {
 307         gmx_fio_fclose(fp);
 308 }
 309
 310 void xvgr_subtitle(FILE *out,const char *subtitle,const output_env_t oenv)
 311 {
 312     char buf[STRLEN];
 313
 314     if (output_env_get_print_xvgr_codes(oenv))
 315     {
 316         fprintf(out,"@ subtitle \"%s\"\n",xvgrstr(subtitle,oenv,buf,STRLEN));
 317     }
 318 }
 319
 320 void xvgr_view(FILE *out,real xmin,real ymin,real xmax,real ymax,
 321                const output_env_t oenv)
 322 {
 323     if (output_env_get_print_xvgr_codes(oenv))
 324     {
 325         fprintf(out,"@ view %g, %g, %g, %g\n",xmin,ymin,xmax,ymax);
 326     }
 327 }
 328
 329 void xvgr_world(FILE *out,real xmin,real ymin,real xmax,real ymax,
 330                 const output_env_t oenv)
 331 {
 332     if (output_env_get_print_xvgr_codes(oenv))
 333     {
 334         fprintf(out,"@ world xmin %g\n"
 335                 "@ world ymin %g\n"
 336                 "@ world xmax %g\n"
 337                 "@ world ymax %g\n",xmin,ymin,xmax,ymax);
 338     }
 339 }
 340
 341 void xvgr_legend(FILE *out,int nsets,const char** setname,
 342                  const output_env_t oenv)
 343 {
 344     int  i;
 345     char buf[STRLEN];
 346
 347     if (output_env_get_print_xvgr_codes(oenv))
 348     {
 349         xvgr_view(out,0.15,0.15,0.75,0.85,oenv);
 350         fprintf(out,"@ legend on\n");
 351         fprintf(out,"@ legend box on\n");
 352         fprintf(out,"@ legend loctype view\n");
 353         fprintf(out,"@ legend %g, %g\n",0.78,0.8);
 354         fprintf(out,"@ legend length %d\n",2);
 355         for(i=0; (i<nsets); i++)
 356         {
 357             if (setname[i]) {
 358                 if (output_env_get_xvg_format(oenv) == exvgXMGR)
 359                 {
 360                     fprintf(out,"@ legend string %d \"%s\"\n",
 361                             i,xvgrstr(setname[i],oenv,buf,STRLEN));
 362                 }
 363                 else
 364                 {
 365                     fprintf(out,"@ s%d legend \"%s\"\n",
 366                             i,xvgrstr(setname[i],oenv,buf,STRLEN));
 367                 }
 368             }
 369         }
 370     }
 371 }
 372
 373 void xvgr_new_dataset(FILE *out, int nr_first, int nsets,
 374                       const char **setname,
 375                       const output_env_t oenv)
 376 {
 377     int i;
 378     char buf[STRLEN];
 379
 380     if (output_env_get_print_xvgr_codes(oenv))
 381     {
 382         fprintf(out,"@\n");
 383         for(i=0; (i<nsets); i++)
 384         {
 385             if (setname[i]) {
 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 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 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     char *p;
 445     int len_remaining=*len; /* remaining amount of allocated bytes in buf */
 446     int curp=0;             /* current position in buf to read into */
 447
 448     do
 449     {
 450         if (len_remaining < 2)
 451         {
 452             if ( *len + STRLEN < maxlen)
 453             {
 454                 /* This line is longer than len characters, let's increase len! */
 455                 *len += STRLEN;
 456                 len_remaining += STRLEN;
 457                 srenew(*ptr, *len);
 458             }
 459             else
 460             {
 461                 /*something is wrong, we'll just keep reading and return NULL*/
 462                 len_remaining = STRLEN;
 463                 curp=0;
 464             }
 465         }
 466         if (fgets(*ptr + curp, len_remaining, fp)==NULL)
 467         {
 468             /* if last line, skip */
 469             return NULL;
 470         }
 471         curp += len_remaining-1; /* overwrite the nul char in next iteration */
 472         len_remaining = 1;
 473     }
 474     while ((strchr(*ptr,'\n') == NULL) && (!feof(fp)));
 475
 476     if ( *len + STRLEN >= maxlen )
 477     {
 478         return NULL; /* this line was too long */
 479     }
 480
 481     if (feof(fp))
 482     {
 483         /* We reached EOF before '\n', skip this last line. */
 484         return NULL;
 485     }
 486     {
 487         /* now remove newline */
 488         int slen = strlen(*ptr);
 489         if ((*ptr)[slen-1] == '\n')
 490         {
 491             (*ptr)[slen-1] = '\0';
 492         }
 493     }
 494
 495     return *ptr;
 496 }
 497
 498 static int wordcount(char *ptr)
 499 {
 500   int i,n,is[2];
 501   int cur=0;
 502 #define prev (1-cur)
 503
 504   if (strlen(ptr) == 0)
 505     return 0;
 506   /* fprintf(stderr,"ptr='%s'\n",ptr); */
 507   n=1;
 508   for(i=0; (ptr[i] != '\0'); i++) {
 509     is[cur] = isspace(ptr[i]);
 510     if ((i > 0)  && (is[cur] && !is[prev]))
 511       n++;
 512     cur=prev;
 513   }
 514   return n;
 515 }
 516
 517 static char *read_xvgr_string(const char *line)
 518 {
 519     const char *ptr0,*ptr1;
 520     char *str;
 521
 522     ptr0 = strchr(line,'"');
 523     if (ptr0 != NULL) {
 524         ptr0++;
 525         ptr1 = strchr(ptr0,'"');
 526         if (ptr1 != NULL)
 527         {
 528             str = strdup(ptr0);
 529             str[ptr1-ptr0] = '\0';
 530         }
 531         else
 532         {
 533             str = strdup("");
 534         }
 535     }
 536     else
 537     {
 538         str = strdup("");
 539     }
 540
 541     return str;
 542 }
 543
 544 int read_xvg_legend(const char *fn,double ***y,int *ny,
 545                     char **subtitle,char ***legend)
 546 {
 547   FILE   *fp;
 548   char   *ptr,*ptr0,*ptr1;
 549   char   *base=NULL;
 550   char   *fmt=NULL;
 551   int    k,line=0,nny,nx,maxx,rval,legend_nalloc,set,nchar;
 552   double lf;
 553   double **yy=NULL;
 554   char  *tmpbuf;
 555   int    len=STRLEN;
 556   *ny  = 0;
 557   nny  = 0;
 558   nx   = 0;
 559   maxx = 0;
 560   fp   = gmx_fio_fopen(fn,"r");
 561
 562   snew(tmpbuf,len);
 563   if (subtitle != NULL) {
 564     *subtitle = NULL;
 565   }
 566   legend_nalloc = 0;
 567   if (legend != NULL) {
 568     *legend = NULL;
 569   }
 570
 571   while ((ptr = fgets3(fp,&tmpbuf,&len,10*STRLEN)) != NULL && ptr[0]!='&') {
 572     line++;
 573     trim(ptr);
 574     if (ptr[0] == '@') {
 575       if (legend != NULL) {
 576         ptr++;
 577         trim(ptr);
 578         set = -1;
 579     if (strncmp(ptr,"subtitle",8) == 0) {
 580           ptr += 8;
 581       if (subtitle != NULL)
 582       {
 583           *subtitle = read_xvgr_string(ptr);
 584       }
 585         } else if (strncmp(ptr,"legend string",13) == 0) {
 586           ptr += 13;
 587           sscanf(ptr,"%d%n",&set,&nchar);
 588           ptr += nchar;
 589         } else if (ptr[0] == 's') {
 590           ptr++;
 591           sscanf(ptr,"%d%n",&set,&nchar);
 592           ptr += nchar;
 593           trim(ptr);
 594           if (strncmp(ptr,"legend",6) == 0) {
 595             ptr += 6;
 596           } else {
 597             set = -1;
 598           }
 599         }
 600         if (set >= 0) {
 601           if (set >= legend_nalloc) {
 602             legend_nalloc = set + 1;
 603             srenew(*legend,legend_nalloc);
 604         (*legend)[set] = read_xvgr_string(ptr);
 605       }
 606     }
 607       }
 608     } else if (ptr[0] != '#') {
 609       if (nny == 0) {
 610         (*ny) = nny = wordcount(ptr);
 611         /* fprintf(stderr,"There are %d columns in your file\n",nny);*/
 612         if (nny == 0)
 613           return 0;
 614         snew(yy,nny);
 615         snew(fmt,3*nny+1);
 616         snew(base,3*nny+1);
 617       }
 618       /* Allocate column space */
 619       if (nx >= maxx) {
 620         maxx+=1024;
 621         for(k=0; (k<nny); k++)
 622           srenew(yy[k],maxx);
 623       }
 624       /* Initiate format string */
 625       fmt[0]  = '\0';
 626       base[0] = '\0';
 627
 628       /* fprintf(stderr,"ptr='%s'\n",ptr);*/
 629       for(k=0; (k<nny); k++) {
 630         strcpy(fmt,base);
 631         strcat(fmt,"%lf");
 632         rval = sscanf(ptr,fmt,&lf);
 633         /* fprintf(stderr,"rval = %d\n",rval);*/
 634         if ((rval == EOF) || (rval == 0))
 635           break;
 636         yy[k][nx] = lf;
 637         srenew(fmt,3*(nny+1)+1);
 638         srenew(base,3*nny+1);
 639         strcat(base,"%*s");
 640       }
 641       if (k != nny) {
 642         fprintf(stderr,"Only %d columns on line %d in file %s\n",
 643                 k,line,fn);
 644         for( ; (k<nny); k++)
 645           yy[k][nx] = 0.0;
 646       }
 647       nx++;
 648     }
 649   }
 650   gmx_fio_fclose(fp);
 651
 652   *y = yy;
 653   sfree(tmpbuf);
 654
 655   if (legend_nalloc > 0) {
 656     if (*ny - 1 > legend_nalloc) {
 657       srenew(*legend,*ny-1);
 658       for(set=legend_nalloc; set<*ny-1; set++) {
 659         (*legend)[set] = NULL;
 660       }
 661     }
 662   }
 663
 664   return nx;
 665 }
 666
 667 int read_xvg(const char *fn,double ***y,int *ny)
 668 {
 669     return read_xvg_legend(fn,y,ny,NULL,NULL);
 670 }
 671
 672 void write_xvg(const char *fn,const char *title,int nx,int ny,real **y,
 673                const char **leg,const output_env_t oenv)
 674 {
 675     FILE *fp;
 676     int  i,j;
 677
 678     fp=xvgropen(fn,title,"X","Y",oenv);
 679     if (leg)
 680         xvgr_legend(fp,ny-1,leg,oenv);
 681     for(i=0; (i<nx); i++) {
 682         for(j=0; (j<ny); j++) {
 683             fprintf(fp,"  %12.5e",y[j][i]);
 684         }
 685         fprintf(fp,"\n");
 686     }
 687     xvgrclose(fp);
 688 }
 689
 690 real **read_xvg_time(const char *fn,
 691                      gmx_bool bHaveT,gmx_bool bTB,real tb,gmx_bool bTE,real te,
 692                      int nsets_in,int *nset,int *nval,real *dt,real **t)
 693 {
 694   FILE   *fp;
 695 #define MAXLINELEN 16384
 696   char line0[MAXLINELEN];
 697   char   *line;
 698   int    t_nalloc,*val_nalloc,a,narg,n,sin,set,nchar;
 699   double dbl,tend=0;
 700   gmx_bool   bEndOfSet,bTimeInRange,bFirstLine=TRUE;
 701   real   **val;
 702
 703   t_nalloc = 0;
 704   *t  = NULL;
 705   val = NULL;
 706   val_nalloc = NULL;
 707   *nset = 0;
 708   *dt = 0;
 709   fp  = gmx_fio_fopen(fn,"r");
 710   for(sin=0; sin<nsets_in; sin++) {
 711     if (nsets_in == 1)
 712       narg = 0;
 713     else
 714       narg = bHaveT ? 2 : 1;
 715     n = 0;
 716     bEndOfSet = FALSE;
 717     while (!bEndOfSet && fgets(line0,MAXLINELEN,fp)) {
 718       line = line0;
 719       /* Remove whitespace */
 720       while (line[0]==' ' || line[0]=='\t')
 721         line++;
 722       bEndOfSet = (line[0] == '&');
 723       if (line[0]!='#' && line[0]!='@' && line[0]!='\n' && !bEndOfSet) {
 724         if (bFirstLine && bHaveT) {
 725           /* Check the first line that should contain data */
 726           a = sscanf(line,"%lf%lf",&dbl,&dbl);
 727           if (a == 0)
 728             gmx_fatal(FARGS,"Expected a number in %s on line:\n%s",fn,line0);
 729           else if (a == 1) {
 730             fprintf(stderr,"Found only 1 number on line, "
 731                     "assuming no time is present.\n");
 732             bHaveT = FALSE;
 733             if (nsets_in > 1)
 734               narg = 1;
 735           }
 736         }
 737
 738         a = 0;
 739         bTimeInRange = TRUE;
 740         while ((a<narg || (nsets_in==1 && n==0)) &&
 741                sscanf(line,"%lf%n",&dbl,&nchar)==1 && bTimeInRange) {
 742           /* Use set=-1 as the time "set" */
 743           if (sin) {
 744             if (!bHaveT || (a>0))
 745               set = sin;
 746             else
 747               set = -1;
 748           } else {
 749             if (!bHaveT)
 750               set = a;
 751             else
 752               set = a-1;
 753           }
 754           if (set==-1 && ((bTB && dbl<tb) || (bTE && dbl>te)))
 755             bTimeInRange = FALSE;
 756
 757           if (bTimeInRange) {
 758             if (n==0) {
 759               if (nsets_in == 1)
 760                 narg++;
 761               if (set >= 0) {
 762                 *nset = set+1;
 763                 srenew(val,*nset);
 764                 srenew(val_nalloc,*nset);
 765                 val_nalloc[set] = 0;
 766                 val[set] = NULL;
 767               }
 768             }
 769             if (set == -1) {
 770               if (sin == 0) {
 771                 if (n >= t_nalloc) {
 772                   t_nalloc = over_alloc_small(n);
 773                   srenew(*t,t_nalloc);
 774                 }
 775                 (*t)[n] = dbl;
 776               }
 777               /* else we should check the time of the next sets with set 0 */
 778             } else {
 779               if (n >= val_nalloc[set]) {
 780                 val_nalloc[set] = over_alloc_small(n);
 781                 srenew(val[set],val_nalloc[set]);
 782               }
 783               val[set][n] = (real)dbl;
 784             }
 785           }
 786           a++;
 787           line += nchar;
 788         }
 789         if (line0[strlen(line0)-1] != '\n') {
 790           fprintf(stderr,"File %s does not end with a newline, ignoring the last line\n",fn);
 791         } else if (bTimeInRange) {
 792           if (a == 0) {
 793             fprintf(stderr,"Ignoring invalid line in %s:\n%s",fn,line0);
 794           } else {
 795             if (a != narg)
 796               fprintf(stderr,"Invalid line in %s:\n%s"
 797                       "Using zeros for the last %d sets\n",
 798                       fn,line0,narg-a);
 799             n++;
 800           }
 801         }
 802         if (a > 0)
 803           bFirstLine = FALSE;
 804       }
 805     }
 806     if (sin==0) {
 807       *nval = n;
 808       if (!bHaveT) {
 809         snew(*t,n);
 810         for(a=0; a<n; a++)
 811           (*t)[a] = a;
 812       }
 813       if (n > 1)
 814         *dt = (real)((*t)[n-1]-(*t)[0])/(n-1.0);
 815       else
 816         *dt = 1;
 817     } else {
 818       if (n < *nval) {
 819         fprintf(stderr,"Set %d is shorter (%d) than the previous set (%d)\n",
 820                 sin+1,n,*nval);
 821         *nval = n;
 822         fprintf(stderr,"Will use only the first %d points of every set\n",
 823                 *nval);
 824       }
 825     }
 826   }
 827   gmx_fio_fclose(fp);
 828
 829   sfree(val_nalloc);
 830
 831   return val;
 832 }
 833