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