src/gromacs/gmxana/gmx_xpm2ps.c

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