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