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