Remove segv in xpm2ps, fixes #611
[alexxy/gromacs.git] / 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     if (NULL == rmat1 || NULL == rmat2)
945     {
946         gmx_fatal(FARGS, "Could not extract real data from %s xpm matrices. Note that, e.g.,\n"
947                          "g_rms and g_mdmat provide such data, but not do_dssp.\n",
948                          (NULL==rmat1 && NULL==rmat2) ? "both" : "one of the" );
949     }
950     rlo=1e38;
951     rhi=-1e38;
952     for(j=0; j<mat1[k].ny; j++)
953       for(i=0; i<mat1[k].nx; i++) {
954         switch (ecombine) {
955         case ecAdd:  rmat1[i][j] += rmat2[i][j]; break;
956         case ecSub:  rmat1[i][j] -= rmat2[i][j]; break;
957         case ecMult: rmat1[i][j] *= rmat2[i][j]; break;
958         case ecDiv:  rmat1[i][j] /= rmat2[i][j]; break;
959         default:
960           gmx_fatal(FARGS,"No such combination rule %d for matrices",ecombine);
961         }
962         rlo = min(rlo, rmat1[i][j]);
963         rhi = max(rhi, rmat1[i][j]);
964       }
965     if (cmin)
966       rlo = *cmin;
967     if (cmax)
968       rhi = *cmax;
969     nlevels = max(mat1[k].nmap,mat2[k].nmap);
970     if (rhi==rlo)
971       fprintf(stderr,
972               "combination results in uniform matrix (%g), no output\n",rhi);
973     /*
974     else if (rlo>=0 || rhi<=0)
975       write_xpm(out, mat1[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, rhi, rhi<=0?red:white, rhi<=0?white:blue, 
979                 &nlevels);
980     else 
981       write_xpm3(out, mat2[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, 0, rhi, red, white, blue, &nlevels);
985     */
986     else
987       write_xpm(out, mat1[k].flags, mat1[k].title, mat1[k].legend, 
988                 mat1[k].label_x, mat1[k].label_y,
989                 mat1[k].nx, mat1[k].ny, mat1[k].axis_x, mat1[k].axis_y, 
990                 rmat1, rlo, rhi, white, black, &nlevels);       
991   }
992   ffclose(out);
993 }
994
995 void do_mat(int nmat,t_matrix *mat,t_matrix *mat2,
996             gmx_bool bFrame,gmx_bool bZeroLine,gmx_bool bDiag,gmx_bool bFirstDiag,gmx_bool bTitle,
997             gmx_bool bTitleOnce,gmx_bool bYonce,int elegend,
998             real size,real boxx,real boxy,
999             const char *epsfile,const char *xpmfile,const char *m2p,
1000             const char *m2pout,int skip, int mapoffset)
1001 {
1002   int      i,j,k;
1003
1004   if (mat2) {
1005     for (k=0; (k<nmat); k++) {
1006       if ((mat2[k].nx!=mat[k].nx) || (mat2[k].ny!=mat[k].ny)) 
1007         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",
1008                     k,mat2[k].nx,mat2[k].ny,mat[k].nx,mat[k].ny);
1009       for (j=0; (j<mat[k].ny); j++)
1010         for (i=bFirstDiag?j+1:j; (i<mat[k].nx); i++)
1011           mat[k].matrix[i][j]=mat2[k].matrix[i][j];
1012     }
1013   }
1014   for(i=0; (i<nmat); i++) 
1015     fprintf(stderr,"Matrix %d is %d x %d\n",i,mat[i].nx,mat[i].ny);
1016
1017   make_axis_labels(nmat, mat);
1018   
1019   if (skip > 1)
1020     prune_mat(nmat,mat,mat2,skip);
1021     
1022   if (bZeroLine)
1023     zero_lines(nmat,mat,mat);
1024   
1025   if (epsfile!=NULL)
1026     ps_mat(epsfile,nmat,mat,mat2,bFrame,bDiag,bFirstDiag,
1027            bTitle,bTitleOnce,bYonce,elegend,
1028            size,boxx,boxy,m2p,m2pout,mapoffset);
1029   if (xpmfile!=NULL)
1030     xpm_mat(xpmfile,nmat,mat,mat2,bDiag,bFirstDiag);
1031 }
1032
1033 void gradient_map(rvec grad, int nmap, t_mapping map[])
1034 {
1035   int i;
1036   real c;
1037   
1038   for(i=0; i<nmap; i++) {
1039     c = i/(nmap-1.0);
1040     map[i].rgb.r = 1-c*(1-grad[XX]);
1041     map[i].rgb.g = 1-c*(1-grad[YY]);
1042     map[i].rgb.b = 1-c*(1-grad[ZZ]);
1043   }
1044 }
1045   
1046 void gradient_mat(rvec grad, int nmat, t_matrix mat[])
1047 {
1048   int m;
1049   
1050   for(m=0; m<nmat; m++)
1051     gradient_map(grad, mat[m].nmap, mat[m].map);
1052 }
1053
1054 void rainbow_map(gmx_bool bBlue, int nmap, t_mapping map[])
1055 {
1056   int i;
1057   real c,r,g,b;
1058
1059   for(i=0; i<nmap; i++) {
1060     c = (map[i].rgb.r + map[i].rgb.g + map[i].rgb.b)/3;
1061     if (c > 1)
1062       c = 1;
1063     if (bBlue)
1064       c = 1 - c;
1065     if (c <= 0.25) { /* 0-0.25 */
1066       r = 0;
1067       g = pow(4*c,0.666);
1068       b = 1;
1069     } else if (c <= 0.5) { /* 0.25-0.5 */
1070       r = 0;
1071       g = 1;
1072       b = pow(2-4*c,0.666);
1073     } else if (c <= 0.75) { /* 0.5-0.75 */
1074       r = pow(4*c-2,0.666);
1075       g = 1;
1076       b = 0;
1077     } else { /* 0.75-1 */
1078       r = 1;
1079       g = pow(4-4*c,0.666);
1080       b = 0;
1081     }
1082     map[i].rgb.r = r;
1083     map[i].rgb.g = g;
1084     map[i].rgb.b = b;
1085   }
1086 }
1087
1088 void rainbow_mat(gmx_bool bBlue, int nmat, t_matrix mat[])
1089 {
1090   int m;
1091   
1092   for(m=0; m<nmat; m++)
1093     rainbow_map(bBlue, mat[m].nmap, mat[m].map);
1094 }
1095   
1096 int gmx_xpm2ps(int argc,char *argv[])
1097 {
1098   const char *desc[] = {
1099     "[TT]xpm2ps[tt] makes a beautiful color plot of an XPixelMap file.",
1100     "Labels and axis can be displayed, when they are supplied",
1101     "in the correct matrix format.",
1102     "Matrix data may be generated by programs such as [TT]do_dssp[tt], [TT]g_rms[tt] or",
1103     "[TT]g_mdmat[tt].[PAR]",
1104     "Parameters are set in the [TT].m2p[tt] file optionally supplied with",
1105     "[TT]-di[tt]. Reasonable defaults are provided. Settings for the [IT]y[it]-axis",
1106     "default to those for the [IT]x[it]-axis. Font names have a defaulting hierarchy:",
1107     "titlefont -> legendfont; titlefont -> (xfont -> yfont -> ytickfont)",
1108     "-> xtickfont, e.g. setting titlefont sets all fonts, setting xfont",
1109     "sets yfont, ytickfont and xtickfont.[PAR]",
1110     "When no [TT].m2p[tt] file is supplied, many settings are taken from",
1111     "command line options. The most important option is [TT]-size[tt],",
1112     "which sets the size of the whole matrix in postscript units.",
1113     "This option can be overridden with the [TT]-bx[tt] and [TT]-by[tt]",
1114     "options (and the corresponding parameters in the [TT].m2p[tt] file),",
1115     "which set the size of a single matrix element.[PAR]",
1116     "With [TT]-f2[tt] a second matrix file can be supplied. Both matrix",
1117     "files will be read simultaneously and the upper left half of the",
1118     "first one ([TT]-f[tt]) is plotted together with the lower right",
1119     "half of the second one ([TT]-f2[tt]). The diagonal will contain",
1120     "values from the matrix file selected with [TT]-diag[tt].",
1121     "Plotting of the diagonal values can be suppressed altogether by",
1122     "setting [TT]-diag[tt] to [TT]none[tt].",
1123     "In this case, a new color map will be generated with",
1124     "a red gradient for negative numbers and a blue for positive.",
1125     "If the color coding and legend labels of both matrices are identical,",
1126     "only one legend will be displayed, else two separate legends are",
1127     "displayed.",
1128     "With [TT]-combine[tt], an alternative operation can be selected",
1129     "to combine the matrices. The output range is automatically set",
1130     "to the actual range of the combined matrix. This can be overridden",
1131     "with [TT]-cmin[tt] and [TT]-cmax[tt].[PAR]",
1132     "[TT]-title[tt] can be set to [TT]none[tt] to suppress the title, or to",
1133     "[TT]ylabel[tt] to show the title in the Y-label position (alongside",
1134     "the [IT]y[it]-axis).[PAR]",
1135     "With the [TT]-rainbow[tt] option, dull grayscale matrices can be turned",
1136     "into attractive color pictures.[PAR]",
1137     "Merged or rainbowed matrices can be written to an XPixelMap file with",
1138     "the [TT]-xpm[tt] option."
1139   };
1140
1141   output_env_t oenv;
1142   const char *fn,*epsfile=NULL,*xpmfile=NULL;
1143   int       i,nmat,nmat2,etitle,elegend,ediag,erainbow,ecombine;
1144   t_matrix *mat=NULL,*mat2=NULL;
1145   gmx_bool      bTitle,bTitleOnce,bDiag,bFirstDiag,bGrad;
1146   static gmx_bool bFrame=TRUE,bZeroLine=FALSE,bYonce=FALSE,bAdd=FALSE;
1147   static real size=400,boxx=0,boxy=0,cmin=0,cmax=0;
1148   static rvec grad={0,0,0};
1149   enum                    { etSel, etTop, etOnce, etYlabel, etNone, etNR };
1150   const char *title[]   = { NULL, "top", "once", "ylabel", "none", NULL };
1151   /* MUST correspond to enum elXxx as defined at top of file */
1152   const char *legend[]  = { NULL, "both", "first", "second", "none", NULL };
1153   enum                    { edSel, edFirst, edSecond, edNone, edNR };
1154   const char *diag[]    = { NULL, "first", "second", "none", NULL };
1155   enum                    { erSel, erNo, erBlue, erRed, erNR };
1156   const char *rainbow[] = { NULL, "no", "blue", "red", NULL };
1157   /* MUST correspond to enum ecXxx as defined at top of file */
1158   const char *combine[] = {
1159     NULL, "halves", "add", "sub", "mult", "div", NULL };
1160   static int skip=1,mapoffset=0;
1161   t_pargs pa[] = {
1162     { "-frame",   FALSE, etBOOL, {&bFrame},
1163       "Display frame, ticks, labels, title and legend" },
1164     { "-title",   FALSE, etENUM, {title},   "Show title at" },
1165     { "-yonce",   FALSE, etBOOL, {&bYonce}, "Show y-label only once" },
1166     { "-legend",  FALSE, etENUM, {legend},  "Show legend" },
1167     { "-diag",    FALSE, etENUM, {diag},    "Diagonal" },
1168     { "-size",    FALSE, etREAL, {&size},
1169       "Horizontal size of the matrix in ps units" },
1170     { "-bx",      FALSE, etREAL, {&boxx},
1171       "Element x-size, overrides [TT]-size[tt] (also y-size when [TT]-by[tt] is not set)" },
1172     { "-by",      FALSE, etREAL, {&boxy},   "Element y-size" },
1173     { "-rainbow", FALSE, etENUM, {rainbow},
1174       "Rainbow colors, convert white to" },
1175     { "-gradient",FALSE, etRVEC, {grad},
1176       "Re-scale colormap to a smooth gradient from white {1,1,1} to {r,g,b}" },
1177     { "-skip",    FALSE, etINT,  {&skip},
1178       "only write out every nr-th row and column" },
1179     { "-zeroline",FALSE, etBOOL, {&bZeroLine},
1180       "insert line in [TT].xpm[tt] matrix where axis label is zero"},
1181     { "-legoffset", FALSE, etINT, {&mapoffset},
1182       "Skip first N colors from [TT].xpm[tt] file for the legend" },
1183     { "-combine", FALSE, etENUM, {combine}, "Combine two matrices" },
1184     { "-cmin",    FALSE, etREAL, {&cmin}, "Minimum for combination output" },
1185     { "-cmax",    FALSE, etREAL, {&cmax}, "Maximum for combination output" }
1186   };
1187 #define NPA asize(pa)
1188   t_filenm  fnm[] = {
1189     { efXPM, "-f",  NULL,      ffREAD },
1190     { efXPM, "-f2", "root2",   ffOPTRD },
1191     { efM2P, "-di", NULL,      ffLIBOPTRD },
1192     { efM2P, "-do", "out",     ffOPTWR },
1193     { efEPS, "-o",  NULL,      ffOPTWR },
1194     { efXPM, "-xpm",NULL,      ffOPTWR }
1195   };
1196 #define NFILE asize(fnm)
1197   
1198   CopyRight(stderr,argv[0]);
1199   parse_common_args(&argc,argv,PCA_CAN_VIEW,
1200                     NFILE,fnm,NPA,pa,
1201                     asize(desc),desc,0,NULL,&oenv);
1202
1203   etitle   = nenum(title);
1204   elegend  = nenum(legend);
1205   ediag    = nenum(diag);
1206   erainbow = nenum(rainbow);
1207   ecombine = nenum(combine);
1208   bGrad    = opt2parg_bSet("-gradient",NPA,pa);
1209   for(i=0; i<DIM; i++)
1210     if (grad[i] < 0 || grad[i] > 1)
1211       gmx_fatal(FARGS, "RGB value %g out of range (0.0-1.0)", grad[i]);
1212   if (!bFrame) {
1213     etitle = etNone;
1214     elegend = elNone;
1215   }
1216
1217   epsfile=ftp2fn_null(efEPS,NFILE,fnm);
1218   xpmfile=opt2fn_null("-xpm",NFILE,fnm);
1219   if ( epsfile==NULL && xpmfile==NULL ) {
1220     if (ecombine!=ecHalves)
1221       xpmfile=opt2fn("-xpm",NFILE,fnm);
1222     else
1223       epsfile=ftp2fn(efEPS,NFILE,fnm);
1224   }
1225   if (ecombine!=ecHalves && epsfile) {
1226     fprintf(stderr,
1227             "WARNING: can only write result of arithmetic combination "
1228             "of two matrices to .xpm file\n"
1229             "         file %s will not be written\n", epsfile);
1230     epsfile = NULL;
1231   }
1232   
1233   bDiag      = ediag!=edNone;
1234   bFirstDiag = ediag!=edSecond;
1235   
1236   fn=opt2fn("-f",NFILE,fnm);
1237   nmat=read_xpm_matrix(fn,&mat);
1238   fprintf(stderr,"There %s %d matri%s in %s\n",(nmat>1)?"are":"is",nmat,(nmat>1)?"ces":"x",fn);
1239   fn=opt2fn_null("-f2",NFILE,fnm);
1240   if (fn) {
1241     nmat2=read_xpm_matrix(fn,&mat2);
1242     fprintf(stderr,"There %s %d matri%s in %s\n",(nmat2>1)?"are":"is",nmat2,(nmat2>1)?"ces":"x",fn);
1243     if (nmat != nmat2) {
1244       fprintf(stderr,"Different number of matrices, using the smallest number.\n");
1245       nmat=nmat2=min(nmat,nmat2);
1246     }
1247   } else {
1248     if (ecombine!=ecHalves)
1249       fprintf(stderr,
1250               "WARNING: arithmetic matrix combination selected (-combine), "
1251               "but no second matrix (-f2) supplied\n"
1252               "         no matrix combination will be performed\n");
1253     ecombine=0;
1254     nmat2=0;
1255   }
1256   bTitle     = etitle==etTop;
1257   bTitleOnce = etitle==etOnce;
1258   if ( etitle==etYlabel )
1259     for (i=0; (i<nmat); i++) {
1260       strcpy(mat[i].label_y, mat[i].title);
1261       if (mat2)
1262         strcpy(mat2[i].label_y, mat2[i].title);
1263     }
1264   if (bGrad) {
1265     gradient_mat(grad,nmat,mat);
1266     if (mat2)
1267       gradient_mat(grad,nmat2,mat2);
1268   }
1269   if (erainbow!=erNo) {
1270     rainbow_mat(erainbow==erBlue,nmat,mat);
1271     if (mat2)
1272       rainbow_mat(erainbow==erBlue,nmat2,mat2);
1273   }
1274
1275   if ((mat2 == NULL) && (elegend!=elNone))
1276     elegend = elFirst;
1277   
1278   if (ecombine && ecombine!=ecHalves)
1279     write_combined_matrix(ecombine, xpmfile, nmat, mat, mat2,
1280                           opt2parg_bSet("-cmin",NPA,pa) ? &cmin : NULL,
1281                           opt2parg_bSet("-cmax",NPA,pa) ? &cmax : NULL);
1282   else
1283     do_mat(nmat,mat,mat2,bFrame,bZeroLine,bDiag,bFirstDiag,
1284            bTitle,bTitleOnce,bYonce,
1285            elegend,size,boxx,boxy,epsfile,xpmfile,
1286            opt2fn_null("-di",NFILE,fnm),opt2fn_null("-do",NFILE,fnm), skip,
1287            mapoffset);
1288   
1289   view_all(oenv,NFILE, fnm);
1290     
1291   thanx(stderr);
1292   
1293   return 0;
1294 }