MADNESS  version 0.9
funcplot.h
Go to the documentation of this file.
1 /*
2  This file is part of MADNESS.
3 
4  Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6  This program is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2 of the License, or
9  (at your option) any later version.
10 
11  This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with this program; if not, write to the Free Software
18  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20  For more information please contact:
21 
22  Robert J. Harrison
23  Oak Ridge National Laboratory
24  One Bethel Valley Road
25  P.O. Box 2008, MS-6367
26 
27  email: harrisonrj@ornl.gov
28  tel: 865-241-3937
29  fax: 865-572-0680
30 
31  $Id$
32 */
33 #ifndef MADNESS_MRA_FUNCPLOT_H__INCLUDED
34 #define MADNESS_MRA_FUNCPLOT_H__INCLUDED
35 
36 #include <madness/constants.h>
37 
47 namespace madness {
49 
66  template <typename T, std::size_t NDIM>
67  void plotdx(const Function<T,NDIM>& f,
68  const char* filename,
69  const Tensor<double>& cell = FunctionDefaults<NDIM>::get_cell(),
70  const std::vector<long>& npt = std::vector<long>(NDIM,201L),
71  bool binary=true);
72 
73 
76  //
83 
96  template<std::size_t NDIM>
97  void plotvtk_begin(World &world, const char *filename,
98  const Vector<double, NDIM> &plotlo, const Vector<double, NDIM> &plothi,
99  const Vector<long, NDIM> &npt, bool binary = false) {
100 
101  PROFILE_FUNC;
102  MADNESS_ASSERT(NDIM>=1 && NDIM<=3); // how do we plot data in more than 3-D?
103 
104  Tensor<double> cell(NDIM, 2);
105  std::size_t i;
106  for(i = 0; i < NDIM; ++i) {
107  cell(i, 0) = plotlo[i];
108  cell(i, 1) = plothi[i];
109  }
110 
111  FILE *f=0;
112  if(world.rank() == 0) {
113  f = fopen(filename, "w");
114  if(!f)
115  MADNESS_EXCEPTION("plotvtk: failed to open the plot file", 0);
116 
117  fprintf(f, "<VTKFile type=\"StructuredGrid\" version=\"0.1\"" \
118  " byte_order=\"LittleEndian\" compressor=\"" \
119  "vtkZLibDataCompressor\">\n");
120  fprintf(f, " <StructuredGrid WholeExtent=\"");
121  for(i = 0; i < NDIM; ++i)
122  fprintf(f, "0 %ld ", npt[i]-1);
123  for(; i < 3; ++i)
124  fprintf(f, "0 0 ");
125  fprintf(f, "\">\n");
126  fprintf(f, " <Piece Extent=\"");
127  for(i = 0; i < NDIM; ++i)
128  fprintf(f, "0 %ld ", npt[i]-1);
129  for(; i < 3; ++i)
130  fprintf(f, "0 0 ");
131  fprintf(f, "\">\n");
132  fprintf(f, " <Points>\n");
133  fprintf(f, " <DataArray NumberOfComponents=\"3\" " \
134  "type=\"Float32\" format=\"ascii\">\n");
135 
136  Vector<double, NDIM> space;
137  for(i = 0; i < NDIM; ++i) {
138  if(npt[i] == 1)
139  space[i] = 0.0;
140  else
141  space[i] = (cell(i, 1) - cell(i, 0)) / (npt[i] - 1);
142  }
143 
144  // go through the grid
145  for(LowDimIndexIterator it(npt); it; ++it) {
146  for(i = 0; i < NDIM; ++i)
147  fprintf(f, "%f ", plotlo[i] + it[i]*space[i]);
148  for(; i < 3; ++i)
149  fprintf(f, "0.0 ");
150  fprintf(f, "\n");
151  }
152 
153  fprintf(f, " </DataArray>\n");
154  fprintf(f, " </Points>\n");
155  fprintf(f, " <PointData>\n");
156  fclose(f);
157  }
158  world.gop.fence();
159  }
160 
163  //
172 
175  template<typename T, std::size_t NDIM>
176  void plotvtk_data(const T &function, const char *fieldname, World &world,
177  const char *filename, const Vector<double, NDIM> &plotlo,
178  const Vector<double, NDIM> &plothi, const Vector<long, NDIM> &npt,
179  bool binary = false) {
180 
181  MADNESS_EXCEPTION("plotvtk only supports madness::functions", 0);
182  }
183 
185 
188  template<typename T, std::size_t NDIM>
189  void plotvtk_data(const Function<T, NDIM> &function, const char *fieldname,
190  World &world, const char *filename, const Vector<double, NDIM> &plotlo,
191  const Vector<double, NDIM> &plothi, const Vector<long, NDIM> &npt,
192  bool binary = false, bool plot_refine = false) {
193 
194  PROFILE_FUNC;
195  MADNESS_ASSERT(NDIM>=1 && NDIM<=3); // no plotting high-D functions, yet...
196 
197  Tensor<double> cell(NDIM, 2);
198  std::size_t i;
199  for(i = 0; i < NDIM; ++i) {
200  cell(i, 0) = plotlo[i];
201  cell(i, 1) = plothi[i];
202  }
203  std::vector<long> numpt(NDIM);
204  for(i = 0; i < NDIM; ++i)
205  numpt[i] = npt[i];
206 
207  world.gop.barrier();
208 
209  function.verify();
210  FILE *f = 0;
211  if(world.rank() == 0) {
212  f = fopen(filename, "a");
213  if(!f)
214  MADNESS_EXCEPTION("plotvtk: failed to open the plot file", 0);
215 
216  fprintf(f, " <DataArray Name=\"%s\" format=\"ascii\" " \
217  "type=\"Float32\" NumberOfComponents=\"1\">\n", fieldname);
218  }
219 
220  world.gop.fence();
221  Tensor<T> tmpr = function.eval_cube(cell, numpt, plot_refine);
222  world.gop.fence();
223 
224  if(world.rank() == 0) {
225  for(LowDimIndexIterator it(numpt); it; ++it) {
226  fprintf(f, "%.6e\n", tmpr(*it));
227  }
228  fprintf(f, " </DataArray>\n");
229  fclose(f);
230  }
231  world.gop.fence();
232  }
233 
235 
241  template<typename T, std::size_t NDIM>
242  void plotvtk_data(const Function<std::complex<T>, NDIM> &function,
243  const char *fieldname, World &world, const char *filename,
244  const Vector<double, NDIM> &plotlo, const Vector<double, NDIM> &plothi,
245  const Vector<long, NDIM> &npt, bool binary = false,
246  bool plot_refine = false) {
247 
248  // this is the same as plotvtk_data for real functions, except the
249  // real and imaginary parts are printed on the same line (needed
250  // to change NumberOfComponents in the XML tag)
251 
252  PROFILE_FUNC;
253  MADNESS_ASSERT(NDIM>=1 && NDIM<=3); // no plotting high-D functions, yet...
254 
255  Tensor<double> cell(NDIM, 2);
256  std::size_t i;
257  for(i = 0; i < NDIM; ++i) {
258  cell(i, 0) = plotlo[i];
259  cell(i, 1) = plothi[i];
260  }
261  std::vector<long> numpt(NDIM);
262  for(i = 0; i < NDIM; ++i)
263  numpt[i] = npt[i];
264 
265  world.gop.barrier();
266 
267  function.verify();
268  FILE *f = 0;
269  if(world.rank() == 0) {
270  f = fopen(filename, "a");
271  if(!f)
272  MADNESS_EXCEPTION("plotvtk: failed to open the plot file", 0);
273 
274  fprintf(f, " <DataArray Name=\"%s\" format=\"ascii\" " \
275  "type=\"Float32\" NumberOfComponents=\"2\">\n", fieldname);
276  }
277 
278  world.gop.fence();
279  Tensor<std::complex<T> > tmpr = function.eval_cube(cell, numpt,
280  plot_refine);
281  world.gop.fence();
282 
283  if(world.rank() == 0) {
284  for(LowDimIndexIterator it(numpt); it; ++it) {
285  fprintf(f, "%.6e %.6e\n", real(tmpr(*it)), imag(tmpr(*it)));
286  }
287  fprintf(f, " </DataArray>\n");
288  fclose(f);
289  }
290  world.gop.fence();
291  }
292 
295  //
299  template<std::size_t NDIM>
300  void plotvtk_end(World &world, const char *filename, bool binary = false) {
301  PROFILE_FUNC;
302  MADNESS_ASSERT(NDIM>=1 && NDIM<=3);
303 
304  FILE *f = 0;
305  if(world.rank() == 0) {
306  f = fopen(filename, "a");
307  if(!f)
308  MADNESS_EXCEPTION("plotvtk: failed to open the plot file", 0);
309 
310  fprintf(f, " </PointData>\n");
311  fprintf(f, " <CellData>\n");
312  fprintf(f, " </CellData>\n");
313  fprintf(f, " </Piece>\n");
314  fprintf(f, " </StructuredGrid>\n");
315  fprintf(f, "</VTKFile>\n");
316  fclose(f);
317  }
318  world.gop.fence();
319  }
320 
321  namespace detail {
322  inline unsigned short htons_x(unsigned short a) {
323  return (a>>8) | (a<<8);
324  }
325  }
326 
328 
332  template <typename T>
333  static void plotpovray(const Function<T,3>& function,
334  const char* filename,
335  const Tensor<double>& cell = FunctionDefaults<3>::get_cell(),
336  const std::vector<long>& npt = std::vector<long>(3,201L))
337  {
338  using detail::htons_x;
339 
340  MADNESS_ASSERT(npt.size() == 3);
341  unsigned short dims[3] = {htons_x(npt[0]),htons_x(npt[1]),htons_x(npt[2])};
342 
343  World& world = const_cast< Function<T,3>& >(function).world();
344  FILE *f=0;
345  if (world.rank() == 0) {
346  f = fopen(filename, "w");
347  if (!f) MADNESS_EXCEPTION("plotdx: failed to open the plot file", 0);
348  fwrite((void*) dims, sizeof(short), 3, f);
349  }
350  Tensor<T> r = function.eval_cube(cell, npt);
351  if (world.rank() == 0) {
352  double rmax = r.max();
353  double rmin = r.min();
354  double rrange = rmax + rmin;
355  double rmean = rrange*0.5;
356  double fac = 65535.0/rrange;
357 
358  printf("plot_povray: %s: min=%.2e(0.0) mean=%.2e(0.5) max=%.2e(1.0) range=%.2e\n",
359  filename,rmin,rmean,rmax,rrange);
360 
361  std::vector<unsigned short> d(npt[0]);
362  for (unsigned int i2=0; i2<npt[2]; ++i2) {
363  for (unsigned int i1=0; i1<npt[1]; ++i1) {
364  for (unsigned int i0=0; i0<npt[0]; ++i0) {
365  d[i0] = (unsigned short)(htons_x((unsigned short)(fac*(r(i0,i1,i2) - rmin))));
366  //printf("%d\n",htons_x(d[i0]));
367  }
368  fwrite((void*) &d[0], sizeof(short), npt[0], f);
369  }
370  }
371 
372  fclose(f);
373  }
374  }
375 
376  static inline void plot_line_print_value(FILE* f, double_complex v) {
377  fprintf(f, " %.14e %.14e ", real(v), imag(v));
378  }
379 
380  static inline void plot_line_print_value(FILE* f, double v) {
381  fprintf(f, " %.14e", v);
382  }
383 
385 
387  template <typename T, std::size_t NDIM>
388  void plot_line(const char* filename, int npt, const Vector<double,NDIM>& lo, const Vector<double,NDIM>& hi,
389  const Function<T,NDIM>& f) {
390  typedef Vector<double,NDIM> coordT;
391  coordT h = (hi - lo)*(1.0/(npt-1));
392 
393  double sum = 0.0;
394  for (std::size_t i=0; i<NDIM; ++i) sum += h[i]*h[i];
395  sum = sqrt(sum);
396 
397  World& world = f.world();
398  f.reconstruct();
399  if (world.rank() == 0) {
400  FILE* file = fopen(filename,"w");
401  if(!file)
402  MADNESS_EXCEPTION("plot_line: failed to open the plot file", 0);
403  for (int i=0; i<npt; ++i) {
404  coordT r = lo + h*double(i);
405  fprintf(file, "%.14e ", i*sum);
406  plot_line_print_value(file, f.eval(r));
407  fprintf(file,"\n");
408  }
409  fclose(file);
410  }
411  world.gop.fence();
412  }
413 
415 
417  template <typename T, typename U, std::size_t NDIM>
418  void plot_line(const char* filename, int npt, const Vector<double,NDIM>& lo, const Vector<double,NDIM>& hi,
419  const Function<T,NDIM>& f, const Function<U,NDIM>& g) {
420  typedef Vector<double,NDIM> coordT;
421  coordT h = (hi - lo)*(1.0/(npt-1));
422 
423  double sum = 0.0;
424  for (std::size_t i=0; i<NDIM; ++i) sum += h[i]*h[i];
425  sum = sqrt(sum);
426 
427  World& world = f.world();
428  f.reconstruct();
429  g.reconstruct();
430  if (world.rank() == 0) {
431  FILE* file = fopen(filename,"w");
432  if(!file)
433  MADNESS_EXCEPTION("plot_line: failed to open the plot file", 0);
434  for (int i=0; i<npt; ++i) {
435  coordT r = lo + h*double(i);
436  fprintf(file, "%.14e ", i*sum);
437  plot_line_print_value(file, f.eval(r));
438  plot_line_print_value(file, g.eval(r));
439  fprintf(file,"\n");
440  }
441  fclose(file);
442  }
443  world.gop.fence();
444  }
445 
446 
448 
450  template <typename T, typename U, typename V, std::size_t NDIM>
451  void plot_line(const char* filename, int npt, const Vector<double,NDIM>& lo, const Vector<double,NDIM>& hi,
452  const Function<T,NDIM>& f, const Function<U,NDIM>& g, const Function<V,NDIM>& a) {
453  typedef Vector<double,NDIM> coordT;
454  coordT h = (hi - lo)*(1.0/(npt-1));
455 
456  double sum = 0.0;
457  for (std::size_t i=0; i<NDIM; ++i) sum += h[i]*h[i];
458  sum = sqrt(sum);
459 
460  World& world = f.world();
461  f.reconstruct();
462  g.reconstruct();
463  a.reconstruct();
464  if (world.rank() == 0) {
465  FILE* file = fopen(filename,"w");
466  if(!file)
467  MADNESS_EXCEPTION("plot_line: failed to open the plot file", 0);
468  for (int i=0; i<npt; ++i) {
469  coordT r = lo + h*double(i);
470  fprintf(file, "%.14e ", i*sum);
471  plot_line_print_value(file, f.eval(r));
472  plot_line_print_value(file, g.eval(r));
473  plot_line_print_value(file, a.eval(r));
474  fprintf(file,"\n");
475  }
476  fclose(file);
477  }
478  world.gop.fence();
479  }
480 
482 
484  template <typename T, typename U, typename V, typename W, std::size_t NDIM>
485  void plot_line(const char* filename, int npt, const Vector<double,NDIM>& lo, const Vector<double,NDIM>& hi,
486  const Function<T,NDIM>& f, const Function<U,NDIM>& g, const Function<V,NDIM>& a, const Function<W,NDIM>& b) {
487  typedef Vector<double,NDIM> coordT;
488  coordT h = (hi - lo)*(1.0/(npt-1));
489 
490  double sum = 0.0;
491  for (std::size_t i=0; i<NDIM; ++i) sum += h[i]*h[i];
492  sum = sqrt(sum);
493 
494  World& world = f.world();
495  f.reconstruct();
496  g.reconstruct();
497  a.reconstruct();
498  b.reconstruct();
499  if (world.rank() == 0) {
500  FILE* file = fopen(filename,"w");
501  for (int i=0; i<npt; ++i) {
502  coordT r = lo + h*double(i);
503  fprintf(file, "%.14e ", i*sum);
504  plot_line_print_value(file, f.eval(r));
505  plot_line_print_value(file, g.eval(r));
506  plot_line_print_value(file, a.eval(r));
507  plot_line_print_value(file, b.eval(r));
508  fprintf(file,"\n");
509  }
510  fclose(file);
511  }
512  world.gop.fence();
513  }
514 
515 
518 
530  template<size_t NDIM>
531  void plot_plane(World& world, const Function<double,NDIM>& function,
532  const std::string name) {
533 
534  if (world.size()>1) return;
535  // determine the ploting plane
536  std::string c1, c2;
537 
538  // zoom factor
539  double zoom=1.0;
540 
541  // output type: mathematica or gnuplot
542  std::string output_type="gnuplot";
543 
544  // number of points in each direction
545  int npoints=200;
546 
547  // the coordinates to be plotted
548  Vector<double,NDIM> coord(0.0);
549  Vector<double,NDIM> origin(0.0);
550 
551  std::ifstream f("input");
552  position_stream(f, "plot");
553  std::string s;
554  while (f >> s) {
555  if (s == "end") {
556  break;
557  } else if (s == "plane") {
558  f >> c1 >> c2;
559  } else if (s == "zoom") {
560  f >> zoom;
561  } else if (s == "output") {
562  f >> output_type;
563  } else if (s == "points") {
564  f >> npoints;
565  } else if (s == "origin") {
566  for (std::size_t i=0; i<NDIM; ++i) f >> origin[i];
567  }
568  }
569  double scale=1.0/zoom;
570  coord=origin;
571 
572  // convert human to mad form
573  int cc1=0, cc2=1;
574  if (c1=="x1") cc1=0;
575  if (c1=="x2") cc1=1;
576  if (c1=="x3") cc1=2;
577  if (c1=="x4") cc1=3;
578  if (c1=="x5") cc1=4;
579  if (c1=="x6") cc1=5;
580  if (c2=="x1") cc2=0;
581  if (c2=="x2") cc2=1;
582  if (c2=="x3") cc2=2;
583  if (c2=="x4") cc2=3;
584  if (c2=="x5") cc2=4;
585  if (c2=="x6") cc2=5;
586 
587  // output file name for the gnuplot data
588  std::string filename="plane_"+c1+c2+"_"+name;
589  // assume a cubic cell
590  double lo=-FunctionDefaults<NDIM>::get_cell_width()[0]*0.5;
591  lo=lo*scale;
592 
593  const double stepsize=FunctionDefaults<NDIM>::get_cell_width()[0]*scale/npoints;
594 
595  if(world.rank() == 0) {
596 
597  // plot 3d plot
598  FILE *f = 0;
599  f=fopen(filename.c_str(), "w");
600  if(!f) MADNESS_EXCEPTION("plot_along: failed to open the plot file", 0);
601 
602  for (int i0=0; i0<npoints; i0++) {
603  for (int i1=0; i1<npoints; i1++) {
604  // plot plane
605  coord[cc1]=lo+origin[cc1]+i0*stepsize;
606  coord[cc2]=lo+origin[cc2]+i1*stepsize;
607 
608  // other electron
609  fprintf(f,"%12.6f %12.6f %12.6f\n",coord[cc1],coord[cc2],
610  function(coord));
611 
612  }
613  // additional blank line between blocks for gnuplot
614  if (output_type=="gnuplot") fprintf(f,"\n");
615  }
616  fclose(f);
617 
618  }
619 
620  // plot mra structure
621  filename="mra_structure_"+c1+c2+"_"+name;
622  function.get_impl()->print_plane(filename.c_str(),cc1,cc2,coord);
623  }
624 
625  template<typename T>
626  static std::string stringify(T arg) {
627  std::ostringstream o;
628  if (!(o << arg))
629  MADNESS_EXCEPTION("stringify<T> failed",1);
630  return o.str();
631  }
632 
633 
634  typedef Vector<double,3> coord_3d;
635  typedef Vector<double,6> coord_6d;
636 
637  // plot along this trajectory
638  template<size_t NDIM>
639  struct trajectory {
640 
641  typedef Vector<double,NDIM> coordT;
642 
643  double lo;
644  double hi;
645  double radius;
646  long npt;
647  coordT start, end;
648  coord_3d el2;
649  coordT (*curve)(const coordT& lo, const coordT& hi, double radius, coord_3d el2, long npt, long ipt);
650 
652 
653  // return a hue number [0,0.7] according to the rank in relation to maxrank,
654  static double hueCode(const int rank) {
655  const double maxrank=40.0;
656  double hue=0.7-(0.7/maxrank)*(rank);
657  return std::max(0.0,hue);
658  }
659 
660 
661  // print a dot of hue color at (x,y) in file f
662  static void print_psdot(FILE *f, double x, double y, double color) {
663  fprintf(f,"\\newhsbcolor{mycolor}{%8.4f 1.0 0.7}\n",color);
664  fprintf(f,"\\psdot[linecolor=mycolor](%12.8f,%12.8f)\n",x,y);
665  }
666 
667 
668  static coord_3d circle2(double lo, double hi, double radius, coord_3d el2, long npt, long ipt) {
669  double stepsize=constants::pi * 2.0 / npt;
670  double phi=ipt*stepsize;
671 
672  // in the xz plane
673  coord_3d coord(0.0);
674  coord[0]=radius * sin(phi);
675  coord[1]=radius * cos(phi);
676  return coord;
677 
678  }
679 
680  static coord_6d circle_6d(const coord_6d& lo, const coord_6d& hi, double radius, coord_3d el2, long npt, long ipt) {
681  double stepsize=constants::pi * 2.0 / npt;
682 
683  // start at phi=1.0
684  double phi=1.0+constants::pi+ipt*stepsize;
685 
686  // in the xz plane
687  coord_6d coord(0.0);
688  coord[0]=radius * sin(phi);
689  coord[1]=radius * cos(phi);
690  coord[2]=0.0;
691  coord[3]=el2[0];
692  coord[4]=el2[1];
693  coord[5]=el2[2];
694 
695  return coord;
696 
697  }
698 
699 
700  // typedef Vector<double,NDIM> (trajectory::circle_6d)(double lo, double hi, double radius, long npt, long ipt) const;
701 
702  trajectory() {}
703 // // ctor for a straight line thru the origin
704 // trajectory(double lo, double hi, long npt) : lo(lo), hi(hi), npt(npt), curve(line) {
705 // }
706 
707  // ctor for circle
708  trajectory(double radius, long npt) : radius(radius), npt(npt), curve(this->circle2) {
709  }
710 
711  // ctor for circle with electron 2 fixed at coord_3d
712  trajectory(double radius, coord_3d el2, long npt) : radius(radius), npt(npt), el2(el2), curve(this->circle_6d) {
713  }
714 
715 
716  static Vector<double, NDIM> line_internal(const coordT& lo, const coordT& hi, double radius, coord_3d el2, long npt, long ipt) {
717  const coordT step=(hi-lo)*(1.0/npt);
718  coordT coord=lo+step*ipt;
719  return coord;
720  }
721 
722 
724 // static trajectory line(const Vector<double,NDIM>& lo, const Vector<double,NDIM>& hi, const long npt) {
725  static trajectory line2(const coordT start, const coordT end, const long npt) {
726  trajectory<NDIM> traj;
727  traj.start=start;
728  traj.end=end;
729  traj.npt=npt;
730  traj.curve=(trajectory::line_internal);
731  return traj;
732  }
733 
735  static trajectory line_xyz(const int xyz, const long npt) {
736  double L=FunctionDefaults<NDIM>::get_cell_width()[0];
737  coordT lo(0.0), hi(0.0);
738  lo[xyz]=-L/2;
739  hi[xyz]=L/2;
740  return trajectory<NDIM>::line2(lo,hi,npt);
741  }
742 
743  Vector<double,NDIM> operator()(int ipt) {
744  return curve(start,end,radius,el2,npt,ipt);
745  }
746 
747  };
748 
749 
750 
751  // plot along a line
752  template<size_t NDIM>
753  void plot_along(World& world, trajectory<NDIM> traj, const Function<double,NDIM>& function, std::string filename) {
754  FILE *f=0;
755  const int npt=traj.npt;
756 
757  const bool psdot=false;
758 
759  if(world.rank() == 0) {
760  f = fopen(filename.c_str(), "w");
761  if(!f) MADNESS_EXCEPTION("plot_along: failed to open the plot file", 0);
762 
763  if (psdot) {
764 
765  fprintf(f,"\\psset{xunit=0.1cm}\n");
766  fprintf(f,"\\psset{yunit=10cm}\n");
767  fprintf(f,"\\begin{pspicture}(0,-0.3)(100,1.0)\n");
768  fprintf(f,"\\pslinewidth=0.05pt\n");
769  }
770 
771  // walk along the line
772  for (int ipt=0; ipt<npt; ipt++) {
773  Vector<double,NDIM> coord=traj(ipt);
774  if (psdot) {
775  long rank=function.evalR(coord);
776  trajectory<NDIM>::print_psdot(f,ipt,function(coord),trajectory<NDIM>::hueCode(rank));
777  } else {
778  fprintf(f,"%4i %12.6f\n",ipt, function(coord));
779  }
780  }
781 
782 
783  if (psdot) fprintf(f,"\\end{pspicture}\n");
784 
785  fclose(f);
786  }
787  world.gop.fence();
788  ;
789  }
790 
791 
792  // plot along a line
793  template<size_t NDIM>
794  void plot_along(World& world, trajectory<NDIM> traj, double (*ff)(const Vector<double,NDIM>&), std::string filename) {
795  FILE *f=0;
796  const int npt=traj.npt;
797 
798  const bool psdot=false;
799 
800  if(world.rank() == 0) {
801  f = fopen(filename.c_str(), "w");
802  if(!f) MADNESS_EXCEPTION("plotvtk: failed to open the plot file", 0);
803 
804  if (psdot) {
805  fprintf(f,"\\psset{xunit=0.05cm}\n");
806  fprintf(f,"\\psset{yunit=100cm}\n");
807  fprintf(f,"\\begin{pspicture}(0,0.25)(100,0.3)\n");
808  fprintf(f,"\\pslinewidth=0.005pt\n");
809  }
810 
811 
812  // walk along the line
813  for (int ipt=0; ipt<npt; ipt++) {
814 
815  Vector<double,NDIM> coord=traj(ipt);
816 // fprintf(f,"%12.6f %12.6f %12.6f %12.6f %12.6f %12.6f %12.6f\n",coord[0],coord[1],coord[2],coord[3],coord[4],coord[5], ff(coord));
817  if (psdot) {
818  // no hue code here
819  // long rank=ff.evalR(coord);
820  trajectory<NDIM>::print_psdot(f,ipt,ff(coord),trajectory<NDIM>::hueCode(0));
821  } else {
822  fprintf(f,"%4i %12.6f\n",ipt, ff(coord));
823  }
824 
825 
826  }
827  if (psdot) fprintf(f,"\\end{pspicture}\n");
828 
829  fclose(f);
830  }
831  world.gop.fence();
832  ;
833  }
834 
835 
836 
837 }
838 
839 /* @} */
840 #endif // MADNESS_MRA_FUNCPLOT_H__INCLUDED
madness::World::gop
WorldGopInterface & gop
Global operations.
Definition: worldfwd.h:462
madness::trajectory::trajectory
trajectory(double radius, long npt)
Definition: funcplot.h:708
madness::g
NDIM const Function< R, NDIM > & g
Definition: mra.h:2179
madness::constants::pi
const double pi
Mathematical constant pi.
Definition: constants.h:44
double_complex
std::complex< double > double_complex
Definition: lineplot.cc:16
madness::imag
double imag(double x)
Definition: complexfun.h:56
madness::trajectory::radius
double radius
Definition: funcplot.h:645
madness::trajectory::el2
coord_3d el2
Definition: funcplot.h:648
NDIM
const int NDIM
Definition: tdse1.cc:44
madness::trajectory
Definition: funcplot.h:639
madness::Function::eval
Future< T > eval(const coordT &xuser) const
Evaluates the function at a point in user coordinates. Possible non-blocking comm.
Definition: mra.h:185
madness::position_stream
std::istream & position_stream(std::istream &f, const std::string &tag)
Definition: position_stream.cc:37
madness::trajectory::trajectory
trajectory(double radius, coord_3d el2, long npt)
Definition: funcplot.h:712
L
const double L
Definition: 3dharmonic.cc:123
madness::plotdx
void plotdx(const Function< T, NDIM > &f, const char *filename, const Tensor< double > &cell=FunctionDefaults< NDIM >::get_cell(), const std::vector< long > &npt=std::vector< long >(NDIM, 201L), bool binary=true)
Writes an OpenDX format file with a cube/slice of points on a uniform grid.
Definition: mraimpl.h:3228
testing::internal::string
::std::string string
Definition: gtest-port.h:872
constants.h
Defines common mathematical and physical constants.
mpfr::cos
const mpreal cos(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2255
madness::LowDimIndexIterator
Definition: indexit.h:180
madness::coord_3d
Vector< double, 3 > coord_3d
Definition: funcplot.h:634
madness::f
NDIM & f
Definition: mra.h:2179
madness::World::size
ProcessID size() const
Returns the number of processes in this world (same as MPI_Comm_size())
Definition: worldfwd.h:533
madness::trajectory::operator()
Vector< double, NDIM > operator()(int ipt)
Definition: funcplot.h:743
madness::coordT
Vector< double, 3 > coordT
Definition: chem/corepotential.cc:57
madness::plot_plane
void plot_plane(World &world, const Function< double, NDIM > &function, const std::string name)
Definition: funcplot.h:531
madness::arg
Tensor< typename Tensor< T >::scalar_type > arg(const Tensor< T > &t)
Return a new tensor holding the argument of each element of t (complex types only) ...
Definition: tensor.h:2429
madness::trajectory::lo
double lo
Definition: funcplot.h:643
madness::trajectory::start
coordT start
Definition: funcplot.h:647
madness::trajectory::hueCode
static double hueCode(const int rank)
some tools for plotting MRA ranks of low order tensors
Definition: funcplot.h:654
madness::trajectory::circle2
static coord_3d circle2(double lo, double hi, double radius, coord_3d el2, long npt, long ipt)
Definition: funcplot.h:668
madness::plot_along
void plot_along(World &world, trajectory< NDIM > traj, const Function< double, NDIM > &function, std::string filename)
Definition: funcplot.h:753
madness::trajectory::coordT
Vector< double, NDIM > coordT
Definition: funcplot.h:641
madness::Function::world
World & world() const
Returns the world.
Definition: mra.h:622
madness::trajectory::line_internal
static Vector< double, NDIM > line_internal(const coordT &lo, const coordT &hi, double radius, coord_3d el2, long npt, long ipt)
Definition: funcplot.h:716
std::tr1::T
const T1 &f1 return GTEST_2_TUPLE_() T(f0, f1)
max
#define max(a, b)
Definition: lda.h:53
madness::scale
void scale(World &world, std::vector< Function< T, NDIM > > &v, const std::vector< Q > &factors, bool fence=true)
Scales inplace a vector of functions by distinct values.
Definition: vmra.h:290
madness::Function::reconstruct
void reconstruct(bool fence=true) const
Reconstructs the function, transforming into scaling function basis. Possible non-blocking comm...
Definition: mra.h:737
a
FLOAT a(int j, FLOAT z)
Definition: y1.cc:86
madness::trajectory::curve
coordT(* curve)(const coordT &lo, const coordT &hi, double radius, coord_3d el2, long npt, long ipt)
Definition: funcplot.h:649
madness::trajectory::print_psdot
static void print_psdot(FILE *f, double x, double y, double color)
Definition: funcplot.h:662
sqrt
tensorT sqrt(const tensorT &s, double tol=1e-8)
Computes matrix square root (not used any more?)
Definition: DFcode/moldft.cc:446
madness::plotvtk_begin
void plotvtk_begin(World &world, const char *filename, const Vector< double, NDIM > &plotlo, const Vector< double, NDIM > &plothi, const Vector< long, NDIM > &npt, bool binary=false)
Definition: funcplot.h:97
madness::World
A parallel world with full functionality wrapping an MPI communicator.
Definition: worldfwd.h:416
madness::WorldGopInterface::fence
void fence()
Synchronizes all processes in communicator AND globally ensures no pending AM or tasks.
Definition: worldgop.cc:52
madness::Function
A multiresolution adaptive numerical function.
Definition: derivative.h:61
madness::coord_6d
Vector< double, 6 > coord_6d
Definition: funcplot.h:635
madness::WorldGopInterface::barrier
void barrier()
Synchronizes all processes in communicator ... does NOT fence pending AM or tasks.
Definition: worldgop.h:657
mpfr::sum
const mpreal sum(const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode)
Definition: mpreal.cc:241
mpfr::sin
const mpreal sin(const mpreal &v, mp_rnd_t rnd_mode)
Definition: mpreal.h:2262
madness::World::rank
ProcessID rank() const
Returns the process rank in this world (same as MPI_Comm_rank()))
Definition: worldfwd.h:526
madness::real
double real(double x)
Definition: complexfun.h:52
madness::trajectory::line_xyz
static trajectory line_xyz(const int xyz, const long npt)
EZ ctor for a line a direction xyz={0,1,2,..,NDIM-1} thru the origin.
Definition: funcplot.h:735
madness::plot_line
void plot_line(const char *filename, int npt, const Vector< double, NDIM > &lo, const Vector< double, NDIM > &hi, const Function< T, NDIM > &f)
Generates ASCII file tabulating f(r) at npoints along line r=lo,...,hi.
Definition: funcplot.h:388
madness::trajectory::line2
static trajectory line2(const coordT start, const coordT end, const long npt)
constructor for a line
Definition: funcplot.h:725
madness::FunctionDefaults::get_cell_width
static const Tensor< double > & get_cell_width()
Returns the width of each user cell dimension.
Definition: funcdefaults.h:391
madness::trajectory::circle_6d
static coord_6d circle_6d(const coord_6d &lo, const coord_6d &hi, double radius, coord_3d el2, long npt, long ipt)
Definition: funcplot.h:680
madness::trajectory::npt
long npt
Definition: funcplot.h:646
MADNESS_EXCEPTION
#define MADNESS_EXCEPTION(msg, value)
Definition: worldexc.h:88
madness::detail::htons_x
unsigned short htons_x(unsigned short a)
Definition: funcplot.h:322
madness::plotvtk_end
void plotvtk_end(World &world, const char *filename, bool binary=false)
Definition: funcplot.h:300
madness::trajectory::end
coordT end
Definition: funcplot.h:647
madness::trajectory::hi
double hi
Definition: funcplot.h:644
madness
Holds machinery to set up Functions/FuncImpls using various Factories and Interfaces.
Definition: chem/atomutil.cc:45
b
FLOAT b(int j, FLOAT z)
Definition: y1.cc:79
madness::trajectory::trajectory
trajectory()
Definition: funcplot.h:702
madness::plotvtk_data
void plotvtk_data(const T &function, const char *fieldname, World &world, const char *filename, const Vector< double, NDIM > &plotlo, const Vector< double, NDIM > &plothi, const Vector< long, NDIM > &npt, bool binary=false)
Definition: funcplot.h:176
PROFILE_FUNC
#define PROFILE_FUNC
Definition: worldprofile.h:198
Generated on Fri Sep 26 2014 10:38:50 for MADNESS by   doxygen 1.8.8