src/gromacs/fileio/xvgr.cpp

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