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