src/tools/gmx_xpm2ps.c

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