src/gromacs/gmxana/gmx_xpm2ps.c

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