src/gromacs/gmxlib/matio.c

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