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