src/gromacs/gmxana/gmx_xpm2ps.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <cmath>
  40 #include <cstdio>
  41 #include <cstdlib>
  42 #include <cstring>
  43
  44 #include <algorithm>
  45
  46 #include "gromacs/commandline/pargs.h"
  47 #include "gromacs/fileio/matio.h"
  48 #include "gromacs/fileio/trxio.h"
  49 #include "gromacs/fileio/writeps.h"
  50 #include "gromacs/gmxana/gmx_ana.h"
  51 #include "gromacs/legacyheaders/macros.h"
  52 #include "gromacs/legacyheaders/typedefs.h"
  53 #include "gromacs/legacyheaders/viewit.h"
  54 #include "gromacs/utility/cstringutil.h"
  55 #include "gromacs/utility/fatalerror.h"
  56 #include "gromacs/utility/futil.h"
  57 #include "gromacs/utility/gmxassert.h"
  58 #include "gromacs/utility/smalloc.h"
  59
  60 #define FUDGE 1.2
  61 #define DDD   2
  62
  63 typedef struct {
  64     real     major;
  65     real     minor;
  66     real     offset;
  67     gmx_bool first;
  68     int      lineatzero;
  69     real     majorticklen;
  70     real     minorticklen;
  71     char     label[STRLEN];
  72     real     fontsize;
  73     char     font[STRLEN];
  74     real     tickfontsize;
  75     char     tickfont[STRLEN];
  76 } t_axisdef;
  77
  78 typedef struct {
  79     int       bw;
  80     real      linewidth;
  81     real      xoffs, yoffs;
  82     gmx_bool  bTitle;
  83     gmx_bool  bTitleOnce;
  84     gmx_bool  bYonce;
  85     real      titfontsize;
  86     char      titfont[STRLEN];
  87     gmx_bool  legend;
  88     real      legfontsize;
  89     char      legfont[STRLEN];
  90     char      leglabel[STRLEN];
  91     char      leg2label[STRLEN];
  92     real      xboxsize;
  93     real      yboxsize;
  94     real      boxspacing;
  95     real      boxlinewidth;
  96     real      ticklinewidth;
  97     real      zerolinewidth;
  98     t_axisdef X, Y;
  99 } t_psrec;
 100
 101 /* MUST correspond to char *legend[] in main() */
 102 enum {
 103     elSel, elBoth, elFirst, elSecond, elNone, elNR
 104 };
 105
 106 /* MUST correspond to char *combine[] in main() */
 107 enum {
 108     ecSel, ecHalves, ecAdd, ecSub, ecMult, ecDiv, ecNR
 109 };
 110
 111 void get_params(const char *mpin, const char *mpout, t_psrec *psr)
 112 {
 113     static const char *gmx_bools[BOOL_NR+1]  = { "no", "yes", NULL };
 114     /* this must correspond to t_rgb *linecolors[] below */
 115     static const char *colors[] = { "none", "black", "white", NULL };
 116     warninp_t          wi;
 117     t_inpfile         *inp;
 118     const char        *tmp;
 119     int                ninp = 0;
 120
 121     wi = init_warning(FALSE, 0);
 122
 123     if (mpin != NULL)
 124     {
 125         inp = read_inpfile(mpin, &ninp, wi);
 126     }
 127     else
 128     {
 129         inp = NULL;
 130     }
 131     ETYPE("black&white",    psr->bw,             gmx_bools);
 132     RTYPE("linewidth",      psr->linewidth,      1.0);
 133     STYPE("titlefont",      psr->titfont,        "Helvetica");
 134     RTYPE("titlefontsize",  psr->titfontsize,    20.0);
 135     ETYPE("legend",         psr->legend,         gmx_bools);
 136     STYPE("legendfont",     psr->legfont,        psr->titfont);
 137     STYPE("legendlabel",    psr->leglabel,       "");
 138     STYPE("legend2label",   psr->leg2label,      psr->leglabel);
 139     RTYPE("legendfontsize", psr->legfontsize,    14.0);
 140     RTYPE("xbox",           psr->xboxsize,       0.0);
 141     RTYPE("ybox",           psr->yboxsize,       0.0);
 142     RTYPE("matrixspacing",  psr->boxspacing,     20.0);
 143     RTYPE("xoffset",        psr->xoffs,          0.0);
 144     RTYPE("yoffset",        psr->yoffs,          psr->xoffs);
 145     RTYPE("boxlinewidth",   psr->boxlinewidth,   psr->linewidth);
 146     RTYPE("ticklinewidth",  psr->ticklinewidth,  psr->linewidth);
 147     RTYPE("zerolinewidth",  psr->zerolinewidth,  psr->ticklinewidth);
 148     ETYPE("x-lineat0value", psr->X.lineatzero,   colors);
 149     RTYPE("x-major",        psr->X.major,        NOTSET);
 150     RTYPE("x-minor",        psr->X.minor,        NOTSET);
 151     RTYPE("x-firstmajor",   psr->X.offset,       0.0);
 152     ETYPE("x-majorat0",     psr->X.first,        gmx_bools);
 153     RTYPE("x-majorticklen", psr->X.majorticklen, 8.0);
 154     RTYPE("x-minorticklen", psr->X.minorticklen, 4.0);
 155     STYPE("x-label",        psr->X.label,        "");
 156     RTYPE("x-fontsize",     psr->X.fontsize,     16.0);
 157     STYPE("x-font",         psr->X.font,         psr->titfont);
 158     RTYPE("x-tickfontsize", psr->X.tickfontsize, 10.0);
 159     STYPE("x-tickfont",     psr->X.tickfont,     psr->X.font);
 160     ETYPE("y-lineat0value", psr->Y.lineatzero,   colors);
 161     RTYPE("y-major",        psr->Y.major,        psr->X.major);
 162     RTYPE("y-minor",        psr->Y.minor,        psr->X.minor);
 163     RTYPE("y-firstmajor",   psr->Y.offset,       psr->X.offset);
 164     ETYPE("y-majorat0",     psr->Y.first,        gmx_bools);
 165     RTYPE("y-majorticklen", psr->Y.majorticklen, psr->X.majorticklen);
 166     RTYPE("y-minorticklen", psr->Y.minorticklen, psr->X.minorticklen);
 167     STYPE("y-label",        psr->Y.label,        psr->X.label);
 168     RTYPE("y-fontsize",     psr->Y.fontsize,     psr->X.fontsize);
 169     STYPE("y-font",         psr->Y.font,         psr->X.font);
 170     RTYPE("y-tickfontsize", psr->Y.tickfontsize, psr->X.tickfontsize);
 171     STYPE("y-tickfont",     psr->Y.tickfont,     psr->Y.font);
 172
 173     if (mpout != NULL)
 174     {
 175         write_inpfile(mpout, ninp, inp, TRUE, wi);
 176     }
 177
 178     done_warning(wi, FARGS);
 179 }
 180
 181 t_rgb black = { 0, 0, 0 };
 182 t_rgb white = { 1, 1, 1 };
 183 t_rgb red   = { 1, 0, 0 };
 184 t_rgb blue  = { 0, 0, 1 };
 185 #define BLACK (&black)
 186 /* this must correspond to *colors[] in get_params */
 187 t_rgb *linecolors[] = { NULL, &black, &white, NULL };
 188
 189 gmx_bool diff_maps(int nmap1, t_mapping *map1, int nmap2, t_mapping *map2)
 190 {
 191     int      i;
 192     gmx_bool bDiff, bColDiff = FALSE;
 193
 194     if (nmap1 != nmap2)
 195     {
 196         bDiff = TRUE;
 197     }
 198     else
 199     {
 200         bDiff = FALSE;
 201         for (i = 0; i < nmap1; i++)
 202         {
 203             if (!matelmt_cmp(map1[i].code, map2[i].code))
 204             {
 205                 bDiff = TRUE;
 206             }
 207             if (std::strcmp(map1[i].desc, map2[i].desc) != 0)
 208             {
 209                 bDiff = TRUE;
 210             }
 211             if ((map1[i].rgb.r != map2[i].rgb.r) ||
 212                 (map1[i].rgb.g != map2[i].rgb.g) ||
 213                 (map1[i].rgb.b != map2[i].rgb.b))
 214             {
 215                 bColDiff = TRUE;
 216             }
 217         }
 218         if (!bDiff && bColDiff)
 219         {
 220             fprintf(stderr, "Warning: two colormaps differ only in RGB value, using one colormap.\n");
 221         }
 222     }
 223
 224     return bDiff;
 225 }
 226
 227 void leg_discrete(t_psdata ps, real x0, real y0, char *label,
 228                   real fontsize, char *font, int nmap, t_mapping map[])
 229 {
 230     int   i;
 231     real  yhh;
 232     real  boxhh;
 233
 234     boxhh = fontsize+DDD;
 235     /* LANDSCAPE */
 236     ps_rgb(ps, BLACK);
 237     ps_strfont(ps, font, fontsize);
 238     yhh = y0+fontsize+3*DDD;
 239     if (std::strlen(label) > 0)
 240     {
 241         ps_ctext(ps, x0, yhh, label, eXLeft);
 242     }
 243     ps_moveto(ps, x0, y0);
 244     for (i = 0; (i < nmap); i++)
 245     {
 246         ps_setorigin(ps);
 247         ps_rgb(ps, &(map[i].rgb));
 248         ps_fillbox(ps, DDD, DDD, DDD+fontsize, boxhh-DDD);
 249         ps_rgb(ps, BLACK);
 250         ps_box(ps, DDD, DDD, DDD+fontsize, boxhh-DDD);
 251         ps_ctext(ps, boxhh+2*DDD, fontsize/3, map[i].desc, eXLeft);
 252         ps_unsetorigin(ps);
 253         ps_moverel(ps, DDD, -fontsize/3);
 254     }
 255 }
 256
 257 void leg_continuous(t_psdata ps, real x0, real x, real y0, char *label,
 258                     real fontsize, char *font,
 259                     int nmap, t_mapping map[],
 260                     int mapoffset)
 261 {
 262     int   i;
 263     real  xx0;
 264     real  yhh, boxxh, boxyh;
 265
 266     boxyh = fontsize;
 267     if (x < 8*fontsize)
 268     {
 269         x = 8*fontsize;
 270     }
 271     boxxh = x/(nmap-mapoffset);
 272     if (boxxh > fontsize)
 273     {
 274         boxxh = fontsize;
 275     }
 276
 277     GMX_RELEASE_ASSERT(map != NULL, "NULL map array provided to leg_continuous()");
 278
 279     /* LANDSCAPE */
 280     xx0 = x0-((nmap-mapoffset)*boxxh)/2.0;
 281
 282     for (i = 0; (i < nmap-mapoffset); i++)
 283     {
 284         ps_rgb(ps, &(map[i+mapoffset].rgb));
 285         ps_fillbox(ps, xx0+i*boxxh, y0, xx0+(i+1)*boxxh, y0+boxyh);
 286     }
 287     ps_strfont(ps, font, fontsize);
 288     ps_rgb(ps, BLACK);
 289     ps_box(ps, xx0, y0, xx0+(nmap-mapoffset)*boxxh, y0+boxyh);
 290
 291     yhh = y0+boxyh+3*DDD;
 292     ps_ctext(ps, xx0+boxxh/2, yhh, map[0].desc, eXCenter);
 293     if (std::strlen(label) > 0)
 294     {
 295         ps_ctext(ps, x0, yhh, label, eXCenter);
 296     }
 297     ps_ctext(ps, xx0+((nmap-mapoffset)*boxxh)
 298              - boxxh/2, yhh, map[nmap-1].desc, eXCenter);
 299 }
 300
 301 void leg_bicontinuous(t_psdata ps, real x0, real x, real y0, char *label1,
 302                       char *label2, real fontsize, char *font,
 303                       int nmap1, t_mapping map1[], int nmap2, t_mapping map2[])
 304 {
 305     real xx1, xx2, x1, x2;
 306
 307     x1  = x/(nmap1+nmap2)*nmap1; /* width of legend 1 */
 308     x2  = x/(nmap1+nmap2)*nmap2; /* width of legend 2 */
 309     xx1 = x0-(x2/2.0)-fontsize;  /* center of legend 1 */
 310     xx2 = x0+(x1/2.0)+fontsize;  /* center of legend 2 */
 311     x1 -= fontsize/2;            /* adjust width */
 312     x2 -= fontsize/2;            /* adjust width */
 313     leg_continuous(ps, xx1, x1, y0, label1, fontsize, font, nmap1, map1, 0);
 314     leg_continuous(ps, xx2, x2, y0, label2, fontsize, font, nmap2, map2, 0);
 315 }
 316
 317 static real box_height(t_matrix *mat, t_psrec *psr)
 318 {
 319     return mat->ny*psr->yboxsize;
 320 }
 321
 322 static real box_dh(t_psrec *psr)
 323 {
 324     return psr->boxspacing;
 325 }
 326
 327 #define IS_ONCE (i == nmat-1)
 328 static real box_dh_top(gmx_bool bOnce, t_psrec *psr)
 329 {
 330     real dh;
 331
 332     if (psr->bTitle || (psr->bTitleOnce && bOnce) )
 333     {
 334         dh = 2*psr->titfontsize;
 335     }
 336     else
 337     {
 338         dh = 0;
 339     }
 340
 341     return dh;
 342 }
 343
 344 static gmx_bool box_do_all_x_maj_ticks(t_psrec *psr)
 345 {
 346     return (psr->boxspacing > (1.5*psr->X.majorticklen));
 347 }
 348
 349 static gmx_bool box_do_all_x_min_ticks(t_psrec *psr)
 350 {
 351     return (psr->boxspacing > (1.5*psr->X.minorticklen));
 352 }
 353
 354 static void draw_boxes(t_psdata ps, real x0, real y0, real w,
 355                        int nmat, t_matrix mat[], t_psrec *psr)
 356 {
 357     char     buf[128];
 358     char    *mylab;
 359     real     xxx;
 360     char   **xtick, **ytick;
 361     real     xx, yy, dy, xx00, yy00, offset_x, offset_y;
 362     int      i, j, x, y, ntx, nty;
 363     size_t   strlength;
 364
 365     /* Only necessary when there will be no y-labels */
 366     strlength = 0;
 367
 368     /* Draw the box */
 369     ps_rgb(ps, BLACK);
 370     ps_linewidth(ps, static_cast<int>(psr->boxlinewidth));
 371     yy00 = y0;
 372     for (i = 0; (i < nmat); i++)
 373     {
 374         dy = box_height(&(mat[i]), psr);
 375         ps_box(ps, x0-1, yy00-1, x0+w+1, yy00+dy+1);
 376         yy00 += dy+box_dh(psr)+box_dh_top(IS_ONCE, psr);
 377     }
 378
 379     /* Draw the ticks on the axes */
 380     ps_linewidth(ps, static_cast<int>(psr->ticklinewidth));
 381     xx00 = x0-1;
 382     yy00 = y0-1;
 383     for (i = 0; (i < nmat); i++)
 384     {
 385         if (mat[i].flags & MAT_SPATIAL_X)
 386         {
 387             ntx      = mat[i].nx + 1;
 388             offset_x = 0.1;
 389         }
 390         else
 391         {
 392             ntx      = mat[i].nx;
 393             offset_x = 0.6;
 394         }
 395         if (mat[i].flags & MAT_SPATIAL_Y)
 396         {
 397             nty      = mat[i].ny + 1;
 398             offset_y = 0.1;
 399         }
 400         else
 401         {
 402             nty      = mat[i].ny;
 403             offset_y = 0.6;
 404         }
 405         snew(xtick, ntx);
 406         for (j = 0; (j < ntx); j++)
 407         {
 408             sprintf(buf, "%g", mat[i].axis_x[j]);
 409             xtick[j] = gmx_strdup(buf);
 410         }
 411         ps_strfont(ps, psr->X.tickfont, psr->X.tickfontsize);
 412         for (x = 0; (x < ntx); x++)
 413         {
 414             xx = xx00 + (x + offset_x)*psr->xboxsize;
 415             if ( ( bRmod(mat[i].axis_x[x], psr->X.offset, psr->X.major) ||
 416                    (psr->X.first && (x == 0))) &&
 417                  ( (i == 0) || box_do_all_x_maj_ticks(psr) ) )
 418             {
 419                 /* Longer tick marks */
 420                 ps_line (ps, xx, yy00, xx, yy00-psr->X.majorticklen);
 421                 /* Plot label on lowest graph only */
 422                 if (i == 0)
 423                 {
 424                     ps_ctext(ps, xx,
 425                              yy00-DDD-psr->X.majorticklen-psr->X.tickfontsize*0.8,
 426                              xtick[x], eXCenter);
 427                 }
 428             }
 429             else if (bRmod(mat[i].axis_x[x], psr->X.offset, psr->X.minor) &&
 430                      ( (i == 0) || box_do_all_x_min_ticks(psr) ) )
 431             {
 432                 /* Shorter tick marks */
 433                 ps_line(ps, xx, yy00, xx, yy00-psr->X.minorticklen);
 434             }
 435             else if (bRmod(mat[i].axis_x[x], psr->X.offset, psr->X.major) )
 436             {
 437                 /* Even shorter marks, only each X.major */
 438                 ps_line(ps, xx, yy00, xx, yy00-(psr->boxspacing/2));
 439             }
 440         }
 441         ps_strfont(ps, psr->Y.tickfont, psr->Y.tickfontsize);
 442         snew(ytick, nty);
 443         for (j = 0; (j < nty); j++)
 444         {
 445             sprintf(buf, "%g", mat[i].axis_y[j]);
 446             ytick[j] = gmx_strdup(buf);
 447         }
 448
 449         for (y = 0; (y < nty); y++)
 450         {
 451             yy = yy00 + (y + offset_y)*psr->yboxsize;
 452             if (bRmod(mat[i].axis_y[y], psr->Y.offset, psr->Y.major) ||
 453                 (psr->Y.first && (y == 0)))
 454             {
 455                 /* Major ticks */
 456                 strlength = std::max(strlength, std::strlen(ytick[y]));
 457                 ps_line (ps, xx00, yy, xx00-psr->Y.majorticklen, yy);
 458                 ps_ctext(ps, xx00-psr->Y.majorticklen-DDD,
 459                          yy-psr->Y.tickfontsize/3.0, ytick[y], eXRight);
 460             }
 461             else if (bRmod(mat[i].axis_y[y], psr->Y.offset, psr->Y.minor) )
 462             {
 463                 /* Minor ticks */
 464                 ps_line(ps, xx00, yy, xx00-psr->Y.minorticklen, yy);
 465             }
 466         }
 467         sfree(xtick);
 468         sfree(ytick);
 469
 470         /* Label on Y-axis */
 471         if (!psr->bYonce || i == nmat/2)
 472         {
 473             if (strlen(psr->Y.label) > 0)
 474             {
 475                 mylab = psr->Y.label;
 476             }
 477             else
 478             {
 479                 mylab = mat[i].label_y;
 480             }
 481             if (strlen(mylab) > 0)
 482             {
 483                 ps_strfont(ps, psr->Y.font, psr->Y.fontsize);
 484                 ps_flip(ps, TRUE);
 485                 xxx = x0-psr->X.majorticklen-psr->X.tickfontsize*strlength-DDD;
 486                 ps_ctext(ps, yy00+box_height(&mat[i], psr)/2.0, 612.5-xxx,
 487                          mylab, eXCenter);
 488                 ps_flip(ps, FALSE);
 489             }
 490         }
 491
 492         yy00 += box_height(&(mat[i]), psr)+box_dh(psr)+box_dh_top(IS_ONCE, psr);
 493     }
 494     /* Label on X-axis */
 495     if (strlen(psr->X.label) > 0)
 496     {
 497         mylab = psr->X.label;
 498     }
 499     else
 500     {
 501         mylab = mat[0].label_x;
 502     }
 503     if (strlen(mylab) > 0)
 504     {
 505         ps_strfont(ps, psr->X.font, psr->X.fontsize);
 506         ps_ctext(ps, x0+w/2, y0-DDD-psr->X.majorticklen-psr->X.tickfontsize*FUDGE-
 507                  psr->X.fontsize, mylab, eXCenter);
 508     }
 509 }
 510
 511 static void draw_zerolines(t_psdata out, real x0, real y0, real w,
 512                            int nmat, t_matrix mat[], t_psrec *psr)
 513 {
 514     real   xx, yy, dy, xx00, yy00;
 515     int    i, x, y;
 516
 517     xx00 = x0-1.5;
 518     yy00 = y0-1.5;
 519     ps_linewidth(out, static_cast<int>(psr->zerolinewidth));
 520     for (i = 0; (i < nmat); i++)
 521     {
 522         dy = box_height(&(mat[i]), psr);
 523         /* mat[i].axis_x and _y were already set by draw_boxes */
 524         if (psr->X.lineatzero)
 525         {
 526             ps_rgb(out, linecolors[psr->X.lineatzero]);
 527             for (x = 0; (x < mat[i].nx); x++)
 528             {
 529                 xx = xx00+(x+0.7)*psr->xboxsize;
 530                 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
 531                 if (x != 0 && x < mat[i].nx-1 &&
 532                     std::abs(mat[i].axis_x[x]) <
 533                     0.1*std::abs(mat[i].axis_x[x+1]-mat[i].axis_x[x]) )
 534                 {
 535                     ps_line (out, xx, yy00, xx, yy00+dy+2);
 536                 }
 537             }
 538         }
 539         if (psr->Y.lineatzero)
 540         {
 541             ps_rgb(out, linecolors[psr->Y.lineatzero]);
 542             for (y = 0; (y < mat[i].ny); y++)
 543             {
 544                 yy = yy00+(y+0.7)*psr->yboxsize;
 545                 /* draw lines whenever tick label almost zero (e.g. next trajectory) */
 546                 if (y != 0 && y < mat[i].ny-1 &&
 547                     std::abs(mat[i].axis_y[y]) <
 548                     0.1*std::abs(mat[i].axis_y[y+1]-mat[i].axis_y[y]) )
 549                 {
 550                     ps_line (out, xx00, yy, xx00+w+2, yy);
 551                 }
 552             }
 553         }
 554         yy00 += box_height(&(mat[i]), psr)+box_dh(psr)+box_dh_top(IS_ONCE, psr);
 555     }
 556 }
 557
 558 static void box_dim(int nmat, t_matrix mat[], t_matrix *mat2, t_psrec *psr,
 559                     int elegend, gmx_bool bFrame,
 560                     real *w, real *h, real *dw, real *dh)
 561 {
 562     int  i, maxytick;
 563     real ww, hh, dww, dhh;
 564
 565     hh       = dww = dhh = 0;
 566     maxytick = 0;
 567
 568     ww = 0;
 569     for (i = 0; (i < nmat); i++)
 570     {
 571         ww       = std::max(ww, mat[i].nx*psr->xboxsize);
 572         hh      += box_height(&(mat[i]), psr);
 573         maxytick = std::max(maxytick, mat[i].nx);
 574     }
 575     if (bFrame)
 576     {
 577         if (mat[0].label_y[0])
 578         {
 579             dww += 2.0*(psr->Y.fontsize+DDD);
 580         }
 581         if (psr->Y.major > 0)
 582         {
 583             dww += psr->Y.majorticklen + DDD +
 584                 psr->Y.tickfontsize*(std::log(static_cast<real>(maxytick))/std::log(10.0));
 585         }
 586         else if (psr->Y.minor > 0)
 587         {
 588             dww += psr->Y.minorticklen;
 589         }
 590
 591         if (mat[0].label_x[0])
 592         {
 593             dhh += psr->X.fontsize+2*DDD;
 594         }
 595         if ( /* fool emacs auto-indent */
 596             (elegend == elBoth && (mat[0].legend[0] || mat2[0].legend[0])) ||
 597             (elegend == elFirst && mat[0].legend[0]) ||
 598             (elegend == elSecond && mat2[0].legend[0]) )
 599         {
 600             dhh += 2*(psr->legfontsize*FUDGE+2*DDD);
 601         }
 602         else
 603         {
 604             dhh += psr->legfontsize*FUDGE+2*DDD;
 605         }
 606         if (psr->X.major > 0)
 607         {
 608             dhh += psr->X.tickfontsize*FUDGE+2*DDD+psr->X.majorticklen;
 609         }
 610         else if (psr->X.minor > 0)
 611         {
 612             dhh += psr->X.minorticklen;
 613         }
 614
 615         hh += (nmat-1)*box_dh(psr);
 616         hh += box_dh_top(TRUE, psr);
 617         if (nmat > 1)
 618         {
 619             hh += (nmat-1)*box_dh_top(FALSE, psr);
 620         }
 621     }
 622     *w  = ww;
 623     *h  = hh;
 624     *dw = dww;
 625     *dh = dhh;
 626 }
 627
 628 int add_maps(t_mapping **newmap,
 629              int nmap1, t_mapping map1[], int nmap2, t_mapping map2[])
 630 {
 631     static char mapper[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<.>/?";
 632     int         nsymbols;
 633     int         nmap, j, k;
 634     t_mapping  *map;
 635
 636     nsymbols = std::strlen(mapper);
 637     nmap     = nmap1+nmap2;
 638     if (nmap > nsymbols*nsymbols)
 639     {
 640         gmx_fatal(FARGS, "Not enough symbols to merge the two colormaps\n");
 641     }
 642     printf("Combining colormaps of %d and %d elements into one of %d elements\n",
 643            nmap1, nmap2, nmap);
 644     snew(map, nmap);
 645     for (j = 0; j < nmap1; j++)
 646     {
 647         map[j].code.c1 = mapper[j % nsymbols];
 648         if (nmap > nsymbols)
 649         {
 650             map[j].code.c2 = mapper[j/nsymbols];
 651         }
 652         map[j].rgb.r = map1[j].rgb.r;
 653         map[j].rgb.g = map1[j].rgb.g;
 654         map[j].rgb.b = map1[j].rgb.b;
 655         map[j].desc  = map1[j].desc;
 656     }
 657     for (j = 0; j < nmap2; j++)
 658     {
 659         k              = j+nmap1;
 660         map[k].code.c1 = mapper[k % nsymbols];
 661         if (nmap > nsymbols)
 662         {
 663             map[k].code.c2 = mapper[k/nsymbols];
 664         }
 665         map[k].rgb.r = map2[j].rgb.r;
 666         map[k].rgb.g = map2[j].rgb.g;
 667         map[k].rgb.b = map2[j].rgb.b;
 668         map[k].desc  = map2[j].desc;
 669     }
 670
 671     *newmap = map;
 672     return nmap;
 673 }
 674
 675 void xpm_mat(const char *outf, int nmat, t_matrix *mat, t_matrix *mat2,
 676              gmx_bool bDiag, gmx_bool bFirstDiag)
 677 {
 678     FILE      *out;
 679     int        i, x, y, col;
 680     int        nmap;
 681     t_mapping *map = NULL;
 682
 683     out = gmx_ffopen(outf, "w");
 684
 685     for (i = 0; i < nmat; i++)
 686     {
 687         if (!mat2 || !diff_maps(mat[i].nmap, mat[i].map, mat2[i].nmap, mat2[i].map))
 688         {
 689             write_xpm_m(out, mat[0]);
 690         }
 691         else
 692         {
 693             nmap = add_maps(&map, mat[i].nmap, mat[i].map, mat2[i].nmap, mat2[i].map);
 694             for (x = 0; (x < mat[i].nx); x++)
 695             {
 696                 for (y = 0; (y < mat[i].nx); y++)
 697                 {
 698                     if ((x < y) || ((x == y) && bFirstDiag)) /* upper left  -> map1 */
 699                     {
 700                         col = mat[i].matrix[x][y];
 701                     }
 702                     else /* lower right -> map2 */
 703                     {
 704                         col = mat[i].nmap+mat[i].matrix[x][y];
 705                     }
 706                     if ((bDiag) || (x != y))
 707                     {
 708                         mat[i].matrix[x][y] = col;
 709                     }
 710                     else
 711                     {
 712                         mat[i].matrix[x][y] = 0;
 713                     }
 714                 }
 715             }
 716             sfree(mat[i].map);
 717             mat[i].nmap = nmap;
 718             mat[i].map  = map;
 719             if (std::strcmp(mat[i].title, mat2[i].title) != 0)
 720             {
 721                 sprintf(mat[i].title+strlen(mat[i].title), " / %s", mat2[i].title);
 722             }
 723             if (std::strcmp(mat[i].legend, mat2[i].legend) != 0)
 724             {
 725                 sprintf(mat[i].legend+strlen(mat[i].legend), " / %s", mat2[i].legend);
 726             }
 727             write_xpm_m(out, mat[i]);
 728         }
 729     }
 730     gmx_ffclose(out);
 731 }
 732
 733 static void tick_spacing(int n, real axis[], real offset, char axisnm,
 734                          real *major, real *minor)
 735 {
 736     real     space;
 737     gmx_bool bTryAgain;
 738     int      i, j, t, f = 0, ten;
 739 #define NFACT 4
 740     real     major_fact[NFACT] = {5, 4, 2, 1};
 741     real     minor_fact[NFACT] = {5, 4, 4, 5};
 742
 743     /* start with interval between 10 matrix points: */
 744     space = std::max(10*axis[1]-axis[0], axis[std::min(10, n-1)]-axis[0]);
 745     /* get power of 10 */
 746     ten       = static_cast<int>(std::ceil(std::log(space)/std::log(10.0))-1);
 747     bTryAgain = TRUE;
 748     for (t = ten+2; t > ten-3 && bTryAgain; t--)
 749     {
 750         for (f = 0; f < NFACT && bTryAgain; f++)
 751         {
 752             space = std::pow(static_cast<real>(10.0), static_cast<real>(t)) * major_fact[f];
 753             /* count how many ticks we would get: */
 754             i = 0;
 755             for (j = 0; j < n; j++)
 756             {
 757                 if (bRmod(axis[j], offset, space) )
 758                 {
 759                     i++;
 760                 }
 761             }
 762             /* do we have a reasonable number of ticks ? */
 763             bTryAgain = (i > std::min(10, n-1)) || (i < 5);
 764         }
 765     }
 766     if (bTryAgain)
 767     {
 768         space = std::max(10*axis[1]-axis[0], axis[std::min(10, n-1)]-axis[0]);
 769         fprintf(stderr, "Auto tick spacing failed for %c-axis, guessing %g\n",
 770                 axisnm, space);
 771     }
 772     *major = space;
 773     *minor = space / minor_fact[(f > 0) ? f-1 : 0];
 774     fprintf(stderr, "Auto tick spacing for %c-axis: major %g, minor %g\n",
 775             axisnm, *major, *minor);
 776 }
 777
 778 void ps_mat(const char *outf, int nmat, t_matrix mat[], t_matrix mat2[],
 779             gmx_bool bFrame, gmx_bool bDiag, gmx_bool bFirstDiag,
 780             gmx_bool bTitle, gmx_bool bTitleOnce, gmx_bool bYonce, int elegend,
 781             real size, real boxx, real boxy, const char *m2p, const char *m2pout,
 782             int mapoffset)
 783 {
 784     const char   *libm2p;
 785     char          buf[256], *legend;
 786     t_psdata      out;
 787     t_psrec       psrec, *psr;
 788     int           W, H;
 789     int           i, x, y, col, leg = 0;
 790     real          x0, y0, xx;
 791     real          w, h, dw, dh;
 792     int           nmap1 = 0, nmap2 = 0, leg_nmap;
 793     t_mapping    *map1  = NULL, *map2 = NULL, *leg_map;
 794     gmx_bool      bMap1, bNextMap1, bDiscrete;
 795
 796     /* memory leak: */
 797     libm2p = m2p ? gmxlibfn(m2p) : m2p;
 798     get_params(libm2p, m2pout, &psrec);
 799
 800     psr = &psrec;
 801
 802     if (psr->X.major <= 0)
 803     {
 804         tick_spacing((mat[0].flags & MAT_SPATIAL_X) ? mat[0].nx + 1 : mat[0].nx,
 805                      mat[0].axis_x, psr->X.offset, 'X',
 806                      &(psr->X.major), &(psr->X.minor) );
 807     }
 808     if (psr->X.minor <= 0)
 809     {
 810         psr->X.minor = psr->X.major / 2;
 811     }
 812     if (psr->Y.major <= 0)
 813     {
 814         tick_spacing((mat[0].flags & MAT_SPATIAL_Y) ? mat[0].ny + 1 : mat[0].ny,
 815                      mat[0].axis_y, psr->Y.offset, 'Y',
 816                      &(psr->Y.major), &(psr->Y.minor) );
 817     }
 818     if (psr->Y.minor <= 0)
 819     {
 820         psr->Y.minor = psr->Y.major / 2;
 821     }
 822
 823     if (boxx > 0)
 824     {
 825         psr->xboxsize = boxx;
 826         psr->yboxsize = boxx;
 827     }
 828     if (boxy > 0)
 829     {
 830         psr->yboxsize = boxy;
 831     }
 832
 833     if (psr->xboxsize == 0)
 834     {
 835         psr->xboxsize = size/mat[0].nx;
 836         printf("Set the x-size of the box to %.3f\n", psr->xboxsize);
 837     }
 838     if (psr->yboxsize == 0)
 839     {
 840         psr->yboxsize = size/mat[0].nx;
 841         printf("Set the y-size of the box to %.3f\n", psr->yboxsize);
 842     }
 843
 844     nmap1 = 0;
 845     for (i = 0; (i < nmat); i++)
 846     {
 847         if (mat[i].nmap > nmap1)
 848         {
 849             nmap1 = mat[i].nmap;
 850             map1  = mat[i].map;
 851             leg   = i+1;
 852         }
 853     }
 854     if (leg != 1)
 855     {
 856         printf("Selected legend of matrix # %d for display\n", leg);
 857     }
 858     if (mat2)
 859     {
 860         nmap2 = 0;
 861         for (i = 0; (i < nmat); i++)
 862         {
 863             if (mat2[i].nmap > nmap2)
 864             {
 865                 nmap2 = mat2[i].nmap;
 866                 map2  = mat2[i].map;
 867                 leg   = i+1;
 868             }
 869         }
 870         if (leg != 1)
 871         {
 872             printf("Selected legend of matrix # %d for second display\n", leg);
 873         }
 874     }
 875     if ( (mat[0].legend[0] == 0) && psr->legend)
 876     {
 877         std::strcpy(mat[0].legend, psr->leglabel);
 878     }
 879
 880     bTitle          = bTitle     && mat[nmat-1].title[0];
 881     bTitleOnce      = bTitleOnce && mat[nmat-1].title[0];
 882     psr->bTitle     = bTitle;
 883     psr->bTitleOnce = bTitleOnce;
 884     psr->bYonce     = bYonce;
 885
 886     /* Set up size of box for nice colors */
 887     box_dim(nmat, mat, mat2, psr, elegend, bFrame, &w, &h, &dw, &dh);
 888
 889     /* Set up bounding box */
 890     W = static_cast<int>(w+dw);
 891     H = static_cast<int>(h+dh);
 892
 893     /* Start box at */
 894     x0 = dw;
 895     y0 = dh;
 896     x  = static_cast<int>(W+psr->xoffs);
 897     y  = static_cast<int>(H+psr->yoffs);
 898     if (bFrame)
 899     {
 900         x += 5*DDD;
 901         y += 4*DDD;
 902     }
 903     out = ps_open(outf, 0, 0, x, y);
 904     ps_linewidth(out, static_cast<int>(psr->linewidth));
 905     ps_init_rgb_box(out, psr->xboxsize, psr->yboxsize);
 906     ps_init_rgb_nbox(out, psr->xboxsize, psr->yboxsize);
 907     ps_translate(out, psr->xoffs, psr->yoffs);
 908
 909     if (bFrame)
 910     {
 911         ps_comment(out, "Here starts the BOX drawing");
 912         draw_boxes(out, x0, y0, w, nmat, mat, psr);
 913     }
 914
 915     for (i = 0; (i < nmat); i++)
 916     {
 917         if (bTitle || (bTitleOnce && i == nmat-1) )
 918         {
 919             /* Print title, if any */
 920             ps_rgb(out, BLACK);
 921             ps_strfont(out, psr->titfont, psr->titfontsize);
 922             if (!mat2 || (std::strcmp(mat[i].title, mat2[i].title) == 0))
 923             {
 924                 std::strcpy(buf, mat[i].title);
 925             }
 926             else
 927             {
 928                 sprintf(buf, "%s / %s", mat[i].title, mat2[i].title);
 929             }
 930             ps_ctext(out, x0+w/2, y0+box_height(&(mat[i]), psr)+psr->titfontsize,
 931                      buf, eXCenter);
 932         }
 933         sprintf(buf, "Here starts the filling of box #%d", i);
 934         ps_comment(out, buf);
 935         for (x = 0; (x < mat[i].nx); x++)
 936         {
 937             int nexty;
 938             int nextcol;
 939
 940             xx = x0+x*psr->xboxsize;
 941             ps_moveto(out, xx, y0);
 942             y     = 0;
 943             bMap1 = (!mat2 || (x < y || (x == y && bFirstDiag)));
 944             if ((bDiag) || (x != y))
 945             {
 946                 col = mat[i].matrix[x][y];
 947             }
 948             else
 949             {
 950                 col = -1;
 951             }
 952             for (nexty = 1; (nexty <= mat[i].ny); nexty++)
 953             {
 954                 bNextMap1 = (!mat2 || (x < nexty || (x == nexty && bFirstDiag)));
 955                 /* TRUE:  upper left  -> map1 */
 956                 /* FALSE: lower right -> map2 */
 957                 if ((nexty == mat[i].ny) || (!bDiag && (x == nexty)))
 958                 {
 959                     nextcol = -1;
 960                 }
 961                 else
 962                 {
 963                     nextcol = mat[i].matrix[x][nexty];
 964                 }
 965                 if ( (nexty == mat[i].ny) || (col != nextcol) || (bMap1 != bNextMap1) )
 966                 {
 967                     if (col >= 0)
 968                     {
 969                         if (bMap1)
 970                         {
 971                             ps_rgb_nbox(out, &(mat[i].map[col].rgb), nexty-y);
 972                         }
 973                         else
 974                         {
 975                             ps_rgb_nbox(out, &(mat2[i].map[col].rgb), nexty-y);
 976                         }
 977                     }
 978                     else
 979                     {
 980                         ps_moverel(out, 0, psr->yboxsize);
 981                     }
 982                     y     = nexty;
 983                     bMap1 = bNextMap1;
 984                     col   = nextcol;
 985                 }
 986             }
 987         }
 988         y0 += box_height(&(mat[i]), psr)+box_dh(psr)+box_dh_top(IS_ONCE, psr);
 989     }
 990
 991     if (psr->X.lineatzero || psr->Y.lineatzero)
 992     {
 993         /* reset y0 for first box */
 994         y0 = dh;
 995         ps_comment(out, "Here starts the zero lines drawing");
 996         draw_zerolines(out, x0, y0, w, nmat, mat, psr);
 997     }
 998
 999     if (elegend != elNone)
1000     {
1001         ps_comment(out, "Now it's legend time!");
1002         ps_linewidth(out, static_cast<int>(psr->linewidth));
1003         if (mat2 == NULL || elegend != elSecond)
1004         {
1005             bDiscrete = mat[0].bDiscrete;
1006             legend    = mat[0].legend;
1007             leg_nmap  = nmap1;
1008             leg_map   = map1;
1009         }
1010         else
1011         {
1012             bDiscrete = mat2[0].bDiscrete;
1013             legend    = mat2[0].legend;
1014             leg_nmap  = nmap2;
1015             leg_map   = map2;
1016         }
1017         if (bDiscrete)
1018         {
1019             leg_discrete(out, psr->legfontsize, DDD, legend,
1020                          psr->legfontsize, psr->legfont, leg_nmap, leg_map);
1021         }
1022         else
1023         {
1024             if (elegend != elBoth)
1025             {
1026                 leg_continuous(out, x0+w/2, w/2, DDD, legend,
1027                                psr->legfontsize, psr->legfont, leg_nmap, leg_map,
1028                                mapoffset);
1029             }
1030             else
1031             {
1032                 leg_bicontinuous(out, x0+w/2, w, DDD, mat[0].legend, mat2[0].legend,
1033                                  psr->legfontsize, psr->legfont, nmap1, map1, nmap2, map2);
1034             }
1035         }
1036         ps_comment(out, "Were there, dude");
1037     }
1038
1039     ps_close(out);
1040 }
1041
1042 void make_axis_labels(int nmat, t_matrix *mat)
1043 {
1044     int i, j;
1045
1046     for (i = 0; (i < nmat); i++)
1047     {
1048         /* Make labels for x axis */
1049         if (mat[i].axis_x == NULL)
1050         {
1051             snew(mat[i].axis_x, mat[i].nx);
1052             for (j = 0; (j < mat[i].nx); j++)
1053             {
1054                 mat[i].axis_x[j] = j;
1055             }
1056         }
1057         /* Make labels for y axis */
1058         if (mat[i].axis_y == NULL)
1059         {
1060             snew(mat[i].axis_y, mat[i].ny);
1061             for (j = 0; (j < mat[i].ny); j++)
1062             {
1063                 mat[i].axis_y[j] = j;
1064             }
1065         }
1066     }
1067 }
1068
1069 void prune_mat(int nmat, t_matrix *mat, t_matrix *mat2, int skip)
1070 {
1071     int i, x, y, xs, ys;
1072
1073     for (i = 0; i < nmat; i++)
1074     {
1075         fprintf(stderr, "converting %dx%d matrix to %dx%d\n",
1076                 mat[i].nx, mat[i].ny,
1077                 (mat[i].nx+skip-1)/skip, (mat[i].ny+skip-1)/skip);
1078         /* walk through matrix */
1079         xs = 0;
1080         for (x = 0; (x < mat[i].nx); x++)
1081         {
1082             if (x % skip == 0)
1083             {
1084                 mat[i].axis_x[xs] = mat[i].axis_x[x];
1085                 if (mat2)
1086                 {
1087                     mat2[i].axis_x[xs] = mat2[i].axis_x[x];
1088                 }
1089                 ys = 0;
1090                 for (y = 0; (y < mat[i].ny); y++)
1091                 {
1092                     if (x == 0)
1093                     {
1094                         mat[i].axis_y[ys] = mat[i].axis_y[y];
1095                         if (mat2)
1096                         {
1097                             mat2[i].axis_y[ys] = mat2[i].axis_y[y];
1098                         }
1099                     }
1100                     if (y % skip == 0)
1101                     {
1102                         mat[i].matrix[xs][ys] = mat[i].matrix[x][y];
1103                         if (mat2)
1104                         {
1105                             mat2[i].matrix[xs][ys] = mat2[i].matrix[x][y];
1106                         }
1107                         ys++;
1108                     }
1109                 }
1110                 xs++;
1111             }
1112         }
1113         /* adjust parameters */
1114         mat[i].nx = (mat[i].nx+skip-1)/skip;
1115         mat[i].ny = (mat[i].ny+skip-1)/skip;
1116         if (mat2)
1117         {
1118             mat2[i].nx = (mat2[i].nx+skip-1)/skip;
1119             mat2[i].ny = (mat2[i].ny+skip-1)/skip;
1120         }
1121     }
1122 }
1123
1124 void zero_lines(int nmat, t_matrix *mat, t_matrix *mat2)
1125 {
1126     int       i, x, y, m;
1127     t_matrix *mats;
1128
1129     for (i = 0; i < nmat; i++)
1130     {
1131         for (m = 0; m < (mat2 ? 2 : 1); m++)
1132         {
1133             if (m == 0)
1134             {
1135                 mats = mat;
1136             }
1137             else
1138             {
1139                 mats = mat2;
1140             }
1141             for (x = 0; x < mats[i].nx-1; x++)
1142             {
1143                 if (std::abs(mats[i].axis_x[x+1]) < 1e-5)
1144                 {
1145                     for (y = 0; y < mats[i].ny; y++)
1146                     {
1147                         mats[i].matrix[x][y] = 0;
1148                     }
1149                 }
1150             }
1151             for (y = 0; y < mats[i].ny-1; y++)
1152             {
1153                 if (std::abs(mats[i].axis_y[y+1]) < 1e-5)
1154                 {
1155                     for (x = 0; x < mats[i].nx; x++)
1156                     {
1157                         mats[i].matrix[x][y] = 0;
1158                     }
1159                 }
1160             }
1161         }
1162     }
1163 }
1164
1165 void write_combined_matrix(int ecombine, const char *fn,
1166                            int nmat, t_matrix *mat1, t_matrix *mat2,
1167                            real *cmin, real *cmax)
1168 {
1169     int        i, j, k, nlevels;
1170     FILE      *out;
1171     real     **rmat1, **rmat2;
1172     real       rhi, rlo;
1173
1174     out = gmx_ffopen(fn, "w");
1175     for (k = 0; k < nmat; k++)
1176     {
1177         if (mat2[k].nx != mat1[k].nx || mat2[k].ny != mat1[k].ny)
1178         {
1179             gmx_fatal(FARGS, "Size of frame %d in 1st (%dx%d) and 2nd matrix (%dx%d) do"
1180                       " not match.\n", k, mat1[k].nx, mat1[k].ny, mat2[k].nx, mat2[k].ny);
1181         }
1182         printf("Combining two %dx%d matrices\n", mat1[k].nx, mat1[k].ny);
1183         rmat1 = matrix2real(&mat1[k], NULL);
1184         rmat2 = matrix2real(&mat2[k], NULL);
1185         if (NULL == rmat1 || NULL == rmat2)
1186         {
1187             gmx_fatal(FARGS, "Could not extract real data from %s xpm matrices. Note that, e.g.,\n"
1188                       "g_rms and g_mdmat provide such data, but not do_dssp.\n",
1189                       (NULL == rmat1 && NULL == rmat2) ? "both" : "one of the" );
1190         }
1191         rlo = 1e38;
1192         rhi = -1e38;
1193         for (j = 0; j < mat1[k].ny; j++)
1194         {
1195             for (i = 0; i < mat1[k].nx; i++)
1196             {
1197                 switch (ecombine)
1198                 {
1199                     case ecAdd:  rmat1[i][j] += rmat2[i][j]; break;
1200                     case ecSub:  rmat1[i][j] -= rmat2[i][j]; break;
1201                     case ecMult: rmat1[i][j] *= rmat2[i][j]; break;
1202                     case ecDiv:  rmat1[i][j] /= rmat2[i][j]; break;
1203                     default:
1204                         gmx_fatal(FARGS, "No such combination rule %d for matrices", ecombine);
1205                 }
1206                 rlo = std::min(rlo, rmat1[i][j]);
1207                 rhi = std::max(rhi, rmat1[i][j]);
1208             }
1209         }
1210         if (cmin)
1211         {
1212             rlo = *cmin;
1213         }
1214         if (cmax)
1215         {
1216             rhi = *cmax;
1217         }
1218         nlevels = std::max(mat1[k].nmap, mat2[k].nmap);
1219         if (rhi == rlo)
1220         {
1221             fprintf(stderr,
1222                     "combination results in uniform matrix (%g), no output\n", rhi);
1223         }
1224         /*
1225            else if (rlo>=0 || rhi<=0)
1226            write_xpm(out, mat1[k].flags, mat1[k].title, mat1[k].legend,
1227             mat1[k].label_x, mat1[k].label_y,
1228             mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y,
1229             rmat1, rlo, rhi, rhi<=0?red:white, rhi<=0?white:blue,
1230             &nlevels);
1231            else
1232            write_xpm3(out, mat2[k].flags, mat1[k].title, mat1[k].legend,
1233              mat1[k].label_x, mat1[k].label_y,
1234              mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y,
1235              rmat1, rlo, 0, rhi, red, white, blue, &nlevels);
1236          */
1237         else
1238         {
1239             write_xpm(out, mat1[k].flags, mat1[k].title, mat1[k].legend,
1240                       mat1[k].label_x, mat1[k].label_y,
1241                       mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y,
1242                       rmat1, rlo, rhi, white, black, &nlevels);
1243         }
1244     }
1245     gmx_ffclose(out);
1246 }
1247
1248 void do_mat(int nmat, t_matrix *mat, t_matrix *mat2,
1249             gmx_bool bFrame, gmx_bool bZeroLine, gmx_bool bDiag, gmx_bool bFirstDiag, gmx_bool bTitle,
1250             gmx_bool bTitleOnce, gmx_bool bYonce, int elegend,
1251             real size, real boxx, real boxy,
1252             const char *epsfile, const char *xpmfile, const char *m2p,
1253             const char *m2pout, int skip, int mapoffset)
1254 {
1255     int      i, j, k;
1256
1257     if (mat2)
1258     {
1259         for (k = 0; (k < nmat); k++)
1260         {
1261             if ((mat2[k].nx != mat[k].nx) || (mat2[k].ny != mat[k].ny))
1262             {
1263                 gmx_fatal(FARGS, "WAKE UP!! Size of frame %d in 2nd matrix file (%dx%d) does not match size of 1st matrix (%dx%d) or the other way around.\n",
1264                           k, mat2[k].nx, mat2[k].ny, mat[k].nx, mat[k].ny);
1265             }
1266             for (j = 0; (j < mat[k].ny); j++)
1267             {
1268                 for (i = bFirstDiag ? j+1 : j; (i < mat[k].nx); i++)
1269                 {
1270                     mat[k].matrix[i][j] = mat2[k].matrix[i][j];
1271                 }
1272             }
1273         }
1274     }
1275     for (i = 0; (i < nmat); i++)
1276     {
1277         fprintf(stderr, "Matrix %d is %d x %d\n", i, mat[i].nx, mat[i].ny);
1278     }
1279
1280     make_axis_labels(nmat, mat);
1281
1282     if (skip > 1)
1283     {
1284         prune_mat(nmat, mat, mat2, skip);
1285     }
1286
1287     if (bZeroLine)
1288     {
1289         zero_lines(nmat, mat, mat);
1290     }
1291
1292     if (epsfile != NULL)
1293     {
1294         ps_mat(epsfile, nmat, mat, mat2, bFrame, bDiag, bFirstDiag,
1295                bTitle, bTitleOnce, bYonce, elegend,
1296                size, boxx, boxy, m2p, m2pout, mapoffset);
1297     }
1298     if (xpmfile != NULL)
1299     {
1300         xpm_mat(xpmfile, nmat, mat, mat2, bDiag, bFirstDiag);
1301     }
1302 }
1303
1304 void gradient_map(rvec grad, int nmap, t_mapping map[])
1305 {
1306     int  i;
1307     real c;
1308
1309     for (i = 0; i < nmap; i++)
1310     {
1311         c            = i/(nmap-1.0);
1312         map[i].rgb.r = 1-c*(1-grad[XX]);
1313         map[i].rgb.g = 1-c*(1-grad[YY]);
1314         map[i].rgb.b = 1-c*(1-grad[ZZ]);
1315     }
1316 }
1317
1318 void gradient_mat(rvec grad, int nmat, t_matrix mat[])
1319 {
1320     int m;
1321
1322     for (m = 0; m < nmat; m++)
1323     {
1324         gradient_map(grad, mat[m].nmap, mat[m].map);
1325     }
1326 }
1327
1328 void rainbow_map(gmx_bool bBlue, int nmap, t_mapping map[])
1329 {
1330     int  i;
1331     real c, r, g, b;
1332
1333     for (i = 0; i < nmap; i++)
1334     {
1335         c = (map[i].rgb.r + map[i].rgb.g + map[i].rgb.b)/3;
1336         if (c > 1)
1337         {
1338             c = 1;
1339         }
1340         if (bBlue)
1341         {
1342             c = 1 - c;
1343         }
1344         if (c <= 0.25) /* 0-0.25 */
1345         {
1346             r = 0;
1347             g = std::pow(4.0*c, 2.0/3.0);
1348             b = 1;
1349         }
1350         else if (c <= 0.5) /* 0.25-0.5 */
1351         {
1352             r = 0;
1353             g = 1;
1354             b = std::pow(2.0-4.0*c, 2.0/3.0);
1355         }
1356         else if (c <= 0.75) /* 0.5-0.75 */
1357         {
1358             r = std::pow(4.0*c-2.0, 2.0/3.0);
1359             g = 1;
1360             b = 0;
1361         }
1362         else /* 0.75-1 */
1363         {
1364             r = 1;
1365             g = std::pow(4.0-4.0*c, 2.0/3.0);
1366             b = 0;
1367         }
1368         map[i].rgb.r = r;
1369         map[i].rgb.g = g;
1370         map[i].rgb.b = b;
1371     }
1372 }
1373
1374 void rainbow_mat(gmx_bool bBlue, int nmat, t_matrix mat[])
1375 {
1376     int m;
1377
1378     for (m = 0; m < nmat; m++)
1379     {
1380         rainbow_map(bBlue, mat[m].nmap, mat[m].map);
1381     }
1382 }
1383
1384 int gmx_xpm2ps(int argc, char *argv[])
1385 {
1386     const char     *desc[] = {
1387         "[THISMODULE] makes a beautiful color plot of an XPixelMap file.",
1388         "Labels and axis can be displayed, when they are supplied",
1389         "in the correct matrix format.",
1390         "Matrix data may be generated by programs such as [gmx-do_dssp], [gmx-rms] or",
1391         "[gmx-mdmat].[PAR]",
1392         "Parameters are set in the [TT].m2p[tt] file optionally supplied with",
1393         "[TT]-di[tt]. Reasonable defaults are provided. Settings for the [IT]y[it]-axis",
1394         "default to those for the [IT]x[it]-axis. Font names have a defaulting hierarchy:",
1395         "titlefont -> legendfont; titlefont -> (xfont -> yfont -> ytickfont)",
1396         "-> xtickfont, e.g. setting titlefont sets all fonts, setting xfont",
1397         "sets yfont, ytickfont and xtickfont.[PAR]",
1398         "When no [TT].m2p[tt] file is supplied, many settings are taken from",
1399         "command line options. The most important option is [TT]-size[tt],",
1400         "which sets the size of the whole matrix in postscript units.",
1401         "This option can be overridden with the [TT]-bx[tt] and [TT]-by[tt]",
1402         "options (and the corresponding parameters in the [TT].m2p[tt] file),",
1403         "which set the size of a single matrix element.[PAR]",
1404         "With [TT]-f2[tt] a second matrix file can be supplied. Both matrix",
1405         "files will be read simultaneously and the upper left half of the",
1406         "first one ([TT]-f[tt]) is plotted together with the lower right",
1407         "half of the second one ([TT]-f2[tt]). The diagonal will contain",
1408         "values from the matrix file selected with [TT]-diag[tt].",
1409         "Plotting of the diagonal values can be suppressed altogether by",
1410         "setting [TT]-diag[tt] to [TT]none[tt].",
1411         "In this case, a new color map will be generated with",
1412         "a red gradient for negative numbers and a blue for positive.",
1413         "If the color coding and legend labels of both matrices are identical,",
1414         "only one legend will be displayed, else two separate legends are",
1415         "displayed.",
1416         "With [TT]-combine[tt], an alternative operation can be selected",
1417         "to combine the matrices. The output range is automatically set",
1418         "to the actual range of the combined matrix. This can be overridden",
1419         "with [TT]-cmin[tt] and [TT]-cmax[tt].[PAR]",
1420         "[TT]-title[tt] can be set to [TT]none[tt] to suppress the title, or to",
1421         "[TT]ylabel[tt] to show the title in the Y-label position (alongside",
1422         "the [IT]y[it]-axis).[PAR]",
1423         "With the [TT]-rainbow[tt] option, dull grayscale matrices can be turned",
1424         "into attractive color pictures.[PAR]",
1425         "Merged or rainbowed matrices can be written to an XPixelMap file with",
1426         "the [TT]-xpm[tt] option."
1427     };
1428
1429     output_env_t    oenv;
1430     const char     *fn, *epsfile = NULL, *xpmfile = NULL;
1431     int             i, nmat, nmat2, etitle, elegend, ediag, erainbow, ecombine;
1432     t_matrix       *mat = NULL, *mat2 = NULL;
1433     gmx_bool        bTitle, bTitleOnce, bDiag, bFirstDiag, bGrad;
1434     static gmx_bool bFrame = TRUE, bZeroLine = FALSE, bYonce = FALSE;
1435     static real     size   = 400, boxx = 0, boxy = 0, cmin = 0, cmax = 0;
1436     static rvec     grad   = {0, 0, 0};
1437     enum                    {
1438         etSel, etTop, etOnce, etYlabel, etNone, etNR
1439     };
1440     const char *title[]   = { NULL, "top", "once", "ylabel", "none", NULL };
1441     /* MUST correspond to enum elXxx as defined at top of file */
1442     const char *legend[]  = { NULL, "both", "first", "second", "none", NULL };
1443     enum                    {
1444         edSel, edFirst, edSecond, edNone, edNR
1445     };
1446     const char *diag[]    = { NULL, "first", "second", "none", NULL };
1447     enum                    {
1448         erSel, erNo, erBlue, erRed, erNR
1449     };
1450     const char *rainbow[] = { NULL, "no", "blue", "red", NULL };
1451     /* MUST correspond to enum ecXxx as defined at top of file */
1452     const char *combine[] = {
1453         NULL, "halves", "add", "sub", "mult", "div", NULL
1454     };
1455     static int  skip = 1, mapoffset = 0;
1456     t_pargs     pa[] = {
1457         { "-frame",   FALSE, etBOOL, {&bFrame},
1458           "Display frame, ticks, labels, title and legend" },
1459         { "-title",   FALSE, etENUM, {title},   "Show title at" },
1460         { "-yonce",   FALSE, etBOOL, {&bYonce}, "Show y-label only once" },
1461         { "-legend",  FALSE, etENUM, {legend},  "Show legend" },
1462         { "-diag",    FALSE, etENUM, {diag},    "Diagonal" },
1463         { "-size",    FALSE, etREAL, {&size},
1464           "Horizontal size of the matrix in ps units" },
1465         { "-bx",      FALSE, etREAL, {&boxx},
1466           "Element x-size, overrides [TT]-size[tt] (also y-size when [TT]-by[tt] is not set)" },
1467         { "-by",      FALSE, etREAL, {&boxy},   "Element y-size" },
1468         { "-rainbow", FALSE, etENUM, {rainbow},
1469           "Rainbow colors, convert white to" },
1470         { "-gradient", FALSE, etRVEC, {grad},
1471           "Re-scale colormap to a smooth gradient from white {1,1,1} to {r,g,b}" },
1472         { "-skip",    FALSE, etINT,  {&skip},
1473           "only write out every nr-th row and column" },
1474         { "-zeroline", FALSE, etBOOL, {&bZeroLine},
1475           "insert line in [REF].xpm[ref] matrix where axis label is zero"},
1476         { "-legoffset", FALSE, etINT, {&mapoffset},
1477           "Skip first N colors from [REF].xpm[ref] file for the legend" },
1478         { "-combine", FALSE, etENUM, {combine}, "Combine two matrices" },
1479         { "-cmin",    FALSE, etREAL, {&cmin}, "Minimum for combination output" },
1480         { "-cmax",    FALSE, etREAL, {&cmax}, "Maximum for combination output" }
1481     };
1482 #define NPA asize(pa)
1483     t_filenm    fnm[] = {
1484         { efXPM, "-f",  NULL,      ffREAD },
1485         { efXPM, "-f2", "root2",   ffOPTRD },
1486         { efM2P, "-di", NULL,      ffLIBOPTRD },
1487         { efM2P, "-do", "out",     ffOPTWR },
1488         { efEPS, "-o",  NULL,      ffOPTWR },
1489         { efXPM, "-xpm", NULL,      ffOPTWR }
1490     };
1491 #define NFILE asize(fnm)
1492
1493     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW,
1494                            NFILE, fnm, NPA, pa,
1495                            asize(desc), desc, 0, NULL, &oenv))
1496     {
1497         return 0;
1498     }
1499
1500     etitle   = nenum(title);
1501     elegend  = nenum(legend);
1502     ediag    = nenum(diag);
1503     erainbow = nenum(rainbow);
1504     ecombine = nenum(combine);
1505     bGrad    = opt2parg_bSet("-gradient", NPA, pa);
1506     for (i = 0; i < DIM; i++)
1507     {
1508         if (grad[i] < 0 || grad[i] > 1)
1509         {
1510             gmx_fatal(FARGS, "RGB value %g out of range (0.0-1.0)", grad[i]);
1511         }
1512     }
1513     if (!bFrame)
1514     {
1515         etitle  = etNone;
1516         elegend = elNone;
1517     }
1518
1519     epsfile = ftp2fn_null(efEPS, NFILE, fnm);
1520     xpmfile = opt2fn_null("-xpm", NFILE, fnm);
1521     if (epsfile == NULL && xpmfile == NULL)
1522     {
1523         if (ecombine != ecHalves)
1524         {
1525             xpmfile = opt2fn("-xpm", NFILE, fnm);
1526         }
1527         else
1528         {
1529             epsfile = ftp2fn(efEPS, NFILE, fnm);
1530         }
1531     }
1532     if (ecombine != ecHalves && epsfile)
1533     {
1534         fprintf(stderr,
1535                 "WARNING: can only write result of arithmetic combination "
1536                 "of two matrices to .xpm file\n"
1537                 "         file %s will not be written\n", epsfile);
1538         epsfile = NULL;
1539     }
1540
1541     bDiag      = ediag != edNone;
1542     bFirstDiag = ediag != edSecond;
1543
1544     fn   = opt2fn("-f", NFILE, fnm);
1545     nmat = read_xpm_matrix(fn, &mat);
1546     fprintf(stderr, "There %s %d matri%s in %s\n", (nmat > 1) ? "are" : "is", nmat, (nmat > 1) ? "ces" : "x", fn);
1547     fn = opt2fn_null("-f2", NFILE, fnm);
1548     if (fn)
1549     {
1550         nmat2 = read_xpm_matrix(fn, &mat2);
1551         fprintf(stderr, "There %s %d matri%s in %s\n", (nmat2 > 1) ? "are" : "is", nmat2, (nmat2 > 1) ? "ces" : "x", fn);
1552         if (nmat != nmat2)
1553         {
1554             fprintf(stderr, "Different number of matrices, using the smallest number.\n");
1555             nmat = nmat2 = std::min(nmat, nmat2);
1556         }
1557     }
1558     else
1559     {
1560         if (ecombine != ecHalves)
1561         {
1562             fprintf(stderr,
1563                     "WARNING: arithmetic matrix combination selected (-combine), "
1564                     "but no second matrix (-f2) supplied\n"
1565                     "         no matrix combination will be performed\n");
1566         }
1567         ecombine = 0;
1568         nmat2    = 0;
1569     }
1570     bTitle     = etitle == etTop;
1571     bTitleOnce = etitle == etOnce;
1572     if (etitle == etYlabel)
1573     {
1574         for (i = 0; (i < nmat); i++)
1575         {
1576             std::strcpy(mat[i].label_y, mat[i].title);
1577             if (mat2)
1578             {
1579                 std::strcpy(mat2[i].label_y, mat2[i].title);
1580             }
1581         }
1582     }
1583     if (bGrad)
1584     {
1585         gradient_mat(grad, nmat, mat);
1586         if (mat2)
1587         {
1588             gradient_mat(grad, nmat2, mat2);
1589         }
1590     }
1591     if (erainbow != erNo)
1592     {
1593         rainbow_mat(erainbow == erBlue, nmat, mat);
1594         if (mat2)
1595         {
1596             rainbow_mat(erainbow == erBlue, nmat2, mat2);
1597         }
1598     }
1599
1600     if ((mat2 == NULL) && (elegend != elNone))
1601     {
1602         elegend = elFirst;
1603     }
1604
1605     if (ecombine && ecombine != ecHalves)
1606     {
1607         write_combined_matrix(ecombine, xpmfile, nmat, mat, mat2,
1608                               opt2parg_bSet("-cmin", NPA, pa) ? &cmin : NULL,
1609                               opt2parg_bSet("-cmax", NPA, pa) ? &cmax : NULL);
1610     }
1611     else
1612     {
1613         do_mat(nmat, mat, mat2, bFrame, bZeroLine, bDiag, bFirstDiag,
1614                bTitle, bTitleOnce, bYonce,
1615                elegend, size, boxx, boxy, epsfile, xpmfile,
1616                opt2fn_null("-di", NFILE, fnm), opt2fn_null("-do", NFILE, fnm), skip,
1617                mapoffset);
1618     }
1619
1620     view_all(oenv, NFILE, fnm);
1621
1622     return 0;
1623 }