src/gromacs/fileio/matio.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <stdio.h>
  42 #include <ctype.h>
  43
  44 #include <algorithm>
  45
  46 #include "gromacs/utility/futil.h"
  47 #include "gromacs/utility/cstringutil.h"
  48 #include "macros.h"
  49 #include "gromacs/utility/smalloc.h"
  50 #include "gromacs/utility/fatalerror.h"
  51 #include "matio.h"
  52 #include "gmxfio.h"
  53 #include "gromacs/math/utilities.h"
  54 #include "copyrite.h"
  55
  56 #include "gromacs/utility/exceptions.h"
  57 #include "gromacs/utility/programcontext.h"
  58
  59 #define round(a) (int)(a+0.5)
  60
  61 static const char mapper[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<.>/?";
  62 #define NMAP (long int)strlen(mapper)
  63
  64 #define MAX_XPM_LINELENGTH 4096
  65
  66 real **mk_matrix(int nx, int ny, gmx_bool b1D)
  67 {
  68     int    i;
  69     real **m;
  70
  71     snew(m, nx);
  72     if (b1D)
  73     {
  74         snew(m[0], nx*ny);
  75     }
  76
  77     for (i = 0; (i < nx); i++)
  78     {
  79         if (b1D)
  80         {
  81             m[i] = &(m[0][i*ny]);
  82         }
  83         else
  84         {
  85             snew(m[i], ny);
  86         }
  87     }
  88
  89     return m;
  90 }
  91
  92 void done_matrix(int nx, real ***m)
  93 {
  94     int i;
  95
  96     for (i = 0; (i < nx); i++)
  97     {
  98         sfree((*m)[i]);
  99     }
 100     sfree(*m);
 101     *m = NULL;
 102 }
 103
 104 void clear_matrix(int nx, int ny, real **m)
 105 {
 106     int x, y;
 107
 108     for (x = 0; x < nx; x++)
 109     {
 110         for (y = 0; y < ny; y++)
 111         {
 112             m[x][y] = 0;
 113         }
 114     }
 115 }
 116
 117 gmx_bool matelmt_cmp(t_xpmelmt e1, t_xpmelmt e2)
 118 {
 119     return (e1.c1 == e2.c1) && (e1.c2 == e2.c2);
 120 }
 121
 122 t_matelmt searchcmap(int n, t_mapping map[], t_xpmelmt c)
 123 {
 124     int i;
 125
 126     for (i = 0; (i < n); i++)
 127     {
 128         if (matelmt_cmp(map[i].code, c))
 129         {
 130             return i;
 131         }
 132     }
 133
 134     return -1;
 135 }
 136
 137 int getcmap(FILE *in, const char *fn, t_mapping **map)
 138 {
 139     int        i, n;
 140     char       line[STRLEN];
 141     char       code[STRLEN], desc[STRLEN];
 142     double     r, g, b;
 143     t_mapping *m;
 144
 145     if (fgets2(line, STRLEN-1, in) == NULL)
 146     {
 147         gmx_fatal(FARGS, "Not enough lines in colormap file %s"
 148                   "(just wanted to read number of entries)", fn);
 149     }
 150     sscanf(line, "%d", &n);
 151     snew(m, n);
 152     for (i = 0; (i < n); i++)
 153     {
 154         if (fgets2(line, STRLEN-1, in) == NULL)
 155         {
 156             gmx_fatal(FARGS, "Not enough lines in colormap file %s"
 157                       "(should be %d, found only %d)", fn, n+1, i);
 158         }
 159         sscanf(line, "%s%s%lf%lf%lf", code, desc, &r, &g, &b);
 160         m[i].code.c1 = code[0];
 161         m[i].code.c2 = 0;
 162         m[i].desc    = strdup(desc);
 163         m[i].rgb.r   = r;
 164         m[i].rgb.g   = g;
 165         m[i].rgb.b   = b;
 166     }
 167     *map = m;
 168
 169     return n;
 170 }
 171
 172 int readcmap(const char *fn, t_mapping **map)
 173 {
 174     FILE      *in;
 175     int        n;
 176
 177     in = libopen(fn);
 178     n  = getcmap(in, fn, map);
 179     gmx_fio_fclose(in);
 180
 181     return n;
 182 }
 183
 184 void printcmap(FILE *out, int n, t_mapping map[])
 185 {
 186     int i;
 187
 188     fprintf(out, "%d\n", n);
 189     for (i = 0; (i < n); i++)
 190     {
 191         fprintf(out, "%c%c  %20s  %10g  %10g  %10g\n",
 192                 map[i].code.c1 ? map[i].code.c1 : ' ',
 193                 map[i].code.c2 ? map[i].code.c2 : ' ',
 194                 map[i].desc, map[i].rgb.r, map[i].rgb.g, map[i].rgb.b);
 195     }
 196 }
 197
 198 void writecmap(const char *fn, int n, t_mapping map[])
 199 {
 200     FILE *out;
 201
 202     out = gmx_fio_fopen(fn, "w");
 203     printcmap(out, n, map);
 204     gmx_fio_fclose(out);
 205 }
 206
 207 static char *fgetline(char **line, int llmax, int *llalloc, FILE *in)
 208 {
 209     char *fg;
 210
 211     if (llmax > *llalloc)
 212     {
 213         srenew(*line, llmax+1);
 214         *llalloc = llmax;
 215     }
 216     fg = fgets(*line, llmax, in);
 217     trim(*line);
 218
 219     return fg;
 220 }
 221
 222 void skipstr(char *line)
 223 {
 224     int i, c;
 225
 226     ltrim(line);
 227     c = 0;
 228     while ((line[c] != ' ') && (line[c] != '\0'))
 229     {
 230         c++;
 231     }
 232     i = c;
 233     while (line[c] != '\0')
 234     {
 235         line[c-i] = line[c];
 236         c++;
 237     }
 238     line[c-i] = '\0';
 239 }
 240
 241 char *line2string(char **line)
 242 {
 243     int i;
 244
 245     if (*line != NULL)
 246     {
 247         while (((*line)[0] != '\"' ) && ( (*line)[0] != '\0' ))
 248         {
 249             (*line)++;
 250         }
 251
 252         if ((*line)[0] != '\"')
 253         {
 254             return NULL;
 255         }
 256         (*line)++;
 257
 258         i = 0;
 259         while (( (*line)[i] != '\"' ) && ( (*line)[i] != '\0' ))
 260         {
 261             i++;
 262         }
 263
 264         if ((*line)[i] != '\"')
 265         {
 266             *line = NULL;
 267         }
 268         else
 269         {
 270             (*line)[i] = 0;
 271         }
 272     }
 273
 274     return *line;
 275 }
 276
 277 void parsestring(char *line, const char *label, char *string)
 278 {
 279     if (strstr(line, label))
 280     {
 281         if (strstr(line, label) < strchr(line, '\"'))
 282         {
 283             line2string(&line);
 284             strcpy(string, line);
 285         }
 286     }
 287 }
 288
 289 void read_xpm_entry(FILE *in, t_matrix *mm)
 290 {
 291     t_mapping   *map;
 292     char        *line_buf = NULL, *line = NULL, *str, buf[256] = {0};
 293     int          i, m, col_len, nch = 0, n_axis_x, n_axis_y, llmax;
 294     int          llalloc = 0;
 295     unsigned int r, g, b;
 296     double       u;
 297     gmx_bool     bGetOnWithIt, bSetLine;
 298     t_xpmelmt    c;
 299
 300     mm->flags      = 0;
 301     mm->title[0]   = 0;
 302     mm->legend[0]  = 0;
 303     mm->label_x[0] = 0;
 304     mm->label_y[0] = 0;
 305     mm->matrix     = NULL;
 306     mm->axis_x     = NULL;
 307     mm->axis_y     = NULL;
 308     mm->bDiscrete  = FALSE;
 309
 310     llmax = STRLEN;
 311
 312     while ((NULL != fgetline(&line_buf, llmax, &llalloc, in)) &&
 313            (strncmp(line_buf, "static", 6) != 0))
 314     {
 315         line = line_buf;
 316         parsestring(line, "title", (mm->title));
 317         parsestring(line, "legend", (mm->legend));
 318         parsestring(line, "x-label", (mm->label_x));
 319         parsestring(line, "y-label", (mm->label_y));
 320         parsestring(line, "type", buf);
 321     }
 322
 323     if (!line || strncmp(line, "static", 6) != 0)
 324     {
 325         gmx_input("Invalid XPixMap");
 326     }
 327
 328     if (buf[0] && (gmx_strcasecmp(buf, "Discrete") == 0))
 329     {
 330         mm->bDiscrete = TRUE;
 331     }
 332
 333     if (debug)
 334     {
 335         fprintf(debug, "%s %s %s %s\n",
 336                 mm->title, mm->legend, mm->label_x, mm->label_y);
 337     }
 338
 339     /* Read sizes */
 340     bGetOnWithIt = FALSE;
 341     while (!bGetOnWithIt && (NULL != fgetline(&line_buf, llmax, &llalloc, in)))
 342     {
 343         line = line_buf;
 344         while (( line[0] != '\"' ) && ( line[0] != '\0' ))
 345         {
 346             line++;
 347         }
 348
 349         if  (line[0] == '\"')
 350         {
 351             line2string(&line);
 352             sscanf(line, "%d %d %d %d", &(mm->nx), &(mm->ny), &(mm->nmap), &nch);
 353             if (nch > 2)
 354             {
 355                 gmx_fatal(FARGS, "Sorry can only read xpm's with at most 2 caracters per pixel\n");
 356             }
 357             if (mm->nx <= 0 || mm->ny <= 0)
 358             {
 359                 gmx_fatal(FARGS, "Dimensions of xpm-file have to be larger than 0\n");
 360             }
 361             llmax        = std::max(STRLEN, mm->nx+10);
 362             bGetOnWithIt = TRUE;
 363         }
 364     }
 365     if (debug)
 366     {
 367         fprintf(debug, "mm->nx %d mm->ny %d mm->nmap %d nch %d\n",
 368                 mm->nx, mm->ny, mm->nmap, nch);
 369     }
 370     if (nch == 0)
 371     {
 372         gmx_fatal(FARGS, "Number of characters per pixel not found in xpm\n");
 373     }
 374
 375     /* Read color map */
 376     snew(map, mm->nmap);
 377     m = 0;
 378     while ((m < mm->nmap) && (NULL != fgetline(&line_buf, llmax, &llalloc, in)))
 379     {
 380         line = strchr(line_buf, '\"');
 381         if  (line)
 382         {
 383             line++;
 384             /* Read xpm color map entry */
 385             map[m].code.c1 = line[0];
 386             if (nch == 1)
 387             {
 388                 map[m].code.c2 = 0;
 389             }
 390             else
 391             {
 392                 map[m].code.c2 = line[1];
 393             }
 394             line += nch;
 395             str   = strchr(line, '#');
 396             if (str)
 397             {
 398                 str++;
 399                 col_len = 0;
 400                 while (isxdigit(str[col_len]))
 401                 {
 402                     col_len++;
 403                 }
 404                 if (col_len == 6)
 405                 {
 406                     sscanf(line, "%*s #%2x%2x%2x", &r, &g, &b);
 407                     map[m].rgb.r = r/255.0;
 408                     map[m].rgb.g = g/255.0;
 409                     map[m].rgb.b = b/255.0;
 410                 }
 411                 else if (col_len == 12)
 412                 {
 413                     sscanf(line, "%*s #%4x%4x%4x", &r, &g, &b);
 414                     map[m].rgb.r = r/65535.0;
 415                     map[m].rgb.g = g/65535.0;
 416                     map[m].rgb.b = b/65535.0;
 417                 }
 418                 else
 419                 {
 420                     gmx_file("Unsupported or invalid colormap in X PixMap");
 421                 }
 422             }
 423             else
 424             {
 425                 str = strchr(line, 'c');
 426                 if (str)
 427                 {
 428                     str += 2;
 429                 }
 430                 else
 431                 {
 432                     gmx_file("Unsupported or invalid colormap in X PixMap");
 433                 }
 434                 fprintf(stderr, "Using white for color \"%s", str);
 435                 map[m].rgb.r = 1;
 436                 map[m].rgb.g = 1;
 437                 map[m].rgb.b = 1;
 438             }
 439             line = strchr(line, '\"');
 440             line++;
 441             line2string(&line);
 442             map[m].desc = strdup(line);
 443             m++;
 444         }
 445     }
 446     if  (m != mm->nmap)
 447     {
 448         gmx_fatal(FARGS, "Number of read colors map entries (%d) does not match the number in the header (%d)", m, mm->nmap);
 449     }
 450     mm->map = map;
 451
 452     /* Read axes, if there are any */
 453     n_axis_x     = 0;
 454     n_axis_y     = 0;
 455     bSetLine     = FALSE;
 456     do
 457     {
 458         if (bSetLine)
 459         {
 460             line = line_buf;
 461         }
 462         bSetLine = TRUE;
 463         if (strstr(line, "x-axis"))
 464         {
 465             line = strstr(line, "x-axis");
 466             skipstr(line);
 467             if (mm->axis_x == NULL)
 468             {
 469                 snew(mm->axis_x, mm->nx + 1);
 470             }
 471             while (sscanf(line, "%lf", &u) == 1)
 472             {
 473                 if (n_axis_x > mm->nx)
 474                 {
 475                     gmx_fatal(FARGS, "Too many x-axis labels in xpm (max %d)", mm->nx);
 476                 }
 477                 else if (n_axis_x == mm->nx)
 478                 {
 479                     mm->flags |= MAT_SPATIAL_X;
 480                 }
 481                 mm->axis_x[n_axis_x] = u;
 482                 n_axis_x++;
 483                 skipstr(line);
 484             }
 485         }
 486         else if (strstr(line, "y-axis"))
 487         {
 488             line = strstr(line, "y-axis");
 489             skipstr(line);
 490             if (mm->axis_y == NULL)
 491             {
 492                 snew(mm->axis_y, mm->ny + 1);
 493             }
 494             while (sscanf(line, "%lf", &u) == 1)
 495             {
 496                 if (n_axis_y > mm->ny)
 497                 {
 498                     gmx_file("Too many y-axis labels in xpm");
 499                 }
 500                 else if (n_axis_y == mm->ny)
 501                 {
 502                     mm->flags |= MAT_SPATIAL_Y;
 503                 }
 504                 mm->axis_y[n_axis_y] = u;
 505                 n_axis_y++;
 506                 skipstr(line);
 507             }
 508         }
 509     }
 510     while ((line[0] != '\"') && (NULL != fgetline(&line_buf, llmax, &llalloc, in)));
 511
 512     /* Read matrix */
 513     snew(mm->matrix, mm->nx);
 514     for (i = 0; i < mm->nx; i++)
 515     {
 516         snew(mm->matrix[i], mm->ny);
 517     }
 518     m        = mm->ny-1;
 519     bSetLine = FALSE;
 520     do
 521     {
 522         if (bSetLine)
 523         {
 524             line = line_buf;
 525         }
 526         bSetLine = TRUE;
 527         if (m%(1+mm->ny/100) == 0)
 528         {
 529             fprintf(stderr, "%3d%%\b\b\b\b", (100*(mm->ny-m))/mm->ny);
 530         }
 531         while ((line[0] != '\"') && (line[0] != '\0'))
 532         {
 533             line++;
 534         }
 535         if (line[0] != '\"')
 536         {
 537             gmx_fatal(FARGS, "Not enough caracters in row %d of the matrix\n", m+1);
 538         }
 539         else
 540         {
 541             line++;
 542             for (i = 0; i < mm->nx; i++)
 543             {
 544                 c.c1 = line[nch*i];
 545                 if (nch == 1)
 546                 {
 547                     c.c2 = 0;
 548                 }
 549                 else
 550                 {
 551                     c.c2 = line[nch*i+1];
 552                 }
 553                 mm->matrix[i][m] = searchcmap(mm->nmap, mm->map, c);
 554             }
 555             m--;
 556         }
 557     }
 558     while ((m >= 0) && (NULL != fgetline(&line_buf, llmax, &llalloc, in)));
 559     if (m >= 0)
 560     {
 561         gmx_incons("Not enough rows in the matrix");
 562     }
 563
 564     sfree(line_buf);
 565 }
 566
 567 int read_xpm_matrix(const char *fnm, t_matrix **matrix)
 568 {
 569     FILE *in;
 570     char *line = NULL;
 571     int   nmat;
 572     int   llalloc = 0;
 573
 574     in = gmx_fio_fopen(fnm, "r");
 575
 576     nmat = 0;
 577     while (NULL != fgetline(&line, STRLEN, &llalloc, in))
 578     {
 579         if (strstr(line, "/* XPM */"))
 580         {
 581             srenew(*matrix, nmat+1);
 582             read_xpm_entry(in, &(*matrix)[nmat]);
 583             nmat++;
 584         }
 585     }
 586     gmx_fio_fclose(in);
 587
 588     if (nmat == 0)
 589     {
 590         gmx_file("Invalid XPixMap");
 591     }
 592
 593     sfree(line);
 594
 595     return nmat;
 596 }
 597
 598 real **matrix2real(t_matrix *matrix, real **mat)
 599 {
 600     t_mapping *map;
 601     double     tmp;
 602     real      *rmap;
 603     int        i, j, nmap;
 604
 605     nmap = matrix->nmap;
 606     map  = matrix->map;
 607     snew(rmap, nmap);
 608
 609     for (i = 0; i < nmap; i++)
 610     {
 611         if ((map[i].desc == NULL) || (sscanf(map[i].desc, "%lf", &tmp) != 1))
 612         {
 613             fprintf(stderr, "Could not convert matrix to reals,\n"
 614                     "color map entry %d has a non-real description: \"%s\"\n",
 615                     i, map[i].desc);
 616             sfree(rmap);
 617             return NULL;
 618         }
 619         rmap[i] = tmp;
 620     }
 621
 622     if (mat == NULL)
 623     {
 624         snew(mat, matrix->nx);
 625         for (i = 0; i < matrix->nx; i++)
 626         {
 627             snew(mat[i], matrix->ny);
 628         }
 629     }
 630     for (i = 0; i < matrix->nx; i++)
 631     {
 632         for (j = 0; j < matrix->ny; j++)
 633         {
 634             mat[i][j] = rmap[matrix->matrix[i][j]];
 635         }
 636     }
 637
 638     sfree(rmap);
 639
 640     fprintf(stderr, "Converted a %dx%d matrix with %d levels to reals\n",
 641             matrix->nx, matrix->ny, nmap);
 642
 643     return mat;
 644 }
 645
 646 void write_xpm_header(FILE *out,
 647                       const char *title, const char *legend,
 648                       const char *label_x, const char *label_y,
 649                       gmx_bool bDiscrete)
 650 {
 651     fprintf(out,  "/* XPM */\n");
 652     try
 653     {
 654         gmx::BinaryInformationSettings settings;
 655         settings.generatedByHeader(true);
 656         settings.linePrefix("/* ");
 657         settings.lineSuffix(" */");
 658         gmx::printBinaryInformation(out, gmx::getProgramContext(),
 659                                     settings);
 660     }
 661     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 662     fprintf(out,  "/* This file can be converted to EPS by the GROMACS program xpm2ps */\n");
 663     fprintf(out,  "/* title:   \"%s\" */\n", title);
 664     fprintf(out,  "/* legend:  \"%s\" */\n", legend);
 665     fprintf(out,  "/* x-label: \"%s\" */\n", label_x);
 666     fprintf(out,  "/* y-label: \"%s\" */\n", label_y);
 667     if (bDiscrete)
 668     {
 669         fprintf(out, "/* type:    \"Discrete\" */\n");
 670     }
 671     else
 672     {
 673         fprintf(out, "/* type:    \"Continuous\" */\n");
 674     }
 675 }
 676
 677 static int calc_nmid(int nlevels, real lo, real mid, real hi)
 678 {
 679     /* Take care that we have at least 1 entry in the mid to hi range
 680      */
 681     return std::min(std::max(0, static_cast<int>(((mid-lo)/(hi-lo))*(nlevels-1))),
 682                     nlevels-1);
 683 }
 684
 685 void write_xpm_map3(FILE *out, int n_x, int n_y, int *nlevels,
 686                     real lo, real mid, real hi,
 687                     t_rgb rlo, t_rgb rmid, t_rgb rhi)
 688 {
 689     int    i, nmid;
 690     real   r, g, b, clev_lo, clev_hi;
 691
 692     if (*nlevels > NMAP*NMAP)
 693     {
 694         fprintf(stderr, "Warning, too many levels (%d) in matrix, using %d only\n",
 695                 *nlevels, (int)(NMAP*NMAP));
 696         *nlevels = NMAP*NMAP;
 697     }
 698     else if (*nlevels < 2)
 699     {
 700         fprintf(stderr, "Warning, too few levels (%d) in matrix, using 2 instead\n",
 701                 *nlevels);
 702         *nlevels = 2;
 703     }
 704     if (!((mid >= lo) && (mid < hi)))
 705     {
 706         gmx_fatal(FARGS, "Lo: %f, Mid: %f, Hi: %f\n", lo, mid, hi);
 707     }
 708
 709     fprintf(out, "static char *gromacs_xpm[] = {\n");
 710     fprintf(out, "\"%d %d   %d %d\",\n",
 711             n_x, n_y, *nlevels, (*nlevels <= NMAP) ? 1 : 2);
 712
 713     nmid    = calc_nmid(*nlevels, lo, mid, hi);
 714     clev_lo = nmid;
 715     clev_hi = (*nlevels - 1 - nmid);
 716     for (i = 0; (i < nmid); i++)
 717     {
 718         r   = rlo.r+(i*(rmid.r-rlo.r)/clev_lo);
 719         g   = rlo.g+(i*(rmid.g-rlo.g)/clev_lo);
 720         b   = rlo.b+(i*(rmid.b-rlo.b)/clev_lo);
 721         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 722                 mapper[i % NMAP],
 723                 (*nlevels <= NMAP) ? ' ' : mapper[i/NMAP],
 724                 (unsigned int)round(255*r),
 725                 (unsigned int)round(255*g),
 726                 (unsigned int)round(255*b),
 727                 ((nmid - i)*lo + i*mid)/clev_lo);
 728     }
 729     for (i = 0; (i < (*nlevels-nmid)); i++)
 730     {
 731         r   = rmid.r+(i*(rhi.r-rmid.r)/clev_hi);
 732         g   = rmid.g+(i*(rhi.g-rmid.g)/clev_hi);
 733         b   = rmid.b+(i*(rhi.b-rmid.b)/clev_hi);
 734         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 735                 mapper[(i+nmid) % NMAP],
 736                 (*nlevels <= NMAP) ? ' ' : mapper[(i+nmid)/NMAP],
 737                 (unsigned int)round(255*r),
 738                 (unsigned int)round(255*g),
 739                 (unsigned int)round(255*b),
 740                 ((*nlevels - 1 - nmid - i)*mid + i*hi)/clev_hi);
 741     }
 742 }
 743
 744 static void pr_simple_cmap(FILE *out, real lo, real hi, int nlevel, t_rgb rlo,
 745                            t_rgb rhi, int i0)
 746 {
 747     int    i;
 748     real   r, g, b, fac;
 749
 750     for (i = 0; (i < nlevel); i++)
 751     {
 752         fac = (i+1.0)/(nlevel);
 753         r   = rlo.r+fac*(rhi.r-rlo.r);
 754         g   = rlo.g+fac*(rhi.g-rlo.g);
 755         b   = rlo.b+fac*(rhi.b-rlo.b);
 756         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 757                 mapper[(i+i0) % NMAP],
 758                 (nlevel <= NMAP) ? ' ' : mapper[(i+i0)/NMAP],
 759                 (unsigned int)round(255*r),
 760                 (unsigned int)round(255*g),
 761                 (unsigned int)round(255*b),
 762                 lo+fac*(hi-lo));
 763     }
 764 }
 765
 766 static void pr_discrete_cmap(FILE *out, int *nlevel, int i0)
 767 {
 768     t_rgb  rgbd[16] = {
 769         { 1.0, 1.0, 1.0 }, /* white */
 770         { 1.0, 0.0, 0.0 }, /* red */
 771         { 1.0, 1.0, 0.0 }, /* yellow */
 772         { 0.0, 0.0, 1.0 }, /* blue */
 773         { 0.0, 1.0, 0.0 }, /* green */
 774         { 1.0, 0.0, 1.0 }, /* purple */
 775         { 1.0, 0.4, 0.0 }, /* orange */
 776         { 0.0, 1.0, 1.0 }, /* cyan */
 777         { 1.0, 0.4, 0.4 }, /* pink */
 778         { 1.0, 1.0, 0.0 }, /* yellow */
 779         { 0.4, 0.4, 1.0 }, /* lightblue */
 780         { 0.4, 1.0, 0.4 }, /* lightgreen */
 781         { 1.0, 0.4, 1.0 }, /* lightpurple */
 782         { 1.0, 0.7, 0.4 }, /* lightorange */
 783         { 0.4, 1.0, 1.0 }, /* lightcyan */
 784         { 0.0, 0.0, 0.0 }  /* black */
 785     };
 786
 787     int    i, n;
 788
 789     *nlevel = std::min(16, *nlevel);
 790     n       = *nlevel;
 791     for (i = 0; (i < n); i++)
 792     {
 793         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%3d\" */,\n",
 794                 mapper[(i+i0) % NMAP],
 795                 (n <= NMAP) ? ' ' : mapper[(i+i0)/NMAP],
 796                 (unsigned int)round(255*rgbd[i].r),
 797                 (unsigned int)round(255*rgbd[i].g),
 798                 (unsigned int)round(255*rgbd[i].b),
 799                 i);
 800     }
 801 }
 802
 803
 804
 805 void write_xpm_map_split(FILE *out, int n_x, int n_y,
 806                          int *nlevel_top, real lo_top, real hi_top,
 807                          t_rgb rlo_top, t_rgb rhi_top,
 808                          gmx_bool bDiscreteColor,
 809                          int *nlevel_bot, real lo_bot, real hi_bot,
 810                          t_rgb rlo_bot, t_rgb rhi_bot)
 811 {
 812     int    ntot;
 813
 814     ntot = *nlevel_top + *nlevel_bot;
 815     if (ntot > NMAP)
 816     {
 817         gmx_fatal(FARGS, "Warning, too many levels (%d) in matrix", ntot);
 818     }
 819
 820     fprintf(out, "static char *gromacs_xpm[] = {\n");
 821     fprintf(out, "\"%d %d   %d %d\",\n", n_x, n_y, ntot, 1);
 822
 823     if (bDiscreteColor)
 824     {
 825         pr_discrete_cmap(out, nlevel_bot, 0);
 826     }
 827     else
 828     {
 829         pr_simple_cmap(out, lo_bot, hi_bot, *nlevel_bot, rlo_bot, rhi_bot, 0);
 830     }
 831
 832     pr_simple_cmap(out, lo_top, hi_top, *nlevel_top, rlo_top, rhi_top, *nlevel_bot);
 833 }
 834
 835
 836 void write_xpm_map(FILE *out, int n_x, int n_y, int *nlevels, real lo, real hi,
 837                    t_rgb rlo, t_rgb rhi)
 838 {
 839     int    i, nlo;
 840     real   invlevel, r, g, b;
 841
 842     if (*nlevels > NMAP*NMAP)
 843     {
 844         fprintf(stderr, "Warning, too many levels (%d) in matrix, using %d only\n",
 845                 *nlevels, (int)(NMAP*NMAP));
 846         *nlevels = NMAP*NMAP;
 847     }
 848     else if (*nlevels < 2)
 849     {
 850         fprintf(stderr, "Warning, too few levels (%d) in matrix, using 2 instead\n", *nlevels);
 851         *nlevels = 2;
 852     }
 853
 854     fprintf(out, "static char *gromacs_xpm[] = {\n");
 855     fprintf(out, "\"%d %d   %d %d\",\n",
 856             n_x, n_y, *nlevels, (*nlevels <= NMAP) ? 1 : 2);
 857
 858     invlevel = 1.0/(*nlevels-1);
 859     for (i = 0; (i < *nlevels); i++)
 860     {
 861         nlo = *nlevels-1-i;
 862         r   = (nlo*rlo.r+i*rhi.r)*invlevel;
 863         g   = (nlo*rlo.g+i*rhi.g)*invlevel;
 864         b   = (nlo*rlo.b+i*rhi.b)*invlevel;
 865         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 866                 mapper[i % NMAP], (*nlevels <= NMAP) ? ' ' : mapper[i/NMAP],
 867                 (unsigned int)round(255*r),
 868                 (unsigned int)round(255*g),
 869                 (unsigned int)round(255*b),
 870                 (nlo*lo+i*hi)*invlevel);
 871     }
 872 }
 873
 874 void write_xpm_axis(FILE *out, const char *axis, gmx_bool bSpatial, int n,
 875                     real *label)
 876 {
 877     int i;
 878
 879     if (label)
 880     {
 881         for (i = 0; i < (bSpatial ? n+1 : n); i++)
 882         {
 883             if (i % 80 == 0)
 884             {
 885                 if (i)
 886                 {
 887                     fprintf(out, "*/\n");
 888                 }
 889                 fprintf(out, "/* %s-axis:  ", axis);
 890             }
 891             fprintf(out, "%g ", label[i]);
 892         }
 893         fprintf(out, "*/\n");
 894     }
 895 }
 896
 897 void write_xpm_data(FILE *out, int n_x, int n_y, real **matrix,
 898                     real lo, real hi, int nlevels)
 899 {
 900     int  i, j, c;
 901     real invlevel;
 902
 903     invlevel = (nlevels-1)/(hi-lo);
 904     for (j = n_y-1; (j >= 0); j--)
 905     {
 906         if (j%(1+n_y/100) == 0)
 907         {
 908             fprintf(stderr, "%3d%%\b\b\b\b", (100*(n_y-j))/n_y);
 909         }
 910         fprintf(out, "\"");
 911         for (i = 0; (i < n_x); i++)
 912         {
 913             c = gmx_nint((matrix[i][j]-lo)*invlevel);
 914             if (c < 0)
 915             {
 916                 c = 0;
 917             }
 918             if (c >= nlevels)
 919             {
 920                 c = nlevels-1;
 921             }
 922             if (nlevels <= NMAP)
 923             {
 924                 fprintf(out, "%c", mapper[c]);
 925             }
 926             else
 927             {
 928                 fprintf(out, "%c%c", mapper[c % NMAP], mapper[c / NMAP]);
 929             }
 930         }
 931         if (j > 0)
 932         {
 933             fprintf(out, "\",\n");
 934         }
 935         else
 936         {
 937             fprintf(out, "\"\n");
 938         }
 939     }
 940 }
 941
 942 void write_xpm_data3(FILE *out, int n_x, int n_y, real **matrix,
 943                      real lo, real mid, real hi, int nlevels)
 944 {
 945     int  i, j, c = 0, nmid;
 946     real invlev_lo, invlev_hi;
 947
 948     nmid      = calc_nmid(nlevels, lo, mid, hi);
 949     invlev_hi = (nlevels-1-nmid)/(hi-mid);
 950     invlev_lo = (nmid)/(mid-lo);
 951
 952     for (j = n_y-1; (j >= 0); j--)
 953     {
 954         if (j%(1+n_y/100) == 0)
 955         {
 956             fprintf(stderr, "%3d%%\b\b\b\b", (100*(n_y-j))/n_y);
 957         }
 958         fprintf(out, "\"");
 959         for (i = 0; (i < n_x); i++)
 960         {
 961             if (matrix[i][j] >= mid)
 962             {
 963                 c = nmid+gmx_nint((matrix[i][j]-mid)*invlev_hi);
 964             }
 965             else if (matrix[i][j] >= lo)
 966             {
 967                 c = gmx_nint((matrix[i][j]-lo)*invlev_lo);
 968             }
 969             else
 970             {
 971                 c = 0;
 972             }
 973
 974             if (c < 0)
 975             {
 976                 c = 0;
 977             }
 978             if (c >= nlevels)
 979             {
 980                 c = nlevels-1;
 981             }
 982             if (nlevels <= NMAP)
 983             {
 984                 fprintf(out, "%c", mapper[c]);
 985             }
 986             else
 987             {
 988                 fprintf(out, "%c%c", mapper[c % NMAP], mapper[c / NMAP]);
 989             }
 990         }
 991         if (j > 0)
 992         {
 993             fprintf(out, "\",\n");
 994         }
 995         else
 996         {
 997             fprintf(out, "\"\n");
 998         }
 999     }
1000 }
1001
1002 void write_xpm_data_split(FILE *out, int n_x, int n_y, real **matrix,
1003                           real lo_top, real hi_top, int nlevel_top,
1004                           real lo_bot, real hi_bot, int nlevel_bot)
1005 {
1006     int  i, j, c;
1007     real invlev_top, invlev_bot;
1008
1009     invlev_top = (nlevel_top-1)/(hi_top-lo_top);
1010     invlev_bot = (nlevel_bot-1)/(hi_bot-lo_bot);
1011
1012     for (j = n_y-1; (j >= 0); j--)
1013     {
1014         if (j % (1+n_y/100) == 0)
1015         {
1016             fprintf(stderr, "%3d%%\b\b\b\b", (100*(n_y-j))/n_y);
1017         }
1018         fprintf(out, "\"");
1019         for (i = 0; (i < n_x); i++)
1020         {
1021             if (i < j)
1022             {
1023                 c = nlevel_bot+round((matrix[i][j]-lo_top)*invlev_top);
1024                 if ((c < nlevel_bot) || (c >= nlevel_bot+nlevel_top))
1025                 {
1026                     gmx_fatal(FARGS, "Range checking i = %d, j = %d, c = %d, bot = %d, top = %d matrix[i,j] = %f", i, j, c, nlevel_bot, nlevel_top, matrix[i][j]);
1027                 }
1028             }
1029             else if (i > j)
1030             {
1031                 c = round((matrix[i][j]-lo_bot)*invlev_bot);
1032                 if ((c < 0) || (c >= nlevel_bot+nlevel_bot))
1033                 {
1034                     gmx_fatal(FARGS, "Range checking i = %d, j = %d, c = %d, bot = %d, top = %d matrix[i,j] = %f", i, j, c, nlevel_bot, nlevel_top, matrix[i][j]);
1035                 }
1036             }
1037             else
1038             {
1039                 c = nlevel_bot;
1040             }
1041
1042             fprintf(out, "%c", mapper[c]);
1043         }
1044         if (j > 0)
1045         {
1046             fprintf(out, "\",\n");
1047         }
1048         else
1049         {
1050             fprintf(out, "\"\n");
1051         }
1052     }
1053 }
1054
1055 void write_xpm_m(FILE *out, t_matrix m)
1056 {
1057     /* Writes a t_matrix struct to .xpm file */
1058
1059     int           i, j;
1060     gmx_bool      bOneChar;
1061     t_xpmelmt     c;
1062
1063     bOneChar = (m.map[0].code.c2 == 0);
1064     write_xpm_header(out, m.title, m.legend, m.label_x, m.label_y,
1065                      m.bDiscrete);
1066     fprintf(out, "static char *gromacs_xpm[] = {\n");
1067     fprintf(out, "\"%d %d   %d %d\",\n", m.nx, m.ny, m.nmap, bOneChar ? 1 : 2);
1068     for (i = 0; (i < m.nmap); i++)
1069     {
1070         fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%s\" */,\n",
1071                 m.map[i].code.c1,
1072                 bOneChar ? ' ' : m.map[i].code.c2,
1073                 (unsigned int)round(m.map[i].rgb.r*255),
1074                 (unsigned int)round(m.map[i].rgb.g*255),
1075                 (unsigned int)round(m.map[i].rgb.b*255), m.map[i].desc);
1076     }
1077     write_xpm_axis(out, "x", m.flags & MAT_SPATIAL_X, m.nx, m.axis_x);
1078     write_xpm_axis(out, "y", m.flags & MAT_SPATIAL_Y, m.ny, m.axis_y);
1079     for (j = m.ny-1; (j >= 0); j--)
1080     {
1081         if (j%(1+m.ny/100) == 0)
1082         {
1083             fprintf(stderr, "%3d%%\b\b\b\b", (100*(m.ny-j))/m.ny);
1084         }
1085         fprintf(out, "\"");
1086         if (bOneChar)
1087         {
1088             for (i = 0; (i < m.nx); i++)
1089             {
1090                 fprintf(out, "%c", m.map[m.matrix[i][j]].code.c1);
1091             }
1092         }
1093         else
1094         {
1095             for (i = 0; (i < m.nx); i++)
1096             {
1097                 c = m.map[m.matrix[i][j]].code;
1098                 fprintf(out, "%c%c", c.c1, c.c2);
1099             }
1100         }
1101         if (j > 0)
1102         {
1103             fprintf(out, "\",\n");
1104         }
1105         else
1106         {
1107             fprintf(out, "\"\n");
1108         }
1109     }
1110 }
1111
1112 void write_xpm3(FILE *out, unsigned int flags,
1113                 const char *title, const char *legend,
1114                 const char *label_x, const char *label_y,
1115                 int n_x, int n_y, real axis_x[], real axis_y[],
1116                 real *matrix[], real lo, real mid, real hi,
1117                 t_rgb rlo, t_rgb rmid, t_rgb rhi, int *nlevels)
1118 {
1119     /* See write_xpm.
1120      * Writes a colormap varying as rlo -> rmid -> rhi.
1121      */
1122
1123     if (hi <= lo)
1124     {
1125         gmx_fatal(FARGS, "hi (%g) <= lo (%g)", hi, lo);
1126     }
1127
1128     write_xpm_header(out, title, legend, label_x, label_y, FALSE);
1129     write_xpm_map3(out, n_x, n_y, nlevels, lo, mid, hi, rlo, rmid, rhi);
1130     write_xpm_axis(out, "x", flags & MAT_SPATIAL_X, n_x, axis_x);
1131     write_xpm_axis(out, "y", flags & MAT_SPATIAL_Y, n_y, axis_y);
1132     write_xpm_data3(out, n_x, n_y, matrix, lo, mid, hi, *nlevels);
1133 }
1134
1135 void write_xpm_split(FILE *out, unsigned int flags,
1136                      const char *title, const char *legend,
1137                      const char *label_x, const char *label_y,
1138                      int n_x, int n_y, real axis_x[], real axis_y[],
1139                      real *matrix[],
1140                      real lo_top, real hi_top, int *nlevel_top,
1141                      t_rgb rlo_top, t_rgb rhi_top,
1142                      real lo_bot, real hi_bot, int *nlevel_bot,
1143                      gmx_bool bDiscreteColor,
1144                      t_rgb rlo_bot, t_rgb rhi_bot)
1145 {
1146     /* See write_xpm.
1147      * Writes a colormap varying as rlo -> rmid -> rhi.
1148      */
1149
1150     if (hi_top <= lo_top)
1151     {
1152         gmx_fatal(FARGS, "hi_top (%g) <= lo_top (%g)", hi_top, lo_top);
1153     }
1154     if (hi_bot <= lo_bot)
1155     {
1156         gmx_fatal(FARGS, "hi_bot (%g) <= lo_bot (%g)", hi_bot, lo_bot);
1157     }
1158     if (bDiscreteColor && (*nlevel_bot >= 16))
1159     {
1160         gmx_impl("Can not plot more than 16 discrete colors");
1161     }
1162
1163     write_xpm_header(out, title, legend, label_x, label_y, FALSE);
1164     write_xpm_map_split(out, n_x, n_y, nlevel_top, lo_top, hi_top, rlo_top, rhi_top,
1165                         bDiscreteColor, nlevel_bot, lo_bot, hi_bot, rlo_bot, rhi_bot);
1166     write_xpm_axis(out, "x", flags & MAT_SPATIAL_X, n_x, axis_x);
1167     write_xpm_axis(out, "y", flags & MAT_SPATIAL_Y, n_y, axis_y);
1168     write_xpm_data_split(out, n_x, n_y, matrix, lo_top, hi_top, *nlevel_top,
1169                          lo_bot, hi_bot, *nlevel_bot);
1170 }
1171
1172 void write_xpm(FILE *out, unsigned int flags,
1173                const char *title, const char *legend,
1174                const char *label_x, const char *label_y,
1175                int n_x, int n_y, real axis_x[], real axis_y[],
1176                real *matrix[], real lo, real hi,
1177                t_rgb rlo, t_rgb rhi, int *nlevels)
1178 {
1179     /* out        xpm file
1180      * title      matrix title
1181      * legend     label for the continuous legend
1182      * label_x    label for the x-axis
1183      * label_y    label for the y-axis
1184      * n_x, n_y   size of the matrix
1185      * axis_x[]   the x-ticklabels
1186      * axis_y[]   the y-ticklables
1187      * *matrix[]  element x,y is matrix[x][y]
1188      * lo         output lower than lo is set to lo
1189      * hi         output higher than hi is set to hi
1190      * rlo        rgb value for level lo
1191      * rhi        rgb value for level hi
1192      * nlevels    number of color levels for the output
1193      */
1194
1195     if (hi <= lo)
1196     {
1197         gmx_fatal(FARGS, "hi (%f) <= lo (%f)", hi, lo);
1198     }
1199
1200     write_xpm_header(out, title, legend, label_x, label_y, FALSE);
1201     write_xpm_map(out, n_x, n_y, nlevels, lo, hi, rlo, rhi);
1202     write_xpm_axis(out, "x", flags & MAT_SPATIAL_X, n_x, axis_x);
1203     write_xpm_axis(out, "y", flags & MAT_SPATIAL_Y, n_y, axis_y);
1204     write_xpm_data(out, n_x, n_y, matrix, lo, hi, *nlevels);
1205 }