src/gmxlib/matio.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  * GROningen Mixture of Alchemy and Childrens' Stories
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <stdio.h>
  40 #include <ctype.h>
  41 #include "sysstuff.h"
  42 #include "futil.h"
  43 #include "string2.h"
  44 #include "macros.h"
  45 #include "smalloc.h"
  46 #include "gmx_fatal.h"
  47 #include "matio.h"
  48 #include "statutil.h"
  49 #include "gmxfio.h"
  50 #include "maths.h"
  51
  52 #define round(a) (int)(a+0.5)
  53
  54 static const char mapper[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<.>/?";
  55 #define NMAP (long int)strlen(mapper)
  56
  57 #define MAX_XPM_LINELENGTH 4096
  58
  59 real **mk_matrix(int nx, int ny, gmx_bool b1D)
  60 {
  61   int  i;
  62   real **m;
  63
  64   snew(m,nx);
  65   if (b1D)
  66     snew(m[0], nx*ny);
  67
  68   for(i=0; (i<nx); i++)
  69     if (b1D)
  70       m[i] = &(m[0][i*ny]);
  71     else
  72       snew(m[i],ny);
  73
  74   return m;
  75 }
  76
  77 void done_matrix(int nx, real ***m)
  78 {
  79   int i;
  80
  81   for(i=0; (i<nx); i++)
  82     sfree((*m)[i]);
  83   sfree(*m);
  84   *m = NULL;
  85 }
  86
  87 void clear_matrix(int nx, int ny, real **m)
  88 {
  89   int x, y;
  90
  91   for(x=0; x<nx; x++)
  92     for(y=0; y<ny; y++)
  93       m[x][y]=0;
  94 }
  95
  96 gmx_bool matelmt_cmp(t_xpmelmt e1, t_xpmelmt e2)
  97 {
  98   return (e1.c1 == e2.c1) && (e1.c2 == e2.c2);
  99 }
 100
 101 t_matelmt searchcmap(int n,t_mapping map[],t_xpmelmt c)
 102 {
 103   int i;
 104
 105   for(i=0; (i<n); i++)
 106     if (matelmt_cmp(map[i].code, c))
 107       return i;
 108
 109   return -1;
 110 }
 111
 112 int getcmap(FILE *in,const char *fn,t_mapping **map)
 113 {
 114   int       i,n;
 115   char      line[STRLEN];
 116   char      code[STRLEN],desc[STRLEN];
 117   double    r,g,b;
 118   t_mapping *m;
 119
 120   if (fgets2(line,STRLEN-1,in) == NULL)
 121     gmx_fatal(FARGS,"Not enough lines in colormap file %s"
 122                 "(just wanted to read number of entries)",fn);
 123   sscanf(line,"%d",&n);
 124   snew(m,n);
 125   for(i=0; (i<n); i++) {
 126     if (fgets2(line,STRLEN-1,in) == NULL)
 127       gmx_fatal(FARGS,"Not enough lines in colormap file %s"
 128                   "(should be %d, found only %d)",fn,n+1,i);
 129     sscanf(line,"%s%s%lf%lf%lf",code,desc,&r,&g,&b);
 130     m[i].code.c1=code[0];
 131     m[i].code.c2=0;
 132     m[i].desc=strdup(desc);
 133     m[i].rgb.r=r;
 134     m[i].rgb.g=g;
 135     m[i].rgb.b=b;
 136   }
 137   *map=m;
 138
 139   return n;
 140 }
 141
 142 int readcmap(const char *fn,t_mapping **map)
 143 {
 144   FILE      *in;
 145   int       n;
 146
 147   in=libopen(fn);
 148   n=getcmap(in,fn,map);
 149   gmx_fio_fclose(in);
 150
 151   return n;
 152 }
 153
 154 void printcmap(FILE *out,int n,t_mapping map[])
 155 {
 156   int i;
 157
 158   fprintf(out,"%d\n",n);
 159   for(i=0; (i<n); i++)
 160     fprintf(out,"%c%c  %20s  %10g  %10g  %10g\n",
 161             map[i].code.c1?map[i].code.c1:' ',
 162             map[i].code.c2?map[i].code.c2:' ',
 163             map[i].desc,map[i].rgb.r,map[i].rgb.g,map[i].rgb.b);
 164 }
 165
 166 void writecmap(const char *fn,int n,t_mapping map[])
 167 {
 168   FILE *out;
 169
 170   out=gmx_fio_fopen(fn,"w");
 171   printcmap(out,n,map);
 172   gmx_fio_fclose(out);
 173 }
 174
 175 void do_wmap(FILE *out,int i0,int imax,
 176              int nlevels,t_rgb rlo,t_rgb rhi,real lo,real hi)
 177 {
 178   int  i,nlo;
 179   real r,g,b;
 180
 181   for(i=0; (i<imax); i++) {
 182     nlo=nlevels-i;
 183     r=(nlo*rlo.r+i*rhi.r)/nlevels;
 184     g=(nlo*rlo.g+i*rhi.g)/nlevels;
 185     b=(nlo*rlo.b+i*rhi.b)/nlevels;
 186     fprintf(out,"%c %10.3g %10g  %10g  %10g\n",
 187             mapper[i+i0],(nlo*lo+i*hi)/nlevels,r,g,b);
 188   }
 189 }
 190
 191 static char *fgetline(char **line,int llmax,int *llalloc,FILE *in)
 192 {
 193   char *fg;
 194
 195   if (llmax > *llalloc)
 196   {
 197       srenew(*line,llmax+1);
 198       *llalloc=llmax;
 199   }
 200   fg=fgets(*line,llmax,in);
 201   trim(*line);
 202
 203   return fg;
 204 }
 205
 206 void skipstr(char *line)
 207 {
 208   int i,c;
 209
 210   ltrim(line);
 211   c=0;
 212   while((line[c] != ' ') && (line[c] != '\0'))
 213     c++;
 214   i=c;
 215   while(line[c] != '\0')
 216     {
 217       line[c-i] = line[c];
 218       c++;
 219     }
 220   line[i-c] = '\0';
 221 }
 222
 223 char *line2string(char **line)
 224 {
 225   int i;
 226
 227   if (*line != NULL) {
 228     while (((*line)[0] != '\"' ) && ( (*line)[0] != '\0' ))
 229       (*line)++;
 230
 231     if ((*line)[0] != '\"')
 232       return NULL;
 233     (*line)++;
 234
 235       i=0;
 236     while (( (*line)[i] != '\"' ) && ( (*line)[i] != '\0' ))
 237       i++;
 238
 239     if ((*line)[i] != '\"')
 240       *line=NULL;
 241     else
 242       (*line)[i] = 0;
 243   }
 244
 245   return *line;
 246 }
 247
 248 void parsestring(char *line,const char *label, char *string)
 249 {
 250   if (strstr(line,label)) {
 251     if (strstr(line,label) < strchr(line,'\"')) {
 252       line2string(&line);
 253       strcpy(string,line);
 254     }
 255   }
 256 }
 257
 258 void read_xpm_entry(FILE *in,t_matrix *mm)
 259 {
 260   t_mapping *map;
 261   char *line=NULL,*str,buf[256];
 262   int i,m,col_len,nch,n_axis_x,n_axis_y,llmax;
 263   int llalloc=0;
 264   unsigned int r,g,b;
 265   double u;
 266   gmx_bool bGetOnWithIt;
 267   t_xpmelmt c;
 268
 269   mm->flags=0;
 270   mm->title[0]=0;
 271   mm->legend[0]=0;
 272   mm->label_x[0]=0;
 273   mm->label_y[0]=0;
 274   mm->matrix=NULL;
 275   mm->axis_x=NULL;
 276   mm->axis_y=NULL;
 277   mm->bDiscrete=FALSE;
 278
 279   llmax = STRLEN;
 280
 281   while ((NULL != fgetline(&line,llmax,&llalloc,in)) &&
 282          (strncmp(line,"static",6) != 0)) {
 283     parsestring(line,"title",(mm->title));
 284     parsestring(line,"legend",(mm->legend));
 285     parsestring(line,"x-label",(mm->label_x));
 286     parsestring(line,"y-label",(mm->label_y));
 287     parsestring(line,"type",buf);
 288   }
 289   if (buf[0] && (gmx_strcasecmp(buf,"Discrete")==0))
 290     mm->bDiscrete=TRUE;
 291
 292   if (debug)
 293     fprintf(debug,"%s %s %s %s\n",
 294             mm->title,mm->legend,mm->label_x,mm->label_y);
 295
 296   if  (strncmp(line,"static",6) != 0)
 297     gmx_input("Invalid XPixMap");
 298   /* Read sizes */
 299   bGetOnWithIt=FALSE;
 300   while (!bGetOnWithIt && (NULL != fgetline(&line,llmax,&llalloc,in))) {
 301     while (( line[0] != '\"' ) && ( line[0] != '\0' ))
 302       line++;
 303
 304     if  ( line[0] == '\"' ) {
 305       line2string(&line);
 306       sscanf(line,"%d %d %d %d",&(mm->nx),&(mm->ny),&(mm->nmap),&nch);
 307       if (nch > 2)
 308         gmx_fatal(FARGS,"Sorry can only read xpm's with at most 2 caracters per pixel\n");
 309       llmax = max(STRLEN,mm->nx+10);
 310       bGetOnWithIt=TRUE;
 311     }
 312   }
 313   if (debug)
 314     fprintf(debug,"mm->nx %d mm->ny %d mm->nmap %d nch %d\n",
 315            mm->nx,mm->ny,mm->nmap,nch);
 316
 317   /* Read color map */
 318   snew(map,mm->nmap);
 319   m=0;
 320   while ((m < mm->nmap) && (NULL != fgetline(&line,llmax,&llalloc,in))) {
 321     line=strchr(line,'\"');
 322     if  (line) {
 323       line++;
 324       /* Read xpm color map entry */
 325       map[m].code.c1 = line[0];
 326       if (nch==1)
 327         map[m].code.c2 = 0;
 328       else
 329         map[m].code.c2 = line[1];
 330       line += nch;
 331       str = strchr(line,'#');
 332       if (str) {
 333         str++;
 334         col_len = 0;
 335         while (isxdigit(str[col_len]))
 336           col_len++;
 337         if (col_len==6) {
 338           sscanf(line,"%*s #%2x%2x%2x",&r,&g,&b);
 339           map[m].rgb.r=r/255.0;
 340           map[m].rgb.g=g/255.0;
 341           map[m].rgb.b=b/255.0;
 342         } else if (col_len==12) {
 343           sscanf(line,"%*s #%4x%4x%4x",&r,&g,&b);
 344           map[m].rgb.r=r/65535.0;
 345           map[m].rgb.g=g/65535.0;
 346           map[m].rgb.b=b/65535.0;
 347         } else
 348           gmx_file("Unsupported or invalid colormap in X PixMap");
 349       } else {
 350         str = strchr(line,'c');
 351         if (str)
 352           str += 2;
 353         else
 354           gmx_file("Unsupported or invalid colormap in X PixMap");
 355         fprintf(stderr,"Using white for color \"%s",str);
 356         map[m].rgb.r = 1;
 357         map[m].rgb.g = 1;
 358         map[m].rgb.b = 1;
 359       }
 360       line=strchr(line,'\"');
 361       line++;
 362       line2string(&line);
 363       map[m].desc = strdup(line);
 364       m++;
 365     }
 366   }
 367   if  ( m != mm->nmap )
 368     gmx_fatal(FARGS,"Number of read colors map entries (%d) does not match the number in the header (%d)",m,mm->nmap);
 369   mm->map = map;
 370
 371   /* Read axes, if there are any */
 372   n_axis_x=0;
 373   n_axis_y=0;
 374   bGetOnWithIt=FALSE;
 375   do {
 376     if (strstr(line,"x-axis")) {
 377       line=strstr(line,"x-axis");
 378       skipstr(line);
 379       if (mm->axis_x==NULL)
 380         snew(mm->axis_x,mm->nx + 1);
 381       while (sscanf(line,"%lf",&u)==1) {
 382         if (n_axis_x > mm->nx) {
 383           gmx_fatal(FARGS,"Too many x-axis labels in xpm (max %d)",mm->nx);
 384         } else if (n_axis_x == mm->nx) {
 385           mm->flags |= MAT_SPATIAL_X;
 386         }
 387         mm->axis_x[n_axis_x] = u;
 388         n_axis_x++;
 389         skipstr(line);
 390       }
 391     }
 392     else if (strstr(line,"y-axis")) {
 393       line=strstr(line,"y-axis");
 394       skipstr(line);
 395       if (mm->axis_y==NULL)
 396         snew(mm->axis_y,mm->ny + 1);
 397       while (sscanf(line,"%lf",&u)==1) {
 398         if (n_axis_y > mm->ny) {
 399           gmx_file("Too many y-axis labels in xpm");
 400         } else if (n_axis_y == mm->ny) {
 401           mm->flags |= MAT_SPATIAL_Y;
 402         }
 403         mm->axis_y[n_axis_y] = u;
 404         n_axis_y++;
 405         skipstr(line);
 406       }
 407     }
 408   } while ((line[0] != '\"') && (NULL != fgetline(&line,llmax,&llalloc,in)));
 409
 410   /* Read matrix */
 411   snew(mm->matrix,mm->nx);
 412   for(i=0; i<mm->nx; i++)
 413     snew(mm->matrix[i],mm->ny);
 414   m=mm->ny-1;
 415   do {
 416     if(m%(1+mm->ny/100)==0)
 417       fprintf(stderr,"%3d%%\b\b\b\b",(100*(mm->ny-m))/mm->ny);
 418     while ((line[0] != '\"') && (line[0] != '\0'))
 419       line++;
 420     if (line[0] != '\"')
 421       gmx_fatal(FARGS,"Not enough caracters in row %d of the matrix\n",m+1);
 422     else {
 423       line++;
 424       for(i=0; i<mm->nx; i++) {
 425         c.c1=line[nch*i];
 426         if (nch==1)
 427           c.c2=0;
 428         else
 429           c.c2=line[nch*i+1];
 430         mm->matrix[i][m]=searchcmap(mm->nmap,mm->map,c);
 431         }
 432       m--;
 433     }
 434   } while ((m>=0) && (NULL != fgetline(&line,llmax,&llalloc,in)));
 435   if (m>=0)
 436     gmx_incons("Not enough rows in the matrix");
 437
 438   /* This code makes me cry. DvdS 2010-07-08 */
 439   /*sfree(line);*/
 440 }
 441
 442 int read_xpm_matrix(const char *fnm,t_matrix **matrix)
 443 {
 444   FILE *in;
 445   char *line=NULL;
 446   int nmat;
 447   int llalloc=0;
 448
 449   in=gmx_fio_fopen(fnm,"r");
 450
 451   nmat=0;
 452   while (NULL != fgetline(&line,STRLEN,&llalloc,in)) {
 453     if (strstr(line,"/* XPM */")) {
 454       srenew(*matrix,nmat+1);
 455       read_xpm_entry(in,&(*matrix)[nmat]);
 456       nmat++;
 457     }
 458   }
 459   gmx_fio_fclose(in);
 460
 461   if (nmat==0)
 462     gmx_file("Invalid XPixMap");
 463
 464   sfree(line);
 465
 466   return nmat;
 467 }
 468
 469 real **matrix2real(t_matrix *matrix,real **mat)
 470 {
 471   t_mapping *map;
 472   double tmp;
 473   real *rmap;
 474   int i,j,nmap;
 475
 476   nmap=matrix->nmap;
 477   map=matrix->map;
 478   snew(rmap,nmap);
 479
 480   for(i=0; i<nmap; i++) {
 481     if ((map[i].desc==NULL) || (sscanf(map[i].desc,"%lf",&tmp)!=1)) {
 482       fprintf(stderr,"Could not convert matrix to reals,\n"
 483               "color map entry %d has a non-real description: \"%s\"\n",
 484               i,map[i].desc);
 485       sfree(rmap);
 486       return NULL;
 487     }
 488     rmap[i]=tmp;
 489   }
 490
 491   if (mat==NULL) {
 492     snew(mat,matrix->nx);
 493     for(i=0; i<matrix->nx; i++)
 494       snew(mat[i],matrix->ny);
 495   }
 496   for(i=0; i<matrix->nx; i++)
 497     for(j=0; j<matrix->ny; j++)
 498       mat[i][j]=rmap[matrix->matrix[i][j]];
 499
 500   sfree(rmap);
 501
 502   fprintf(stderr,"Converted a %dx%d matrix with %d levels to reals\n",
 503           matrix->nx,matrix->ny,nmap);
 504
 505   return mat;
 506 }
 507
 508 void write_xpm_header(FILE *out,
 509                       const char *title,const char *legend,
 510                       const char *label_x,const char *label_y,
 511                       gmx_bool bDiscrete)
 512 {
 513   fprintf(out,  "/* XPM */\n");
 514   fprintf(out,  "/* Generated by %s */\n",Program());
 515   fprintf(out,  "/* This file can be converted to EPS by the GROMACS program xpm2ps */\n");
 516   fprintf(out,  "/* title:   \"%s\" */\n",title);
 517   fprintf(out,  "/* legend:  \"%s\" */\n",legend);
 518   fprintf(out,  "/* x-label: \"%s\" */\n",label_x);
 519   fprintf(out,  "/* y-label: \"%s\" */\n",label_y);
 520   if (bDiscrete)
 521     fprintf(out,"/* type:    \"Discrete\" */\n");
 522   else
 523     fprintf(out,"/* type:    \"Continuous\" */\n");
 524 }
 525
 526 static int calc_nmid(int nlevels,real lo,real mid,real hi)
 527 {
 528   /* Take care that we have at least 1 entry in the mid to hi range
 529    */
 530   return min(max(0,((mid-lo)/(hi-lo))*(nlevels-1)),nlevels-1);
 531 }
 532
 533 void write_xpm_map3(FILE *out,int n_x,int n_y,int *nlevels,
 534                     real lo,real mid,real hi,
 535                     t_rgb rlo,t_rgb rmid,t_rgb rhi)
 536 {
 537   int    i,nmid;
 538   real   r,g,b,clev_lo,clev_hi;
 539
 540   if (*nlevels > NMAP*NMAP) {
 541     fprintf(stderr,"Warning, too many levels (%d) in matrix, using %d only\n",
 542             *nlevels,(int)(NMAP*NMAP));
 543     *nlevels=NMAP*NMAP;
 544   }
 545   else if (*nlevels < 2) {
 546     fprintf(stderr,"Warning, too few levels (%d) in matrix, using 2 instead\n",
 547             *nlevels);
 548     *nlevels=2;
 549   }
 550   if (!((mid >= lo) && (mid < hi)))
 551     gmx_fatal(FARGS,"Lo: %f, Mid: %f, Hi: %f\n",lo,mid,hi);
 552
 553   fprintf(out,"static char *gromacs_xpm[] = {\n");
 554   fprintf(out,"\"%d %d   %d %d\",\n",
 555           n_x,n_y,*nlevels,(*nlevels <= NMAP) ? 1 : 2);
 556
 557   nmid    = calc_nmid(*nlevels,lo,mid,hi);
 558   clev_lo = nmid;
 559   clev_hi = (*nlevels - 1 - nmid);
 560   for(i=0; (i<nmid); i++) {
 561     r   = rlo.r+(i*(rmid.r-rlo.r)/clev_lo);
 562     g   = rlo.g+(i*(rmid.g-rlo.g)/clev_lo);
 563     b   = rlo.b+(i*(rmid.b-rlo.b)/clev_lo);
 564     fprintf(out,"\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 565             mapper[i % NMAP],
 566             (*nlevels <= NMAP) ? ' ' : mapper[i/NMAP],
 567             (unsigned int)round(255*r),
 568             (unsigned int)round(255*g),
 569             (unsigned int)round(255*b),
 570             ((nmid - i)*lo + i*mid)/clev_lo);
 571   }
 572   for(i=0; (i<(*nlevels-nmid)); i++) {
 573     r   = rmid.r+(i*(rhi.r-rmid.r)/clev_hi);
 574     g   = rmid.g+(i*(rhi.g-rmid.g)/clev_hi);
 575     b   = rmid.b+(i*(rhi.b-rmid.b)/clev_hi);
 576     fprintf(out,"\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 577             mapper[(i+nmid) % NMAP],
 578             (*nlevels <= NMAP) ? ' ' : mapper[(i+nmid)/NMAP],
 579             (unsigned int)round(255*r),
 580             (unsigned int)round(255*g),
 581             (unsigned int)round(255*b),
 582             ((*nlevels - 1 - nmid - i)*mid + i*hi)/clev_hi);
 583   }
 584 }
 585
 586 static void pr_simple_cmap(FILE *out,real lo,real hi,int nlevel,t_rgb rlo,
 587                            t_rgb rhi,int i0)
 588 {
 589   int    i;
 590   real   r,g,b,fac;
 591
 592   for(i=0; (i<nlevel); i++) {
 593     fac = (i+1.0)/(nlevel);
 594     r   = rlo.r+fac*(rhi.r-rlo.r);
 595     g   = rlo.g+fac*(rhi.g-rlo.g);
 596     b   = rlo.b+fac*(rhi.b-rlo.b);
 597     fprintf(out,"\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 598             mapper[(i+i0) % NMAP],
 599             (nlevel <= NMAP) ? ' ' : mapper[(i+i0)/NMAP],
 600             (unsigned int)round(255*r),
 601             (unsigned int)round(255*g),
 602             (unsigned int)round(255*b),
 603             lo+fac*(hi-lo));
 604   }
 605 }
 606
 607 static void pr_discrete_cmap(FILE *out,int *nlevel,int i0)
 608 {
 609   t_rgb  rgbd[16] = {
 610     { 1.0, 1.0, 1.0 }, /* white */
 611     { 1.0, 0.0, 0.0 }, /* red */
 612     { 1.0, 1.0, 0.0 }, /* yellow */
 613     { 0.0, 0.0, 1.0 }, /* blue */
 614     { 0.0, 1.0, 0.0 }, /* green */
 615     { 1.0, 0.0, 1.0 }, /* purple */
 616     { 1.0, 0.4, 0.0 }, /* orange */
 617     { 0.0, 1.0, 1.0 }, /* cyan */
 618     { 1.0, 0.4, 0.4 }, /* pink */
 619     { 1.0, 1.0, 0.0 }, /* yellow */
 620     { 0.4, 0.4, 1.0 }, /* lightblue */
 621     { 0.4, 1.0, 0.4 }, /* lightgreen */
 622     { 1.0, 0.4, 1.0 }, /* lightpurple */
 623     { 1.0, 0.7, 0.4 }, /* lightorange */
 624     { 0.4, 1.0, 1.0 }, /* lightcyan */
 625     { 0.0, 0.0, 0.0 }  /* black */
 626   };
 627
 628   int    i,n;
 629
 630   *nlevel = min(16,*nlevel);
 631   n = *nlevel;
 632   for(i=0; (i<n); i++) {
 633     fprintf(out,"\"%c%c c #%02X%02X%02X \" /* \"%3d\" */,\n",
 634             mapper[(i+i0) % NMAP],
 635             (n <= NMAP) ? ' ' : mapper[(i+i0)/NMAP],
 636             (unsigned int)round(255*rgbd[i].r),
 637             (unsigned int)round(255*rgbd[i].g),
 638             (unsigned int)round(255*rgbd[i].b),
 639             i);
 640   }
 641 }
 642
 643
 644
 645 void write_xpm_map_split(FILE *out,int n_x,int n_y,
 646                          int *nlevel_top,real lo_top,real hi_top,
 647                          t_rgb rlo_top,t_rgb rhi_top,
 648                          gmx_bool bDiscreteColor,
 649                          int *nlevel_bot,real lo_bot,real hi_bot,
 650                          t_rgb rlo_bot,t_rgb rhi_bot)
 651 {
 652   int    i,ntot;
 653   real   r,g,b,fac;
 654
 655   ntot = *nlevel_top + *nlevel_bot;
 656   if (ntot > NMAP)
 657     gmx_fatal(FARGS,"Warning, too many levels (%d) in matrix",ntot);
 658
 659   fprintf(out,"static char *gromacs_xpm[] = {\n");
 660   fprintf(out,"\"%d %d   %d %d\",\n",n_x,n_y,ntot,1);
 661
 662   if (bDiscreteColor)
 663     pr_discrete_cmap(out,nlevel_bot,0);
 664   else
 665     pr_simple_cmap(out,lo_bot,hi_bot,*nlevel_bot,rlo_bot,rhi_bot,0);
 666
 667   pr_simple_cmap(out,lo_top,hi_top,*nlevel_top,rlo_top,rhi_top,*nlevel_bot);
 668 }
 669
 670
 671 void write_xpm_map(FILE *out,int n_x, int n_y,int *nlevels,real lo,real hi,
 672                    t_rgb rlo,t_rgb rhi)
 673 {
 674   int    i,nlo;
 675   real   invlevel,r,g,b;
 676
 677   if (*nlevels > NMAP*NMAP) {
 678     fprintf(stderr,"Warning, too many levels (%d) in matrix, using %d only\n",
 679             *nlevels,(int)(NMAP*NMAP));
 680     *nlevels=NMAP*NMAP;
 681   }
 682   else if (*nlevels < 2) {
 683     fprintf(stderr,"Warning, too few levels (%d) in matrix, using 2 instead\n",*nlevels);
 684     *nlevels=2;
 685   }
 686
 687   fprintf(out,"static char *gromacs_xpm[] = {\n");
 688   fprintf(out,"\"%d %d   %d %d\",\n",
 689           n_x,n_y,*nlevels,(*nlevels <= NMAP) ? 1 : 2);
 690
 691   invlevel=1.0/(*nlevels-1);
 692   for(i=0; (i<*nlevels); i++) {
 693     nlo=*nlevels-1-i;
 694     r=(nlo*rlo.r+i*rhi.r)*invlevel;
 695     g=(nlo*rlo.g+i*rhi.g)*invlevel;
 696     b=(nlo*rlo.b+i*rhi.b)*invlevel;
 697     fprintf(out,"\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
 698             mapper[i % NMAP],(*nlevels <= NMAP) ? ' ' : mapper[i/NMAP],
 699             (unsigned int)round(255*r),
 700             (unsigned int)round(255*g),
 701             (unsigned int)round(255*b),
 702             (nlo*lo+i*hi)*invlevel);
 703   }
 704 }
 705
 706 void write_xpm_axis(FILE *out, const char *axis, gmx_bool bSpatial, int n,
 707                     real *label)
 708 {
 709   int i;
 710
 711   if (label) {
 712     for(i=0;i<(bSpatial ? n+1 : n);i++) {
 713       if (i % 80 == 0) {
 714         if (i)
 715           fprintf(out,"*/\n");
 716         fprintf(out,"/* %s-axis:  ",axis);
 717       }
 718       fprintf(out,"%g ",label[i]);
 719     }
 720     fprintf(out,"*/\n");
 721   }
 722 }
 723
 724 void write_xpm_data(FILE *out, int n_x, int n_y, real **matrix,
 725                     real lo, real hi, int nlevels)
 726 {
 727   int i,j,c;
 728   real invlevel;
 729
 730   invlevel=(nlevels-1)/(hi-lo);
 731   for(j=n_y-1; (j>=0); j--) {
 732     if(j%(1+n_y/100)==0)
 733       fprintf(stderr,"%3d%%\b\b\b\b",(100*(n_y-j))/n_y);
 734     fprintf(out,"\"");
 735     for(i=0; (i<n_x); i++) {
 736       c=gmx_nint((matrix[i][j]-lo)*invlevel);
 737       if (c<0) c=0;
 738       if (c>=nlevels) c=nlevels-1;
 739       if (nlevels <= NMAP)
 740         fprintf(out,"%c",mapper[c]);
 741       else
 742         fprintf(out,"%c%c",mapper[c % NMAP],mapper[c / NMAP]);
 743     }
 744     if (j > 0)
 745       fprintf(out,"\",\n");
 746     else
 747       fprintf(out,"\"\n");
 748   }
 749 }
 750
 751 void write_xpm_data3(FILE *out,int n_x,int n_y,real **matrix,
 752                     real lo,real mid,real hi,int nlevels)
 753 {
 754   int  i,j,c=0,nmid;
 755   real invlev_lo,invlev_hi;
 756
 757   nmid = calc_nmid(nlevels,lo,mid,hi);
 758   invlev_hi=(nlevels-1-nmid)/(hi-mid);
 759   invlev_lo=(nmid)/(mid-lo);
 760
 761   for(j=n_y-1; (j>=0); j--) {
 762     if(j%(1+n_y/100)==0)
 763       fprintf(stderr,"%3d%%\b\b\b\b",(100*(n_y-j))/n_y);
 764     fprintf(out,"\"");
 765     for(i=0; (i<n_x); i++) {
 766       if (matrix[i][j] >= mid)
 767         c=nmid+gmx_nint((matrix[i][j]-mid)*invlev_hi);
 768       else if (matrix[i][j] >= lo)
 769         c=gmx_nint((matrix[i][j]-lo)*invlev_lo);
 770       else
 771         c = 0;
 772
 773       if (c<0)
 774         c=0;
 775       if (c>=nlevels)
 776         c=nlevels-1;
 777       if (nlevels <= NMAP)
 778         fprintf(out,"%c",mapper[c]);
 779       else
 780         fprintf(out,"%c%c",mapper[c % NMAP],mapper[c / NMAP]);
 781     }
 782     if (j > 0)
 783       fprintf(out,"\",\n");
 784     else
 785       fprintf(out,"\"\n");
 786   }
 787 }
 788
 789 void write_xpm_data_split(FILE *out,int n_x,int n_y,real **matrix,
 790                           real lo_top,real hi_top,int nlevel_top,
 791                           real lo_bot,real hi_bot,int nlevel_bot)
 792 {
 793   int  i,j,c;
 794   real invlev_top,invlev_bot;
 795
 796   invlev_top=(nlevel_top-1)/(hi_top-lo_top);
 797   invlev_bot=(nlevel_bot-1)/(hi_bot-lo_bot);
 798
 799   for(j=n_y-1; (j>=0); j--) {
 800     if(j % (1+n_y/100)==0)
 801       fprintf(stderr,"%3d%%\b\b\b\b",(100*(n_y-j))/n_y);
 802     fprintf(out,"\"");
 803     for(i=0; (i<n_x); i++) {
 804       if (i < j) {
 805         c = nlevel_bot+round((matrix[i][j]-lo_top)*invlev_top);
 806         if ((c < nlevel_bot) || (c >= nlevel_bot+nlevel_top))
 807           gmx_fatal(FARGS,"Range checking i = %d, j = %d, c = %d, bot = %d, top = %d matrix[i,j] = %f",i,j,c,nlevel_bot,nlevel_top,matrix[i][j]);
 808       }
 809       else if (i > j) {
 810         c = round((matrix[i][j]-lo_bot)*invlev_bot);
 811         if ((c < 0) || (c >= nlevel_bot+nlevel_bot))
 812           gmx_fatal(FARGS,"Range checking i = %d, j = %d, c = %d, bot = %d, top = %d matrix[i,j] = %f",i,j,c,nlevel_bot,nlevel_top,matrix[i][j]);
 813       }
 814       else
 815         c = nlevel_bot;
 816
 817       fprintf(out,"%c",mapper[c]);
 818     }
 819     if (j > 0)
 820       fprintf(out,"\",\n");
 821     else
 822       fprintf(out,"\"\n");
 823   }
 824 }
 825
 826 void write_xpm_m(FILE *out, t_matrix m)
 827 {
 828   /* Writes a t_matrix struct to .xpm file */
 829
 830   int       i,j;
 831   gmx_bool      bOneChar;
 832   t_xpmelmt c;
 833
 834   bOneChar=(m.map[0].code.c2 == 0);
 835   write_xpm_header(out,m.title,m.legend,m.label_x,m.label_y,
 836                    m.bDiscrete);
 837   fprintf(out,"static char *gromacs_xpm[] = {\n");
 838   fprintf(out,"\"%d %d   %d %d\",\n",m.nx,m.ny,m.nmap,bOneChar ? 1 : 2);
 839   for(i=0; (i<m.nmap); i++)
 840     fprintf(out,"\"%c%c c #%02X%02X%02X \" /* \"%s\" */,\n",
 841             m.map[i].code.c1,
 842             bOneChar ? ' ' : m.map[i].code.c2,
 843             (unsigned int)round(m.map[i].rgb.r*255),
 844             (unsigned int)round(m.map[i].rgb.g*255),
 845             (unsigned int)round(m.map[i].rgb.b*255),m.map[i].desc);
 846   write_xpm_axis(out,"x",m.flags & MAT_SPATIAL_X,m.nx,m.axis_x);
 847   write_xpm_axis(out,"y",m.flags & MAT_SPATIAL_Y,m.ny,m.axis_y);
 848   for(j=m.ny-1; (j>=0); j--) {
 849     if(j%(1+m.ny/100)==0)
 850       fprintf(stderr,"%3d%%\b\b\b\b",(100*(m.ny-j))/m.ny);
 851     fprintf(out,"\"");
 852     if (bOneChar)
 853       for(i=0; (i<m.nx); i++)
 854         fprintf(out,"%c",m.map[m.matrix[i][j]].code.c1);
 855     else
 856       for(i=0; (i<m.nx); i++) {
 857         c=m.map[m.matrix[i][j]].code;
 858         fprintf(out,"%c%c",c.c1,c.c2);
 859       }
 860     if (j > 0)
 861       fprintf(out,"\",\n");
 862     else
 863       fprintf(out,"\"\n");
 864   }
 865 }
 866
 867 void write_xpm3(FILE *out,unsigned int flags,
 868                 const char *title,const char *legend,
 869                 const char *label_x,const char *label_y,
 870                 int n_x,int n_y,real axis_x[],real axis_y[],
 871                 real *matrix[],real lo,real mid,real hi,
 872                 t_rgb rlo,t_rgb rmid,t_rgb rhi,int *nlevels)
 873 {
 874   /* See write_xpm.
 875    * Writes a colormap varying as rlo -> rmid -> rhi.
 876    */
 877
 878   if (hi <= lo)
 879     gmx_fatal(FARGS,"hi (%g) <= lo (%g)",hi,lo);
 880
 881   write_xpm_header(out,title,legend,label_x,label_y,FALSE);
 882   write_xpm_map3(out,n_x,n_y,nlevels,lo,mid,hi,rlo,rmid,rhi);
 883   write_xpm_axis(out,"x",flags & MAT_SPATIAL_X,n_x,axis_x);
 884   write_xpm_axis(out,"y",flags & MAT_SPATIAL_Y,n_y,axis_y);
 885   write_xpm_data3(out,n_x,n_y,matrix,lo,mid,hi,*nlevels);
 886 }
 887
 888 void write_xpm_split(FILE *out,unsigned int flags,
 889                      const char *title,const char *legend,
 890                      const char *label_x,const char *label_y,
 891                      int n_x,int n_y,real axis_x[],real axis_y[],
 892                      real *matrix[],
 893                      real lo_top,real hi_top,int *nlevel_top,
 894                      t_rgb rlo_top,t_rgb rhi_top,
 895                      real lo_bot,real hi_bot,int *nlevel_bot,
 896                      gmx_bool bDiscreteColor,
 897                      t_rgb rlo_bot,t_rgb rhi_bot)
 898 {
 899   /* See write_xpm.
 900    * Writes a colormap varying as rlo -> rmid -> rhi.
 901    */
 902
 903   if (hi_top <= lo_top)
 904     gmx_fatal(FARGS,"hi_top (%g) <= lo_top (%g)",hi_top,lo_top);
 905   if (hi_bot <= lo_bot)
 906     gmx_fatal(FARGS,"hi_bot (%g) <= lo_bot (%g)",hi_bot,lo_bot);
 907   if (bDiscreteColor && (*nlevel_bot >= 16))
 908     gmx_impl("Can not plot more than 16 discrete colors");
 909
 910   write_xpm_header(out,title,legend,label_x,label_y,FALSE);
 911   write_xpm_map_split(out,n_x,n_y,nlevel_top,lo_top,hi_top,rlo_top,rhi_top,
 912                       bDiscreteColor,nlevel_bot,lo_bot,hi_bot,rlo_bot,rhi_bot);
 913   write_xpm_axis(out,"x",flags & MAT_SPATIAL_X,n_x,axis_x);
 914   write_xpm_axis(out,"y",flags & MAT_SPATIAL_Y,n_y,axis_y);
 915   write_xpm_data_split(out,n_x,n_y,matrix,lo_top,hi_top,*nlevel_top,
 916                        lo_bot,hi_bot,*nlevel_bot);
 917 }
 918
 919 void write_xpm(FILE *out,unsigned int flags,
 920                const char *title,const char *legend,
 921                const char *label_x,const char *label_y,
 922                int n_x,int n_y,real axis_x[],real axis_y[],
 923                real *matrix[],real lo,real hi,
 924                t_rgb rlo,t_rgb rhi,int *nlevels)
 925 {
 926   /* out        xpm file
 927    * title      matrix title
 928    * legend     label for the continuous legend
 929    * label_x    label for the x-axis
 930    * label_y    label for the y-axis
 931    * n_x, n_y   size of the matrix
 932    * axis_x[]   the x-ticklabels
 933    * axis_y[]   the y-ticklables
 934    * *matrix[]  element x,y is matrix[x][y]
 935    * lo         output lower than lo is set to lo
 936    * hi         output higher than hi is set to hi
 937    * rlo        rgb value for level lo
 938    * rhi        rgb value for level hi
 939    * nlevels    number of color levels for the output
 940    */
 941
 942   if (hi <= lo)
 943     gmx_fatal(FARGS,"hi (%f) <= lo (%f)",hi,lo);
 944
 945   write_xpm_header(out,title,legend,label_x,label_y,FALSE);
 946   write_xpm_map(out,n_x,n_y,nlevels,lo,hi,rlo,rhi);
 947   write_xpm_axis(out,"x",flags & MAT_SPATIAL_X,n_x,axis_x);
 948   write_xpm_axis(out,"y",flags & MAT_SPATIAL_Y,n_y,axis_y);
 949   write_xpm_data(out,n_x,n_y,matrix,lo,hi,*nlevels);
 950 }
 951