src/gromacs/gmxlib/matio.cpp

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