Tidy: modernize-use-nullptr
[gromacs.git] / src / gromacs / fileio / matio.cpp
blob0fb307ca73b7dd7883c5abd92df01d82548847c4
1 /*
2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2017, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
37 #include "gmxpre.h"
39 #include "matio.h"
41 #include <cctype>
42 #include <cmath>
43 #include <cstdio>
44 #include <cstring>
46 #include <algorithm>
48 #include "gromacs/fileio/gmxfio.h"
49 #include "gromacs/math/utilities.h"
50 #include "gromacs/utility/binaryinformation.h"
51 #include "gromacs/utility/cstringutil.h"
52 #include "gromacs/utility/exceptions.h"
53 #include "gromacs/utility/fatalerror.h"
54 #include "gromacs/utility/futil.h"
55 #include "gromacs/utility/programcontext.h"
56 #include "gromacs/utility/smalloc.h"
58 static const char mapper[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789!@#$%^&*()-_=+{}|;:',<.>/?";
59 #define NMAP (long int)strlen(mapper)
61 #define MAX_XPM_LINELENGTH 4096
63 real **mk_matrix(int nx, int ny, gmx_bool b1D)
65 int i;
66 real **m;
68 snew(m, nx);
69 if (b1D)
71 snew(m[0], nx*ny);
74 for (i = 0; (i < nx); i++)
76 if (b1D)
78 m[i] = &(m[0][i*ny]);
80 else
82 snew(m[i], ny);
86 return m;
89 void done_matrix(int nx, real ***m)
91 int i;
93 for (i = 0; (i < nx); i++)
95 sfree((*m)[i]);
97 sfree(*m);
98 *m = nullptr;
101 gmx_bool matelmt_cmp(t_xpmelmt e1, t_xpmelmt e2)
103 return (e1.c1 == e2.c1) && (e1.c2 == e2.c2);
106 t_matelmt searchcmap(int n, t_mapping map[], t_xpmelmt c)
108 int i;
110 for (i = 0; (i < n); i++)
112 if (matelmt_cmp(map[i].code, c))
114 return i;
118 return -1;
121 int getcmap(FILE *in, const char *fn, t_mapping **map)
123 int i, n;
124 char line[STRLEN];
125 char code[STRLEN], desc[STRLEN];
126 double r, g, b;
127 t_mapping *m;
129 if (fgets2(line, STRLEN-1, in) == nullptr)
131 gmx_fatal(FARGS, "Not enough lines in colormap file %s"
132 "(just wanted to read number of entries)", fn);
134 sscanf(line, "%d", &n);
135 snew(m, n);
136 for (i = 0; (i < n); i++)
138 if (fgets2(line, STRLEN-1, in) == nullptr)
140 gmx_fatal(FARGS, "Not enough lines in colormap file %s"
141 "(should be %d, found only %d)", fn, n+1, i);
143 sscanf(line, "%s%s%lf%lf%lf", code, desc, &r, &g, &b);
144 m[i].code.c1 = code[0];
145 m[i].code.c2 = 0;
146 m[i].desc = gmx_strdup(desc);
147 m[i].rgb.r = r;
148 m[i].rgb.g = g;
149 m[i].rgb.b = b;
151 *map = m;
153 return n;
156 int readcmap(const char *fn, t_mapping **map)
158 FILE *in;
159 int n;
161 in = libopen(fn);
162 n = getcmap(in, fn, map);
163 gmx_ffclose(in);
165 return n;
168 void printcmap(FILE *out, int n, t_mapping map[])
170 int i;
172 fprintf(out, "%d\n", n);
173 for (i = 0; (i < n); i++)
175 fprintf(out, "%c%c %20s %10g %10g %10g\n",
176 map[i].code.c1 ? map[i].code.c1 : ' ',
177 map[i].code.c2 ? map[i].code.c2 : ' ',
178 map[i].desc, map[i].rgb.r, map[i].rgb.g, map[i].rgb.b);
182 void writecmap(const char *fn, int n, t_mapping map[])
184 FILE *out;
186 out = gmx_fio_fopen(fn, "w");
187 printcmap(out, n, map);
188 gmx_fio_fclose(out);
191 static char *fgetline(char **line, int llmax, int *llalloc, FILE *in)
193 char *fg;
195 if (llmax > *llalloc)
197 srenew(*line, llmax+1);
198 *llalloc = llmax;
200 fg = fgets(*line, llmax, in);
201 trim(*line);
203 return fg;
206 static void skipstr(char *line)
208 int i, c;
210 ltrim(line);
211 c = 0;
212 while ((line[c] != ' ') && (line[c] != '\0'))
214 c++;
216 i = c;
217 while (line[c] != '\0')
219 line[c-i] = line[c];
220 c++;
222 line[c-i] = '\0';
225 static char *line2string(char **line)
227 int i;
229 if (*line != nullptr)
231 while (((*line)[0] != '\"' ) && ( (*line)[0] != '\0' ))
233 (*line)++;
236 if ((*line)[0] != '\"')
238 return nullptr;
240 (*line)++;
242 i = 0;
243 while (( (*line)[i] != '\"' ) && ( (*line)[i] != '\0' ))
245 i++;
248 if ((*line)[i] != '\"')
250 *line = nullptr;
252 else
254 (*line)[i] = 0;
258 return *line;
261 static void parsestring(char *line, const char *label, char *string)
263 if (strstr(line, label))
265 if (strstr(line, label) < strchr(line, '\"'))
267 line2string(&line);
268 strcpy(string, line);
273 static void read_xpm_entry(FILE *in, t_matrix *mm)
275 t_mapping *map;
276 char *line_buf = nullptr, *line = nullptr, *str, buf[256] = {0};
277 int i, m, col_len, nch = 0, n_axis_x, n_axis_y, llmax;
278 int llalloc = 0;
279 unsigned int r, g, b;
280 double u;
281 gmx_bool bGetOnWithIt, bSetLine;
282 t_xpmelmt c;
284 mm->flags = 0;
285 mm->title[0] = 0;
286 mm->legend[0] = 0;
287 mm->label_x[0] = 0;
288 mm->label_y[0] = 0;
289 mm->matrix = nullptr;
290 mm->axis_x = nullptr;
291 mm->axis_y = nullptr;
292 mm->bDiscrete = FALSE;
294 llmax = STRLEN;
296 while ((nullptr != fgetline(&line_buf, llmax, &llalloc, in)) &&
297 (std::strncmp(line_buf, "static", 6) != 0))
299 line = line_buf;
300 parsestring(line, "title", (mm->title));
301 parsestring(line, "legend", (mm->legend));
302 parsestring(line, "x-label", (mm->label_x));
303 parsestring(line, "y-label", (mm->label_y));
304 parsestring(line, "type", buf);
307 if (!line || strncmp(line, "static", 6) != 0)
309 gmx_input("Invalid XPixMap");
312 if (buf[0] && (gmx_strcasecmp(buf, "Discrete") == 0))
314 mm->bDiscrete = TRUE;
317 if (debug)
319 fprintf(debug, "%s %s %s %s\n",
320 mm->title, mm->legend, mm->label_x, mm->label_y);
323 /* Read sizes */
324 bGetOnWithIt = FALSE;
325 while (!bGetOnWithIt && (nullptr != fgetline(&line_buf, llmax, &llalloc, in)))
327 line = line_buf;
328 while (( line[0] != '\"' ) && ( line[0] != '\0' ))
330 line++;
333 if (line[0] == '\"')
335 line2string(&line);
336 sscanf(line, "%d %d %d %d", &(mm->nx), &(mm->ny), &(mm->nmap), &nch);
337 if (nch > 2)
339 gmx_fatal(FARGS, "Sorry can only read xpm's with at most 2 caracters per pixel\n");
341 if (mm->nx <= 0 || mm->ny <= 0)
343 gmx_fatal(FARGS, "Dimensions of xpm-file have to be larger than 0\n");
345 llmax = std::max(STRLEN, mm->nx+10);
346 bGetOnWithIt = TRUE;
349 if (debug)
351 fprintf(debug, "mm->nx %d mm->ny %d mm->nmap %d nch %d\n",
352 mm->nx, mm->ny, mm->nmap, nch);
354 if (nch == 0)
356 gmx_fatal(FARGS, "Number of characters per pixel not found in xpm\n");
359 /* Read color map */
360 snew(map, mm->nmap);
361 m = 0;
362 while ((m < mm->nmap) && (nullptr != fgetline(&line_buf, llmax, &llalloc, in)))
364 line = std::strchr(line_buf, '\"');
365 if (line)
367 line++;
368 /* Read xpm color map entry */
369 map[m].code.c1 = line[0];
370 if (nch == 1)
372 map[m].code.c2 = 0;
374 else
376 map[m].code.c2 = line[1];
378 line += nch;
379 str = std::strchr(line, '#');
380 if (str)
382 str++;
383 col_len = 0;
384 while (std::isxdigit(str[col_len]))
386 col_len++;
388 if (col_len == 6)
390 sscanf(line, "%*s #%2x%2x%2x", &r, &g, &b);
391 map[m].rgb.r = r/255.0;
392 map[m].rgb.g = g/255.0;
393 map[m].rgb.b = b/255.0;
395 else if (col_len == 12)
397 sscanf(line, "%*s #%4x%4x%4x", &r, &g, &b);
398 map[m].rgb.r = r/65535.0;
399 map[m].rgb.g = g/65535.0;
400 map[m].rgb.b = b/65535.0;
402 else
404 gmx_file("Unsupported or invalid colormap in X PixMap");
407 else
409 str = std::strchr(line, 'c');
410 if (str)
412 str += 2;
414 else
416 gmx_file("Unsupported or invalid colormap in X PixMap");
418 fprintf(stderr, "Using white for color \"%s", str);
419 map[m].rgb.r = 1;
420 map[m].rgb.g = 1;
421 map[m].rgb.b = 1;
423 line = std::strchr(line, '\"');
424 line++;
425 line2string(&line);
426 map[m].desc = gmx_strdup(line);
427 m++;
430 if (m != mm->nmap)
432 gmx_fatal(FARGS, "Number of read colors map entries (%d) does not match the number in the header (%d)", m, mm->nmap);
434 mm->map = map;
436 /* Read axes, if there are any */
437 n_axis_x = 0;
438 n_axis_y = 0;
439 bSetLine = FALSE;
442 if (bSetLine)
444 line = line_buf;
446 bSetLine = TRUE;
447 if (strstr(line, "x-axis"))
449 line = std::strstr(line, "x-axis");
450 skipstr(line);
451 if (mm->axis_x == nullptr)
453 snew(mm->axis_x, mm->nx + 1);
455 while (sscanf(line, "%lf", &u) == 1)
457 if (n_axis_x > mm->nx)
459 gmx_fatal(FARGS, "Too many x-axis labels in xpm (max %d)", mm->nx);
461 else if (n_axis_x == mm->nx)
463 mm->flags |= MAT_SPATIAL_X;
465 mm->axis_x[n_axis_x] = u;
466 n_axis_x++;
467 skipstr(line);
470 else if (std::strstr(line, "y-axis"))
472 line = std::strstr(line, "y-axis");
473 skipstr(line);
474 if (mm->axis_y == nullptr)
476 snew(mm->axis_y, mm->ny + 1);
478 while (sscanf(line, "%lf", &u) == 1)
480 if (n_axis_y > mm->ny)
482 gmx_file("Too many y-axis labels in xpm");
484 else if (n_axis_y == mm->ny)
486 mm->flags |= MAT_SPATIAL_Y;
488 mm->axis_y[n_axis_y] = u;
489 n_axis_y++;
490 skipstr(line);
494 while ((line[0] != '\"') && (nullptr != fgetline(&line_buf, llmax, &llalloc, in)));
496 /* Read matrix */
497 snew(mm->matrix, mm->nx);
498 for (i = 0; i < mm->nx; i++)
500 snew(mm->matrix[i], mm->ny);
502 m = mm->ny-1;
503 bSetLine = FALSE;
506 if (bSetLine)
508 line = line_buf;
510 bSetLine = TRUE;
511 if (m%(1+mm->ny/100) == 0)
513 fprintf(stderr, "%3d%%\b\b\b\b", (100*(mm->ny-m))/mm->ny);
515 while ((line[0] != '\"') && (line[0] != '\0'))
517 line++;
519 if (line[0] != '\"')
521 gmx_fatal(FARGS, "Not enough caracters in row %d of the matrix\n", m+1);
523 else
525 line++;
526 for (i = 0; i < mm->nx; i++)
528 c.c1 = line[nch*i];
529 if (nch == 1)
531 c.c2 = 0;
533 else
535 c.c2 = line[nch*i+1];
537 mm->matrix[i][m] = searchcmap(mm->nmap, mm->map, c);
539 m--;
542 while ((m >= 0) && (nullptr != fgetline(&line_buf, llmax, &llalloc, in)));
543 if (m >= 0)
545 gmx_incons("Not enough rows in the matrix");
548 sfree(line_buf);
551 int read_xpm_matrix(const char *fnm, t_matrix **mat)
553 FILE *in;
554 char *line = nullptr;
555 int nmat;
556 int llalloc = 0;
558 in = gmx_fio_fopen(fnm, "r");
560 nmat = 0;
561 while (nullptr != fgetline(&line, STRLEN, &llalloc, in))
563 if (std::strstr(line, "/* XPM */"))
565 srenew(*mat, nmat+1);
566 read_xpm_entry(in, &(*mat)[nmat]);
567 nmat++;
570 gmx_fio_fclose(in);
572 if (nmat == 0)
574 gmx_file("Invalid XPixMap");
577 sfree(line);
579 return nmat;
582 real **matrix2real(t_matrix *in, real **out)
584 t_mapping *map;
585 double tmp;
586 real *rmap;
587 int i, j, nmap;
589 nmap = in->nmap;
590 map = in->map;
591 snew(rmap, nmap);
593 for (i = 0; i < nmap; i++)
595 if ((map[i].desc == nullptr) || (sscanf(map[i].desc, "%lf", &tmp) != 1))
597 fprintf(stderr, "Could not convert matrix to reals,\n"
598 "color map entry %d has a non-real description: \"%s\"\n",
599 i, map[i].desc);
600 sfree(rmap);
601 return nullptr;
603 rmap[i] = tmp;
606 if (out == nullptr)
608 snew(out, in->nx);
609 for (i = 0; i < in->nx; i++)
611 snew(out[i], in->ny);
614 for (i = 0; i < in->nx; i++)
616 for (j = 0; j < in->ny; j++)
618 out[i][j] = rmap[in->matrix[i][j]];
622 sfree(rmap);
624 fprintf(stderr, "Converted a %dx%d matrix with %d levels to reals\n",
625 in->nx, in->ny, nmap);
627 return out;
630 static void write_xpm_header(FILE *out,
631 const char *title, const char *legend,
632 const char *label_x, const char *label_y,
633 gmx_bool bDiscrete)
635 fprintf(out, "/* XPM */\n");
638 gmx::BinaryInformationSettings settings;
639 settings.generatedByHeader(true);
640 settings.linePrefix("/* ");
641 settings.lineSuffix(" */");
642 gmx::printBinaryInformation(out, gmx::getProgramContext(),
643 settings);
645 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
646 fprintf(out, "/* This file can be converted to EPS by the GROMACS program xpm2ps */\n");
647 fprintf(out, "/* title: \"%s\" */\n", title);
648 fprintf(out, "/* legend: \"%s\" */\n", legend);
649 fprintf(out, "/* x-label: \"%s\" */\n", label_x);
650 fprintf(out, "/* y-label: \"%s\" */\n", label_y);
651 if (bDiscrete)
653 fprintf(out, "/* type: \"Discrete\" */\n");
655 else
657 fprintf(out, "/* type: \"Continuous\" */\n");
661 static int calc_nmid(int nlevels, real lo, real mid, real hi)
663 /* Take care that we have at least 1 entry in the mid to hi range
665 return std::min(std::max(0, static_cast<int>(((mid-lo)/(hi-lo))*(nlevels-1))),
666 nlevels-1);
669 static void write_xpm_map3(FILE *out, int n_x, int n_y, int *nlevels,
670 real lo, real mid, real hi,
671 t_rgb rlo, t_rgb rmid, t_rgb rhi)
673 int i, nmid;
674 real r, g, b, clev_lo, clev_hi;
676 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;
682 else if (*nlevels < 2)
684 fprintf(stderr, "Warning, too few levels (%d) in matrix, using 2 instead\n",
685 *nlevels);
686 *nlevels = 2;
688 if (!((mid >= lo) && (mid < hi)))
690 gmx_fatal(FARGS, "Lo: %f, Mid: %f, Hi: %f\n", lo, mid, hi);
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);
697 nmid = calc_nmid(*nlevels, lo, mid, hi);
698 clev_lo = nmid;
699 clev_hi = (*nlevels - 1 - nmid);
700 for (i = 0; (i < nmid); i++)
702 r = rlo.r+(i*(rmid.r-rlo.r)/clev_lo);
703 g = rlo.g+(i*(rmid.g-rlo.g)/clev_lo);
704 b = rlo.b+(i*(rmid.b-rlo.b)/clev_lo);
705 fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
706 mapper[i % NMAP],
707 (*nlevels <= NMAP) ? ' ' : mapper[i/NMAP],
708 static_cast<unsigned int>(round(255*r)),
709 static_cast<unsigned int>(round(255*g)),
710 static_cast<unsigned int>(round(255*b)),
711 ((nmid - i)*lo + i*mid)/clev_lo);
713 for (i = 0; (i < (*nlevels-nmid)); i++)
715 r = rmid.r+(i*(rhi.r-rmid.r)/clev_hi);
716 g = rmid.g+(i*(rhi.g-rmid.g)/clev_hi);
717 b = rmid.b+(i*(rhi.b-rmid.b)/clev_hi);
718 fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
719 mapper[(i+nmid) % NMAP],
720 (*nlevels <= NMAP) ? ' ' : mapper[(i+nmid)/NMAP],
721 static_cast<unsigned int>(round(255*r)),
722 static_cast<unsigned int>(round(255*g)),
723 static_cast<unsigned int>(round(255*b)),
724 ((*nlevels - 1 - nmid - i)*mid + i*hi)/clev_hi);
728 static void pr_simple_cmap(FILE *out, real lo, real hi, int nlevel, t_rgb rlo,
729 t_rgb rhi, int i0)
731 int i;
732 real r, g, b, fac;
734 for (i = 0; (i < nlevel); i++)
736 fac = (i+1.0)/(nlevel);
737 r = rlo.r+fac*(rhi.r-rlo.r);
738 g = rlo.g+fac*(rhi.g-rlo.g);
739 b = rlo.b+fac*(rhi.b-rlo.b);
740 fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
741 mapper[(i+i0) % NMAP],
742 (nlevel <= NMAP) ? ' ' : mapper[(i+i0)/NMAP],
743 static_cast<unsigned int>(round(255*r)),
744 static_cast<unsigned int>(round(255*g)),
745 static_cast<unsigned int>(round(255*b)),
746 lo+fac*(hi-lo));
750 static void pr_discrete_cmap(FILE *out, int *nlevel, int i0)
752 t_rgb rgbd[16] = {
753 { 1.0, 1.0, 1.0 }, /* white */
754 { 1.0, 0.0, 0.0 }, /* red */
755 { 1.0, 1.0, 0.0 }, /* yellow */
756 { 0.0, 0.0, 1.0 }, /* blue */
757 { 0.0, 1.0, 0.0 }, /* green */
758 { 1.0, 0.0, 1.0 }, /* purple */
759 { 1.0, 0.4, 0.0 }, /* orange */
760 { 0.0, 1.0, 1.0 }, /* cyan */
761 { 1.0, 0.4, 0.4 }, /* pink */
762 { 1.0, 1.0, 0.0 }, /* yellow */
763 { 0.4, 0.4, 1.0 }, /* lightblue */
764 { 0.4, 1.0, 0.4 }, /* lightgreen */
765 { 1.0, 0.4, 1.0 }, /* lightpurple */
766 { 1.0, 0.7, 0.4 }, /* lightorange */
767 { 0.4, 1.0, 1.0 }, /* lightcyan */
768 { 0.0, 0.0, 0.0 } /* black */
771 int i, n;
773 *nlevel = std::min(16, *nlevel);
774 n = *nlevel;
775 for (i = 0; (i < n); i++)
777 fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%3d\" */,\n",
778 mapper[(i+i0) % NMAP],
779 (n <= NMAP) ? ' ' : mapper[(i+i0)/NMAP],
780 static_cast<unsigned int>(round(255*rgbd[i].r)),
781 static_cast<unsigned int>(round(255*rgbd[i].g)),
782 static_cast<unsigned int>(round(255*rgbd[i].b)),
789 static void write_xpm_map_split(FILE *out, int n_x, int n_y,
790 int *nlevel_top, real lo_top, real hi_top,
791 t_rgb rlo_top, t_rgb rhi_top,
792 gmx_bool bDiscreteColor,
793 int *nlevel_bot, real lo_bot, real hi_bot,
794 t_rgb rlo_bot, t_rgb rhi_bot)
796 int ntot;
798 ntot = *nlevel_top + *nlevel_bot;
799 if (ntot > NMAP)
801 gmx_fatal(FARGS, "Warning, too many levels (%d) in matrix", ntot);
804 fprintf(out, "static char *gromacs_xpm[] = {\n");
805 fprintf(out, "\"%d %d %d %d\",\n", n_x, n_y, ntot, 1);
807 if (bDiscreteColor)
809 pr_discrete_cmap(out, nlevel_bot, 0);
811 else
813 pr_simple_cmap(out, lo_bot, hi_bot, *nlevel_bot, rlo_bot, rhi_bot, 0);
816 pr_simple_cmap(out, lo_top, hi_top, *nlevel_top, rlo_top, rhi_top, *nlevel_bot);
820 static void write_xpm_map(FILE *out, int n_x, int n_y, int *nlevels,
821 real lo, real hi, t_rgb rlo, t_rgb rhi)
823 int i, nlo;
824 real invlevel, r, g, b;
826 if (*nlevels > NMAP*NMAP)
828 fprintf(stderr, "Warning, too many levels (%d) in matrix, using %d only\n",
829 *nlevels, (int)(NMAP*NMAP));
830 *nlevels = NMAP*NMAP;
832 else if (*nlevels < 2)
834 fprintf(stderr, "Warning, too few levels (%d) in matrix, using 2 instead\n", *nlevels);
835 *nlevels = 2;
838 fprintf(out, "static char *gromacs_xpm[] = {\n");
839 fprintf(out, "\"%d %d %d %d\",\n",
840 n_x, n_y, *nlevels, (*nlevels <= NMAP) ? 1 : 2);
842 invlevel = 1.0/(*nlevels-1);
843 for (i = 0; (i < *nlevels); i++)
845 nlo = *nlevels-1-i;
846 r = (nlo*rlo.r+i*rhi.r)*invlevel;
847 g = (nlo*rlo.g+i*rhi.g)*invlevel;
848 b = (nlo*rlo.b+i*rhi.b)*invlevel;
849 fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%.3g\" */,\n",
850 mapper[i % NMAP], (*nlevels <= NMAP) ? ' ' : mapper[i/NMAP],
851 static_cast<unsigned int>(round(255*r)),
852 static_cast<unsigned int>(round(255*g)),
853 static_cast<unsigned int>(round(255*b)),
854 (nlo*lo+i*hi)*invlevel);
858 static void write_xpm_axis(FILE *out, const char *axis, gmx_bool bSpatial,
859 int n, real *label)
861 int i;
863 if (label)
865 for (i = 0; i < (bSpatial ? n+1 : n); i++)
867 if (i % 80 == 0)
869 if (i)
871 fprintf(out, "*/\n");
873 fprintf(out, "/* %s-axis: ", axis);
875 fprintf(out, "%g ", label[i]);
877 fprintf(out, "*/\n");
881 static void write_xpm_data(FILE *out, int n_x, int n_y, real **mat,
882 real lo, real hi, int nlevels)
884 int i, j, c;
885 real invlevel;
887 invlevel = (nlevels-1)/(hi-lo);
888 for (j = n_y-1; (j >= 0); j--)
890 if (j%(1+n_y/100) == 0)
892 fprintf(stderr, "%3d%%\b\b\b\b", (100*(n_y-j))/n_y);
894 fprintf(out, "\"");
895 for (i = 0; (i < n_x); i++)
897 c = std::round((mat[i][j]-lo)*invlevel);
898 if (c < 0)
900 c = 0;
902 if (c >= nlevels)
904 c = nlevels-1;
906 if (nlevels <= NMAP)
908 fprintf(out, "%c", mapper[c]);
910 else
912 fprintf(out, "%c%c", mapper[c % NMAP], mapper[c / NMAP]);
915 if (j > 0)
917 fprintf(out, "\",\n");
919 else
921 fprintf(out, "\"\n");
926 static void write_xpm_data3(FILE *out, int n_x, int n_y, real **mat,
927 real lo, real mid, real hi, int nlevels)
929 int i, j, c = 0, nmid;
930 real invlev_lo, invlev_hi;
932 nmid = calc_nmid(nlevels, lo, mid, hi);
933 invlev_hi = (nlevels-1-nmid)/(hi-mid);
934 invlev_lo = (nmid)/(mid-lo);
936 for (j = n_y-1; (j >= 0); j--)
938 if (j%(1+n_y/100) == 0)
940 fprintf(stderr, "%3d%%\b\b\b\b", (100*(n_y-j))/n_y);
942 fprintf(out, "\"");
943 for (i = 0; (i < n_x); i++)
945 if (mat[i][j] >= mid)
947 c = nmid+std::round((mat[i][j]-mid)*invlev_hi);
949 else if (mat[i][j] >= lo)
951 c = std::round((mat[i][j]-lo)*invlev_lo);
953 else
955 c = 0;
958 if (c < 0)
960 c = 0;
962 if (c >= nlevels)
964 c = nlevels-1;
966 if (nlevels <= NMAP)
968 fprintf(out, "%c", mapper[c]);
970 else
972 fprintf(out, "%c%c", mapper[c % NMAP], mapper[c / NMAP]);
975 if (j > 0)
977 fprintf(out, "\",\n");
979 else
981 fprintf(out, "\"\n");
986 static void write_xpm_data_split(FILE *out, int n_x, int n_y, real **mat,
987 real lo_top, real hi_top, int nlevel_top,
988 real lo_bot, real hi_bot, int nlevel_bot)
990 int i, j, c;
991 real invlev_top, invlev_bot;
993 invlev_top = (nlevel_top-1)/(hi_top-lo_top);
994 invlev_bot = (nlevel_bot-1)/(hi_bot-lo_bot);
996 for (j = n_y-1; (j >= 0); j--)
998 if (j % (1+n_y/100) == 0)
1000 fprintf(stderr, "%3d%%\b\b\b\b", (100*(n_y-j))/n_y);
1002 fprintf(out, "\"");
1003 for (i = 0; (i < n_x); i++)
1005 if (i < j)
1007 c = nlevel_bot+round((mat[i][j]-lo_top)*invlev_top);
1008 if ((c < nlevel_bot) || (c >= nlevel_bot+nlevel_top))
1010 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, mat[i][j]);
1013 else if (i > j)
1015 c = round((mat[i][j]-lo_bot)*invlev_bot);
1016 if ((c < 0) || (c >= nlevel_bot+nlevel_bot))
1018 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, mat[i][j]);
1021 else
1023 c = nlevel_bot;
1026 fprintf(out, "%c", mapper[c]);
1028 if (j > 0)
1030 fprintf(out, "\",\n");
1032 else
1034 fprintf(out, "\"\n");
1039 void write_xpm_m(FILE *out, t_matrix m)
1041 /* Writes a t_matrix struct to .xpm file */
1043 int i, j;
1044 gmx_bool bOneChar;
1045 t_xpmelmt c;
1047 bOneChar = (m.map[0].code.c2 == 0);
1048 write_xpm_header(out, m.title, m.legend, m.label_x, m.label_y,
1049 m.bDiscrete);
1050 fprintf(out, "static char *gromacs_xpm[] = {\n");
1051 fprintf(out, "\"%d %d %d %d\",\n", m.nx, m.ny, m.nmap, bOneChar ? 1 : 2);
1052 for (i = 0; (i < m.nmap); i++)
1054 fprintf(out, "\"%c%c c #%02X%02X%02X \" /* \"%s\" */,\n",
1055 m.map[i].code.c1,
1056 bOneChar ? ' ' : m.map[i].code.c2,
1057 static_cast<unsigned int>(round(m.map[i].rgb.r*255)),
1058 static_cast<unsigned int>(round(m.map[i].rgb.g*255)),
1059 static_cast<unsigned int>(round(m.map[i].rgb.b*255)), m.map[i].desc);
1061 write_xpm_axis(out, "x", m.flags & MAT_SPATIAL_X, m.nx, m.axis_x);
1062 write_xpm_axis(out, "y", m.flags & MAT_SPATIAL_Y, m.ny, m.axis_y);
1063 for (j = m.ny-1; (j >= 0); j--)
1065 if (j%(1+m.ny/100) == 0)
1067 fprintf(stderr, "%3d%%\b\b\b\b", (100*(m.ny-j))/m.ny);
1069 fprintf(out, "\"");
1070 if (bOneChar)
1072 for (i = 0; (i < m.nx); i++)
1074 fprintf(out, "%c", m.map[m.matrix[i][j]].code.c1);
1077 else
1079 for (i = 0; (i < m.nx); i++)
1081 c = m.map[m.matrix[i][j]].code;
1082 fprintf(out, "%c%c", c.c1, c.c2);
1085 if (j > 0)
1087 fprintf(out, "\",\n");
1089 else
1091 fprintf(out, "\"\n");
1096 void write_xpm3(FILE *out, unsigned int flags,
1097 const char *title, const char *legend,
1098 const char *label_x, const char *label_y,
1099 int n_x, int n_y, real axis_x[], real axis_y[],
1100 real *mat[], real lo, real mid, real hi,
1101 t_rgb rlo, t_rgb rmid, t_rgb rhi, int *nlevels)
1103 /* See write_xpm.
1104 * Writes a colormap varying as rlo -> rmid -> rhi.
1107 if (hi <= lo)
1109 gmx_fatal(FARGS, "hi (%g) <= lo (%g)", hi, lo);
1112 write_xpm_header(out, title, legend, label_x, label_y, FALSE);
1113 write_xpm_map3(out, n_x, n_y, nlevels, lo, mid, hi, rlo, rmid, rhi);
1114 write_xpm_axis(out, "x", flags & MAT_SPATIAL_X, n_x, axis_x);
1115 write_xpm_axis(out, "y", flags & MAT_SPATIAL_Y, n_y, axis_y);
1116 write_xpm_data3(out, n_x, n_y, mat, lo, mid, hi, *nlevels);
1119 void write_xpm_split(FILE *out, unsigned int flags,
1120 const char *title, const char *legend,
1121 const char *label_x, const char *label_y,
1122 int n_x, int n_y, real axis_x[], real axis_y[],
1123 real *mat[],
1124 real lo_top, real hi_top, int *nlevel_top,
1125 t_rgb rlo_top, t_rgb rhi_top,
1126 real lo_bot, real hi_bot, int *nlevel_bot,
1127 gmx_bool bDiscreteColor,
1128 t_rgb rlo_bot, t_rgb rhi_bot)
1130 /* See write_xpm.
1131 * Writes a colormap varying as rlo -> rmid -> rhi.
1134 if (hi_top <= lo_top)
1136 gmx_fatal(FARGS, "hi_top (%g) <= lo_top (%g)", hi_top, lo_top);
1138 if (hi_bot <= lo_bot)
1140 gmx_fatal(FARGS, "hi_bot (%g) <= lo_bot (%g)", hi_bot, lo_bot);
1142 if (bDiscreteColor && (*nlevel_bot >= 16))
1144 gmx_impl("Can not plot more than 16 discrete colors");
1147 write_xpm_header(out, title, legend, label_x, label_y, FALSE);
1148 write_xpm_map_split(out, n_x, n_y, nlevel_top, lo_top, hi_top, rlo_top, rhi_top,
1149 bDiscreteColor, nlevel_bot, lo_bot, hi_bot, rlo_bot, rhi_bot);
1150 write_xpm_axis(out, "x", flags & MAT_SPATIAL_X, n_x, axis_x);
1151 write_xpm_axis(out, "y", flags & MAT_SPATIAL_Y, n_y, axis_y);
1152 write_xpm_data_split(out, n_x, n_y, mat, lo_top, hi_top, *nlevel_top,
1153 lo_bot, hi_bot, *nlevel_bot);
1156 void write_xpm(FILE *out, unsigned int flags,
1157 const char *title, const char *legend,
1158 const char *label_x, const char *label_y,
1159 int n_x, int n_y, real axis_x[], real axis_y[],
1160 real *mat[], real lo, real hi,
1161 t_rgb rlo, t_rgb rhi, int *nlevels)
1163 /* out xpm file
1164 * title matrix title
1165 * legend label for the continuous legend
1166 * label_x label for the x-axis
1167 * label_y label for the y-axis
1168 * n_x, n_y size of the matrix
1169 * axis_x[] the x-ticklabels
1170 * axis_y[] the y-ticklables
1171 * *matrix[] element x,y is matrix[x][y]
1172 * lo output lower than lo is set to lo
1173 * hi output higher than hi is set to hi
1174 * rlo rgb value for level lo
1175 * rhi rgb value for level hi
1176 * nlevels number of color levels for the output
1179 if (hi <= lo)
1181 gmx_fatal(FARGS, "hi (%f) <= lo (%f)", hi, lo);
1184 write_xpm_header(out, title, legend, label_x, label_y, FALSE);
1185 write_xpm_map(out, n_x, n_y, nlevels, lo, hi, rlo, rhi);
1186 write_xpm_axis(out, "x", flags & MAT_SPATIAL_X, n_x, axis_x);
1187 write_xpm_axis(out, "y", flags & MAT_SPATIAL_Y, n_y, axis_y);
1188 write_xpm_data(out, n_x, n_y, mat, lo, hi, *nlevels);