#include "config.h"
#include "graph_drawer.h"
#include "point_drawer.h"
#include "exc.H"
#include "terminal.h"
#include "plottable.h"
#include "fgraph.h"
#include "axis.h"
#include "frame.h"
#include "warning.h"
#include "ignore.h"
#include "color_legend.h"
#include "mframe.h"
#include "global.h"
#include "blop_gts.h"
#include <iostream>
#include <cstdio>
#include <map>
#include <limits>
#include <stdlib.h>

using namespace std;

namespace blop
{
    inline double sq(double a) { return a*a; }

    double linear(double x,double x1,double y1,double x2,double y2)
    {
	return y1 + (y2-y1)/(x2-x1)*(x-x1);
    }

    color color_mapping::the_color_;

    color_mapping::color_mapping()
    {
	p2f_ = def;
	p2f_cint_ = 0;
	fname_ = "color_mapping::def";
	start_color_ = white;
	end_color_ =  black;
    }

    void color_mapping::set(const color &start_color, const color &end_color)
    {
	p2f_ = 0;
	p2f_cint_ = 0;
	fname_ = "this_should_not_happen";
	start_color_ = start_color;
	end_color_ = end_color;
    }

    color color_mapping::def(double value,double mini,double maxi)
    {
	double norm = (value - mini) / (maxi - mini);
	if(norm < 0)
	{
	    warning::print("Too small value","color_mapping::def(...)");
	    norm = 0;
	}
	if(norm > 1)
	{
	    warning::print("Too large value","color_mapping::def");
	    norm = 1;
	}

	double red = 0,green = 0,blue = 0;

	if(norm<0.33)
	{
	    blue = linear(norm,0,0.3,0.33,1);
	}
	if(0.33<=norm && norm<0.8)
	{
	    blue = linear(norm,0.33,1,0.8,0);
	}

	
	if(0.15 < norm && norm < 0.8)
	{
	    red = linear(norm,0.15,0,0.8,1);
	}
	if(norm >= 0.8) red = 1;

	if(norm > 0.66)
	{
	    green = linear(norm,0.66,0,1,1);
	}
	
	return color(red,green,blue);
    }


    color color_mapping::map(double val, double mini, double maxi) const
    {
	if(p2f_) return p2f_(val,mini,maxi);
	else if(p2f_cint_) return map_interpreted(val,mini,maxi);

	double f = (val-mini)/(maxi-mini);
	color result;
	result.red(start_color_.red() + f * (end_color_.red() - start_color_.red()));
	result.green(start_color_.green() + f * (end_color_.green() - start_color_.green()));
	result.blue(start_color_.blue() + f * (end_color_.blue() - start_color_.blue()));
	return result;

	//return def(val, mini, maxi);
    }


    function graph_drawer::get_x(const plottable *g) const
    {
	function result = g->x_hint();
	if(!result.initialized()) result = _1;
	return result;
    }
    function graph_drawer::get_y(const plottable *g) const
    {
	function result = g->y_hint();
	if(!result.initialized()) result = _2;
	return result;
    }

    void set_ranges_default(plottable *g,axis *xaxis,axis *yaxis,int i1,int i2)
	TRY
    {
	if(!g) err("Plottable to draw not set");
	if(!xaxis) err("Xaxis not set");
	if(!yaxis) err("Yaxis not set");


	if(g->empty())
	{                                                                       
	    warning::print("Warning: empty plottable!",
			   "void set_ranges_default(plottable*, axis*, axis*, int, int)");
	    return;                                                             
	}

	if(g->xmin()!=unset) xaxis->extend_range(g->xmin());
	if(g->xmax()!=unset) xaxis->extend_range(g->xmax());
	if(g->ymin()!=unset) yaxis->extend_range(g->ymin());
	if(g->ymax()!=unset) yaxis->extend_range(g->ymax());

	// if all ranges are defined, then simply return at this point
	if(g->xmin()!=unset && g->xmax()!=unset &&
	   g->ymin()!=unset && g->ymax()!=unset) return;

	double xmin=unset;
	double xmax=unset;
	double ymin=unset;
	double ymax=unset;
	function getx=g->x_hint();
	if (!getx.initialized()) getx=ARG(i1);
	function gety=g->y_hint();
	if (!gety.initialized()) gety=ARG(i2);

	for ( unsigned int i=0 ; i<g->size() ; ++i )
	{
	    const var xorig_var=getx.eval(*g->get(i));
	    if( ignore::it(xorig_var) ) continue;
	    const double xorig=xorig_var.dbl();
	    if( xaxis->logscale() && xorig<=0.0 ) continue;
	    if( g->xmin() != unset && xorig<g->xmin() ) continue;
	    if( g->xmax() != unset && xorig>g->xmax() ) continue;
	    if(xmin == unset || xorig<xmin) xmin = xorig;
	    if(xmax == unset || xorig>xmax) xmax = xorig;

	    const var yorig_var=gety.eval(*g->get(i));
	    if ( ignore::it(yorig_var) ) continue;
	    const double yorig=yorig_var.dbl();
	    if( yaxis->logscale() && yorig<=0.0 ) continue;
	    if( g->ymin()!=unset && yorig<g->ymin()) continue;
	    if( g->ymax()!=unset && yorig>g->ymax()) continue;
	    if(ymin == unset || yorig < ymin) ymin = yorig;
	    if(ymax == unset || yorig > ymax) ymax = yorig;

	}

	if (g->xmin()==unset && xmin!=unset) xaxis->extend_range(xmin);
	if (g->xmax()==unset && xmax!=unset) xaxis->extend_range(xmax);
	if (g->ymin()==unset && ymin!=unset) yaxis->extend_range(ymin);
	if (g->ymax()==unset && ymax!=unset) yaxis->extend_range(ymax);	
    }
    CATCH("set_ranges_default(...)")

    // ---------- graphd_dots ---------------------------------------------------------

    void dots::set_ranges(plottable *g, axis *x, axis *y)
    {
	set_ranges_default(g,x,y,1,2);
    }

    void dots::draw(plottable *g,frame *f,terminal *term)
    {
	if(!g) err("Plottable to draw not set");
	if(g->empty()) return;

	function getx = g->x_hint();
	function gety = g->y_hint();
	if(!getx.initialized()) getx = _1;
	if(!gety.initialized()) gety = _2;
	
	term->set_color(g->pointcolor());
	term->reset_transformation();

	axis *xaxis =
	    (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis =
	    (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());

	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    if(ignore::it(getx.eval(*(g->get(i))))) continue;
	    if(ignore::it(gety.eval(*(g->get(i))))) continue;

	    double xorig = getx.eval(*(g->get(i))).dbl();
	    if ( xaxis->logscale() && !(xorig>0.0) ) continue;
	    if(g->xmin() != unset &&  xorig < g->xmin()) continue;
	    if(g->xmax() != unset &&  xorig > g->xmax()) continue;
	    double x = xaxis->map_point(xorig);
	    if(x<0 || 1<x) continue;

	    double yorig = gety.eval(*(g->get(i))).dbl();
	    if ( yaxis->logscale() && !(yorig>0.0) ) continue;
	    if(g->ymin() != unset && yorig < g->ymin()) continue;
	    if(g->ymax() != unset && yorig > g->ymax()) continue;
	    double y = yaxis->map_point(yorig);
	    if(y<0 || 1<y) continue;

	    term->draw_dot(terminal::coord(terminal::id(x,1),terminal::id(y,2)));
	}
    }

    void dots::draw_sample(const length &x,
			   const length &y,
			   const length &size,
			   const plottable *style,
			   terminal *t)
    {
	t->set_color(style->pointcolor());
	for(int i=0; i<15; ++i)
	{
	    length l1 = !x + (drand48()-0.5)*!size;
	    length l2 = !y + ((drand48()-0.5)*0.8)*EM;
	    l1.specialize(t);
	    l2.specialize(t);
	    t->draw_dot(terminal::coord(l1.termspecific_id(), l2.termspecific_id()));
	}
    }
    
    graph_drawer *dots::clone() const
    {
	return new dots;
    }

    // ---------- lines ---------------------------------------------------------------

    void lines::prepare_for_draw(plottable *g, frame *, int count)
    {
	if(count==1) g->linewidth().register_me();
    }

    void lines::set_ranges(plottable *g,axis *xaxis,axis *yaxis)
	TRY
    {
	set_ranges_default(g,xaxis,yaxis,1,2);
    }
    CATCH("lines::set_ranges(...)")

    graph_drawer *lines::clone() const
    {
	return new lines(*this);
    }

    void draw_lines_raw(const vector<double> &x, const vector<double> &y,
			double xmin_a, double xmax_a, double ymin_a, double ymax_a,
			terminal *term,bool fill, bool break_line = true)
    {
	const int n = ::min(x.size(),y.size());

	vector<terminal::coord> c;

	if(term->clip(terminal::coord(terminal::id(xmin_a,1), terminal::id(ymin_a,2)),
		      terminal::coord(terminal::id(xmax_a,1), terminal::id(ymax_a,2))))
	{
	    for(int i=0; i<n; ++i)
	    {
		if(x[i] == unset || y[i] == unset)
		{
                    // If the line should be broken at empty (unset) values, then flush the
                    // current stack of points, and start it again. Otherwise just skip this point.
                    if(break_line)
                    {
                        if(!fill && !c.empty())
                        {
                            if(c.size()>1)
                            {
                                if(fill) term->fill_polygon(c);
                                else term->draw_lines(c);
                            }
                            else
                            {
                                warning::print("Separated single point skipped when drawing a line",
                                               "draw_lines_raw(...)");
                            }
                            c.clear();
                        }
                    }
		    continue;
		}
		c.push_back(terminal::coord(terminal::id(x[i],1),terminal::id(y[i],2)));
	    }
	    if(fill) term->fill_polygon(c);
	    else term->draw_lines(c);
	    term->noclip();
	    return;
	}

	warning::print("This terminal does not support clipping. Trying own algorithm, "
		       "no guarantee for complicated filled polygons ...","draw_lines_raw(...)");

	for(int i=1; i<n; ++i)
	{
	    if(x[i] == unset ||
	       isinf(x[i]) ||
	       isnan(x[i]) ||
	       y[i] == unset ||
	       isinf(y[i]) ||
	       isnan(y[i]))
	    {
		if(!c.empty())
		{
		    if(fill) term->fill_polygon(c);
		    else term->draw_lines(c);
		    c.clear();
		    ++i;
		}
	    }
	    else
	    {
		double x2 = x[i];
		double y2 = y[i];

		double x1 = x[i-1];
		double y1 = y[i-1];

		double xf_1 = (xmin_a - x1)/(x2 - x1);
		double xf_2 = (xmax_a - x1)/(x2 - x1);
		double yf_1 = (ymin_a - y1)/(y2 - y1);
		double yf_2 = (ymax_a - y1)/(y2 - y1);
		double f_min = ::max(::min(xf_1,xf_2),::min(yf_1,yf_2));
		double f_max = ::min(::max(xf_1,xf_2),::max(yf_1,yf_2));

		if( (f_min<=0 && 1<=f_max) ||
		    (0<=f_min && f_min<=1 && 1<=f_max) ||
		    (f_min<=0 && 0<=f_max && f_max<=1) ||
		    (0<=f_min && f_min<=1 && 0<=f_max && f_max<=1) )
		{
		    const double x_in = x1+::max(f_min,0.0)*(x2-x1);
		    const double y_in = y1+::max(f_min,0.0)*(y2-y1);
		    const double x_out = x1+::min(f_max,1.0)*(x2-x1);
		    const double y_out = y1+::min(f_max,1.0)*(y2-y1);

		    if(f_min>0.00001 || c.empty())
		    {
			if(!c.empty())
			{
			    if(!fill) cerr<<"Strange ....."<<endl;
			    if(::fabs(c[c.size()-1].x.relcoord - x_in) > 0.001 &&
			       ::fabs(c[c.size()-1].y.relcoord - y_in) > 0.001)
			    {
				double x_extra=0, y_extra=0;
				if(x1 > xmax_a) x_extra = xmax_a;
				else if(x1 < xmin_a) x_extra = xmin_a;
				else cerr<<"I don't understsand this 1"<<endl;
				if(y1 > ymax_a) y_extra = ymax_a;
				else if(y1 < ymin_a) y_extra = ymin_a;
				else cerr<<"I don't understsand this 2"<<endl;
				c.push_back(terminal::coord(terminal::id(x_extra,1),
							    terminal::id(y_extra,2)));
			    }
			}
			c.push_back(terminal::coord(terminal::id(x_in,1),
						    terminal::id(y_in,2)));
		    }
		    c.push_back(terminal::coord(terminal::id(x_out,1),
						terminal::id(y_out,2)));
		    if(f_max<0.999999)
		    {
			if(!fill)
			{
			    term->draw_lines(c);
			    c.clear();
			}
		    }
		}
	    }
	}
	if(!c.empty())
	{
	    if(fill) term->fill_polygon(c);
	    else term->draw_lines(c);
	}

    }

    void lines::draw(plottable *g,frame *f,terminal *term)
	TRY
    {
	if(!g) err("Plottable to draw not set");
	if(g->size() < 2)  return;

	function getx = g->x_hint();
	function gety = g->y_hint();
	if(!getx.initialized()) getx = _1;
	if(!gety.initialized()) gety = _2;

	axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());

	double xmin = xaxis->min();
	if(g->xmin() != unset && g->xmin() > xmin) xmin = g->xmin();
	double xmax = xaxis->max();
	if(g->xmax() != unset && g->xmax() < xmax) xmax = g->xmax();

	double ymin = yaxis->min();
	if(g->ymin() != unset && g->ymin() > ymin) ymin = g->ymin();
	double ymax = yaxis->max();
	if(g->ymax() != unset && g->ymax() < ymax) ymax = g->ymax();

	vector<double> xx(g->size());
	vector<double> yy(g->size());
	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    const var xorig = getx.eval(*(g->get(i)));
	    const var yorig = gety.eval(*(g->get(i)));
	    if(ignore::it(xorig) || (xaxis->logscale() && !(xorig.dbl()>0.0))) xx[i] = unset;
	    else xx[i] = xaxis->map_point(xorig.dbl());
	    if(ignore::it(yorig) || (yaxis->logscale() && !(yorig.dbl()>0.0))) yy[i] = unset;
	    else yy[i] = yaxis->map_point(yorig.dbl());
	}	

	term->reset_transformation();
	if(g->fill())
	{
	    term->set_color(g->fillcolor());
	    draw_lines_raw(xx, yy,
			   xaxis->map_point(xmin), xaxis->map_point(xmax),
			   yaxis->map_point(ymin), yaxis->map_point(ymax), term, true, break_line_);
	}

	term->set_color(g->linecolor());
	term->set_linestyle(g->linestyle());
	term->set_linewidth(g->linewidth().termspecific_id());
	draw_lines_raw(xx, yy,
		       xaxis->map_point(xmin), xaxis->map_point(xmax),
		       yaxis->map_point(ymin), yaxis->map_point(ymax), term, false, break_line_);

    }
    CATCH("graphd_line::draw(terminal *)")

    void lines::draw_sample(const length &x,
			    const length &y,
			    const length &size,
			    const plottable *s,
			    terminal *t)
    {
	t->set_linewidth(s->linewidth().termspecific_id());
	t->set_color(s->linecolor());
	t->set_linestyle(s->linestyle());

	length l1 = !x - 0.5*!size;
	length l2 = !x + 0.5*!size;

	l1.specialize(t);
	l2.specialize(t);

	t->draw_line(terminal::coord(l1.termspecific_id(), y.termspecific_id()),
		     terminal::coord(l2.termspecific_id(), y.termspecific_id()));
    }

    lines lines;

    // ---------- splines -------------------------------------------------------------

    graph_drawer *splines::clone() const
    {
	return new splines(*this);
    }

    void splines::draw_sample(const length &x,const length &y,
			      const length &size,
			      const plottable *g,terminal *t)
    {
	if(g->fill())
	{
	    length x1 = -0.5*!size;
	    length x2 =  0.5*!size;
	    length y1 = -0.5*EX;
	    length y2 = 0.5*EX;
	    
	    x1.specialize(t);
	    x2.specialize(t);
	    y1.specialize(t);
	    y2.specialize(t);
	    
	    vector<terminal::coord> c;
	    
	    c.push_back(terminal::coord(x1.termspecific_id(),y1.termspecific_id()));
	    c.push_back(terminal::coord(x2.termspecific_id(),y1.termspecific_id()));
	    c.push_back(terminal::coord(x2.termspecific_id(),y2.termspecific_id()));
	    c.push_back(terminal::coord(x1.termspecific_id(),y2.termspecific_id()));
	    c.push_back(terminal::coord(x1.termspecific_id(),y1.termspecific_id()));
	    t->translate(x.termspecific_id(),y.termspecific_id());
	    t->set_color(g->fillcolor());
	    t->fill_polygon(c);
	    t->reset_transformation();
	}
	else
	{
	    lines::draw_sample(x,y,size,g,t);
	}
    }


    void splines::draw(plottable *g,frame *f,terminal *term)
    {
	if(!g) err("Plottable to draw not set");
	if(g->size() < 2)
	{
	    warning::print("Can't draw plottable with 1 points","splines::draw(...)");
	    return;
	}
	if(g->size() < 3)
	{
	    warning::print("Plottable with 2 points is drawn with a line","splines::draw(...)");
	    lines::draw(g,f,term);
	    return;
	}

	function getx = g->x_hint();
	function gety = g->y_hint();
	if(!getx.initialized()) getx = _1;
	if(!gety.initialized()) gety = _2;

	if(g->fill()) term->set_color(g->fillcolor());
	else
	{
	    term->set_color(g->linecolor());
	    term->set_linestyle(g->linestyle());
	    term->set_linewidth(g->linewidth().termspecific_id());
	}
	term->reset_transformation();

	axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());

	double xmin = xaxis->min();
	if(g->xmin() != unset && g->xmin() > xmin) xmin = g->xmin();
	double xmax = xaxis->max();
	if(g->xmax() != unset && g->xmax() < xmax) xmax = g->xmax();

	double ymin = yaxis->min();
	if(g->ymin() != unset && g->ymin() > ymin) ymin = g->ymin();
	double ymax = yaxis->max();
	if(g->ymax() != unset && g->ymax() < ymax) ymax = g->ymax();

	vector<double> xa(g->size());
	vector<double> ya(g->size());

	const int n=g->size();
	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    const var xorig = getx.eval(*(g->get(i)));
	    const var yorig = gety.eval(*(g->get(i)));
	    if(ignore::it(xorig) || (xaxis->logscale() && !(xorig.dbl()>0.0))) xa[i] = unset;
	    else xa[i] = xaxis->map_point(xorig.dbl());
	    if(ignore::it(yorig) || (yaxis->logscale() && !(yorig.dbl()>0.0))) ya[i] = unset;
	    else ya[i] = yaxis->map_point(yorig.dbl());
	}
	function the_spline = function::interpolate_spline(xa,ya);

	vector<double> xx,yy;

	const double step = (xa.back()-xa.front())/100;

	double current_x = xa.front();
	if(g->fill())
	{
	    xx.push_back(current_x);
	    yy.push_back(0);
	}

	xx.push_back(current_x);
	yy.push_back(the_spline(current_x));

	for(int i=1; i<n; ++i)
	{
	    if(xa[i] == unset)
	    {
		if(!xx.empty())
		{
		    xx.push_back(xx.back());
		    yy.push_back(0);
		}
		xx.push_back(unset);
		yy.push_back(unset);
		if(i+1 < n)
		{
		    current_x = xa[i+1];
		    xx.push_back(current_x);
		    yy.push_back(0);
		}
	    }
	    else
	    {
		for(; current_x<xa[i]; current_x += step)
		{
		    if(current_x>xa[i]-0.5*step) current_x = xa[i]-0.5*step;
		    xx.push_back(current_x);
		    yy.push_back(the_spline(current_x));
		}
		current_x = xa[i];
		xx.push_back(current_x);
		yy.push_back(the_spline(current_x));
	    }
	}

	if(g->fill())
	{
	    xx.push_back(xx.back());
	    yy.push_back(0);
	}

	draw_lines_raw(xx, yy,
		       xaxis->map_point(xmin),
		       xaxis->map_point(xmax),
		       yaxis->map_point(ymin),
		       yaxis->map_point(ymax),
		       term, g->fill());

    }

    // ---------- graphd_symbol -------------------------------------------------------

    graph_drawer *points::clone() const
    {
	return new points(*this);
    }

    void points::prepare_for_draw(plottable *g, frame *, int count)
    {
	if(count==1)
	{
	    point_drawer_ = g->pointtype();
	    if(point_drawer_ == 0) point_drawer_ = new square;
	    point_drawer_->prepare_for_draw(g->pointsize());
	    g->linewidth().register_me();
	}
    }

    void points::set_ranges(plottable *g,axis *xaxis,axis *yaxis)
	TRY
    {
	set_ranges_default(g,xaxis,yaxis,1,2);
    }
    CATCH("points::set_ranges(...)")


    void points::draw(plottable *g,frame *f,terminal *term)
    TRY
    {
	if(!g) err("Plottable to draw not set");
	if(g->empty())  return;

	function getx = g->x_hint();
	function gety = g->y_hint();
	if(!getx.initialized()) getx = _1;
	if(!gety.initialized()) gety = _2;

	axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());

	double xmin = xaxis->min();
	if(g->xmin()!=unset && g->xmin()>xmin) xmin = g->xmin();
	double xmax = xaxis->max();
	if(g->xmax()!=unset && g->xmax()<xmax) xmax = g->xmax();

	double ymin = yaxis->min();
	if(g->ymin()!=unset && g->ymin()>ymin) ymin = g->ymin();
	double ymax = yaxis->max();
	if(g->ymax()!=unset && g->ymax()<ymax) ymax = g->ymax();

	term->set_color(g->pointcolor());
	term->set_linewidth(g->linewidth().termspecific_id());
	term->set_linestyle(g->linestyle());

	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    const var xorig = getx.eval(*(g->get(i))).dbl();
	    const var yorig = gety.eval(*(g->get(i))).dbl();
	    if(ignore::it(xorig) || ignore::it(yorig)) continue;
	    if(xaxis->logscale() && !(xorig.dbl()>0.0)) continue;
	    if(yaxis->logscale() && !(yorig.dbl()>0.0)) continue;
	    if(xorig.dbl() < xmin) continue;
	    if(xorig.dbl() > xmax) continue;
	    if(yorig.dbl() < ymin) continue;
	    if(yorig.dbl() > ymax) continue;

	    double x = xaxis->map_point(xorig.dbl());
	    double y = yaxis->map_point(yorig.dbl());
	    if(0<=x && x<=1 && 0<=y && y<=1)
	    {
		term->translate(terminal::id(x,1),terminal::id(y,2));
		point_drawer_->draw(term);
		term->reset_transformation();
	    }
	}
    }
    CATCH("graphd_symbol::draw(...)")

    void points::draw_sample(const length &x,
			     const length &y,
			     const length &size,
			     const plottable *s,terminal *t)
    {
	t->set_color(s->pointcolor());
	t->set_linewidth(s->linewidth().termspecific_id());
	t->set_linestyle(s->linestyle());
	t->translate(x.termspecific_id(),y.termspecific_id());
	point_drawer_->draw(t);
	t->reset_transformation();
    }


    // ---------- linespoints ---------------------------------------------------------

    graph_drawer *linespoints::clone() const
    {
	return new linespoints(*this);
    }

    void linespoints::prepare_for_draw(plottable *g, frame *, int count)
    {
	if(count==1)
	{
	    point_drawer_ = g->pointtype();
	    if(point_drawer_ == 0) point_drawer_ = new square;
	    point_drawer_->prepare_for_draw(g->pointsize());
	    g->linewidth().register_me();
	}
    }

    void linespoints::set_ranges(plottable *g,axis *xaxis,axis *yaxis)
	TRY
    {
	set_ranges_default(g,xaxis,yaxis,1,2);
    }
    CATCH("linespoints::set_ranges(...)")


    void linespoints::draw(plottable *g,frame *f,terminal *term)
    TRY
    {
	if(!g) err("Plottable to draw not set");
	if(g->empty())  return;

	function getx = g->x_hint();
	function gety = g->y_hint();
	if(!getx.initialized()) getx = _1;
	if(!gety.initialized()) gety = _2;

	axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());

	double xmin = xaxis->min();
	if(g->xmin() != unset && g->xmin() > xmin) xmin = g->xmin();
	double xmax = xaxis->max();
	if(g->xmax() != unset && g->xmax() < xmax) xmax = g->xmax();

	double ymin = yaxis->min();
	if(g->ymin() != unset && g->ymin() > ymin) ymin = g->ymin();
	double ymax = yaxis->max();
	if(g->ymax() != unset && g->ymax() < ymax) ymax = g->ymax();

	term->set_color(g->pointcolor());
	term->set_linewidth(g->linewidth().termspecific_id());
	term->set_linestyle(g->linestyle());

	vector<double> xx(g->size());
	vector<double> yy(g->size());

	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    const var xorig = getx.eval(*(g->get(i)));
	    const var yorig = gety.eval(*(g->get(i)));
	    double x,y;

	    if(ignore::it(xorig) || (xaxis->logscale() && !(xorig.dbl()>0.0))) x=unset;
	    else x = xaxis->map_point(xorig.dbl());
	    if(ignore::it(yorig) || (yaxis->logscale() && !(yorig.dbl()>0.0))) y=unset;
	    else y = yaxis->map_point(yorig.dbl());

	    xx[i] = x;
	    yy[i] = y;

	    if(xorig.dbl() < xmin) continue;
	    if(xorig.dbl() > xmax) continue;
	    if(yorig.dbl() < ymin) continue;
	    if(yorig.dbl() > ymax) continue;

	    if(0<=x && x<=1 && 0<=y && y<=1 && x!=unset && y!=unset)
	    {
		term->translate(terminal::id(x,1),terminal::id(y,2));
		point_drawer_->draw(term);
		term->reset_transformation();
	    }
	}
//	term->reset_transformation();
	draw_lines_raw(xx, yy,
		       xaxis->map_point(xmin), xaxis->map_point(xmax),
		       yaxis->map_point(ymin), yaxis->map_point(ymax), term, false);
    }
    CATCH("linespoints::draw(...)")

    void linespoints::draw_sample(const length &x,
				  const length &y,
				  const length &size,
				  const plottable *s,terminal *t)
    {
	t->set_color(s->pointcolor());
	t->set_linewidth(s->linewidth().termspecific_id());
	t->set_linestyle(s->linestyle());
	t->translate(x.termspecific_id(),y.termspecific_id());

	point_drawer_->draw(t);

	length l1 = - 0.5*!size;
	length l2 =   0.5*!size;
	l1.specialize(t);
	l2.specialize(t);
	t->draw_line(terminal::coord(l1.termspecific_id(), terminal::ZERO),
		     terminal::coord(l2.termspecific_id(), terminal::ZERO));
	t->reset_transformation();
    }

    // --------------------------------------------------------------------------------

    graph_drawer *cpoints::clone() const
    {
	return new cpoints(*this);
    }

    void cpoints::draw(plottable *g,frame *f,terminal *term)
	TRY
    {
	if(!g) err("Plottable to draw not set");
	if(g->empty())  return;

	function getx = g->x_hint();
	function gety = g->y_hint();
	function getr = _3; // red
	function getg = _4; // green
	function getb = _5; // blue

	if(!getx.initialized()) getx = _1;
	if(!gety.initialized()) gety = _2;

	axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());

	double xmin = xaxis->min();
	if(g->xmin()!=unset && g->xmin()>xmin) xmin = g->xmin();
	double xmax = xaxis->max();
	if(g->xmax()!=unset && g->xmax()<xmax) xmax = g->xmax();

	double ymin = yaxis->min();
	if(g->ymin()!=unset && g->ymin()>ymin) ymin = g->ymin();
	double ymax = yaxis->max();
	if(g->ymax()!=unset && g->ymax()<ymax) ymax = g->ymax();

	term->set_color(g->pointcolor());
	term->set_linewidth(g->linewidth().termspecific_id());
	term->set_linestyle(g->linestyle());

	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    const var xorig = getx.eval(*(g->get(i)));
	    const var yorig = gety.eval(*(g->get(i)));
	    if(ignore::it(xorig) || ignore::it(yorig)) continue;
	    if(xaxis->logscale() && !(xorig.dbl()>0.0)) continue;
	    if(yaxis->logscale() && !(yorig.dbl()>0.0)) continue;

	    if(xorig.dbl() < xmin) continue;
	    if(xorig.dbl() > xmax) continue;
	    if(yorig.dbl() < ymin) continue;
	    if(yorig.dbl() > ymax) continue;

	    double x = xaxis->map_point(xorig.dbl());
	    double y = yaxis->map_point(yorig.dbl());
	    if(0<=x && x<=1 && 0<=y && y<=1)
	    {
		double red   = getr.eval(*(g->get(i))).dbl(); 
		double green = getg.eval(*(g->get(i))).dbl(); 
		double blue  = getb.eval(*(g->get(i))).dbl(); 

		term->set_color(color(red,green,blue));
		term->translate(terminal::id(x,1),terminal::id(y,2));
		point_drawer_->draw(term);
		term->reset_transformation();
	    }
	}
    }
    CATCH("cpoints::draw(...)")

    cpoints cpoints;


    // --------------------------------------------------------------------------------

    graph_drawer *spoints::clone() const
    {
	return new spoints(*this);
    }

    void spoints::draw(plottable *g,frame *f,terminal *term)
	TRY
    {
	if(!g) err("Plottable to draw not set");
	if(g->empty())  return;

	function getx = g->x_hint();
	function gety = g->y_hint();
	function getz = g->z_hint();
	if(!getx.initialized()) getx = _1;
	if(!gety.initialized()) gety = _2;
	if(!getz.initialized()) getz = _3;

	axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());

	double xmin = xaxis->min();
	if(g->xmin()!=unset && g->xmin()>xmin) xmin = g->xmin();
	double xmax = xaxis->max();
	if(g->xmax()!=unset && g->xmax()<xmax) xmax = g->xmax();

	double ymin = yaxis->min();
	if(g->ymin()!=unset && g->ymin()>ymin) ymin = g->ymin();
	double ymax = yaxis->max();
	if(g->ymax()!=unset && g->ymax()<ymax) ymax = g->ymax();

	term->set_color(g->pointcolor());
	term->set_linewidth(g->linewidth().termspecific_id());
	term->set_linestyle(g->linestyle());

	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    const var xorig = getx.eval(*(g->get(i)));
	    const var yorig = gety.eval(*(g->get(i)));
	    const var zorig = getz.eval(*(g->get(i)));

	    if(ignore::it(xorig) ||
	       ignore::it(yorig) ||
	       ignore::it(zorig)) continue;
	    if(xaxis->logscale() && !(xorig.dbl()>0.0)) continue;
	    if(yaxis->logscale() && !(yorig.dbl()>0.0)) continue;
	    if(xorig.dbl() < xmin) continue;
	    if(xorig.dbl() > xmax) continue;
	    if(yorig.dbl() < ymin) continue;
	    if(yorig.dbl() > ymax) continue;

	    double x = xaxis->map_point(xorig.dbl());
	    double y = yaxis->map_point(yorig.dbl());
	    if(0<=x && x<=1 && 0<=y && y<=1)
	    {
		length l = !g->pointsize() * zorig.dbl();
		term->translate(terminal::id(x,1),terminal::id(y,2));
		point_drawer_->draw(term,l);
		term->reset_transformation();
	    }
	}
    }
    CATCH("spoints::draw(...)")

    spoints spoints;



    // ---------- histo ---------------------------------------------------------------

    graph_drawer *histo::clone() const
    {
	return new histo(*this);
    }

    void histo::set_ranges(plottable *g,axis *xaxis,axis *yaxis)
	TRY
    {
	if(!g) err("Plottable to draw not set");
	if(!xaxis) err("Xaxis not set");
	if(!yaxis) err("Yaxis not set");
	if(g->empty()) return;

	function getx = g->x_hint();
	function gety = g->y_hint();
	if(!getx.initialized()) getx = _1;
	if(!gety.initialized()) gety = _2;

	function getx1 = g->x1_hint();
	function getx2 = g->x2_hint();

	yaxis->autoextend_min_soft(false);


	if(g->xmin() != unset && (xaxis->logscale() == false || g->xmin()>0))
	{
	    xaxis->extend_range(g->xmin());
	}
	else
	{
	    if(g->xmin()!=unset && g->xmin() <= 0 && xaxis->logscale())
	    {
		warning::print("xmin set to negative value with logscale, ignoring xmin");
	    }

	    double xfirst=unset;
	    double xsec=unset;
	    int indfirst=0;
	    int indsec=1;
	    for ( unsigned int i=0 ; i<g->size() ; ++i )
	    {
		var xorig=getx.eval(*(g->get(i)));
		if(ignore::it(xorig) || (xaxis->logscale() && !(xorig.dbl()>0.0))) continue;
		if( xfirst==unset) { xfirst=xorig.dbl(); indfirst=i; continue; }
		xsec=xorig.dbl(); 
		indsec=i;
		break; 
	    }
	    if (getx1.initialized())
	    {
		xaxis->extend_range(getx1.eval(*(g->get(indfirst))).dbl());
	    }
	    else if ( xsec!=unset )
	    {
		xaxis->extend_range(xfirst-0.5*(xsec - xfirst)/(indsec-indfirst));
	    }
	    else if ( xfirst!=unset )
	    {
		xaxis->extend_range(xfirst- 0.5);
	    }
	    else
	    {
		warning::print("Could not set range with histogram style");
	    }
	}

	if(g->xmax() != unset && (xaxis->logscale() == false || g->xmax()>0))
	{
	    xaxis->extend_range(g->xmax());
	}
	else
	{
	    if(g->xmax()!=unset && g->xmax()<=0 && xaxis->logscale())
	    {
		warning::print("xmax set to negative value with logscale, ignoring xmax");
	    }
	    double xfirst=unset;
	    double xsec=unset;
	    int indfirst=g->size()-1;
	    int indsec=g->size()-2;
	    for ( unsigned int i=g->size()-1 ; i>=0 ; --i )
	    {
		var xorig=getx.eval(*(g->get(i)));
		if(ignore::it(xorig) || (xaxis->logscale() && !(xorig.dbl()>0.0))) continue;
		if ( xfirst==unset && xsec==unset ) { xfirst=xorig.dbl(); indfirst=i; continue; }
		if ( xfirst!=unset && xsec==unset ) { xsec=xorig.dbl(); indsec=i;  break; }
	    }
	    if (getx2.initialized())
	    {
		xaxis->extend_range(getx2.eval(*(g->get(indfirst))).dbl());
	    }
	    if ( xsec!=unset )
	    {
		xaxis->extend_range(xfirst+0.5*(xfirst - xsec)/(indfirst-indsec));
	    }
	    else if ( xfirst!=unset )
	    {
		xaxis->extend_range(xfirst + 0.5);
	    }
	    else
	    {
		warning::print("Could not set x range with histogram style");
	    }
	}

	bool ymin_set = false, ymax_set = false;
	if(g->ymin()!=unset)
	{
	    if(g->ymin()>0 || yaxis->logscale()==false)
	    {
		yaxis->extend_range(g->ymin());
		ymin_set = true;
	    }
	    else
	    {
		warning::print("ymin set to negative value with logscale y, ignoring");
	    }
	}
	if(g->ymax()!=unset)
	{
	    if(g->ymax()>0 || yaxis->logscale()==false)
	    {
		yaxis->extend_range(g->ymax());
		ymax_set = true;
	    }
	    else
	    {
		warning::print("ymax set to negative value with logscale y, ignoring");
	    }
	}
	if(ymin_set && ymax_set) return;

	double min_y = unset;
	double max_y = unset;
	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    const var xorig_var = getx.eval(*(g->get(i)));
	    const var yorig_var = gety.eval(*(g->get(i)));
	    if(ignore::it(xorig_var) || ignore::it(yorig_var)) continue;
	    const double xorig = xorig_var.dbl();
	    const double yorig = yorig_var.dbl();
	    
	    if(xaxis->logscale() && xorig<=0)        continue;
	    if(g->xmin()!=unset  && xorig<g->xmin()) continue;
	    if(g->xmax()!=unset  && xorig>g->xmax()) continue;

	    if(yaxis->logscale() && yorig<=0) continue;
	    if(min_y == unset || yorig<min_y) min_y = yorig;
	    if(max_y == unset || yorig>max_y) max_y = yorig;
	}

	// have not set ymin before yet
	if(ymin_set == false)
	{
	    if(!yaxis->logscale())
	    {
		if(min_y!=unset) yaxis->extend_range(min_y);
		yaxis->extend_range(0);
	    }
	    else if(min_y > 0 && min_y!=unset) yaxis->extend_range(min_y * 0.5);
	}

	if(ymax_set == false)
	{
	    if(max_y != unset) yaxis->extend_range(max_y);
	}

    }
    CATCH("histo::set_ranges(plottable *,axis *,axis *)")

    
    void histo::draw(plottable *g,frame *f,terminal *term)
    TRY
    {
	if(!g) err("Plottable to draw not set");
	if(g->empty())
	{
	    warning::print("Empty graph","histo::draw(...)");
	    return;
	}
	if(g->size() == 1)
	{
	    warning::print("Can't yet draw a histo with 1 entries","histo::draw(...)");
	    return;
	}

	function getx = g->x_hint();
	function gety = g->y_hint();
	if(!getx.initialized()) getx = _1;
	if(!gety.initialized()) gety = _2;

	function getx1 = g->x1_hint();
	function getx2 = g->x2_hint();

	if(g->fill()) term->set_color(g->fillcolor());
	else term->set_color(g->linecolor());

	term->set_linestyle(g->linestyle());
	term->set_linewidth(g->linewidth().termspecific_id());
	term->reset_transformation();

	axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());
	
	vector<double>  xxx,yyy;

	double xmin = xaxis->min();
	if(g->xmin() != unset && g->xmin() > xmin) xmin = g->xmin();
	double xmax = xaxis->max();
	if(g->xmax() != unset && g->xmax() < xmax) xmax = g->xmax();

	double xmin_a = xaxis->map_point(xmin);
	double xmax_a = xaxis->map_point(xmax);
	double ymin_a = 0;
	if(g->ymin() != unset) ymin_a = ::max(yaxis->map_point(g->ymin()), 0.0);
	double ymax_a = 1;
	if(g->ymax() != unset) ymax_a = ::min(yaxis->map_point(g->ymax()), 1.0);

	// the first entry:

	for(plottable::size_type index = 0; index < g->size(); ++index)
	{
	    const var xorig_var = getx.eval(*(g->get(index)));
	    const var yorig_var = gety.eval(*(g->get(index)));
	    if(ignore::it(xorig_var) || ignore::it(yorig_var)) continue;

	    double x1=0, x2=0;

	    if(getx1.initialized() && getx2.initialized())
	    {
		x1 = getx1.eval(*(g->get(index))).dbl();
		x2 = getx2.eval(*(g->get(index))).dbl();
	    }
	    else
	    {
		plottable::size_type next_index = index+1;
		if (next_index<g->size())
		{
		    var xnew=getx.eval(*(g->get(next_index)));
		    while(next_index<g->size() && (ignore::it(xnew) || (xaxis->logscale() && !(xnew.dbl()>0.0)))) { ++next_index; xnew=getx.eval(*(g->get(next_index))); }
		}

		int prev_index = index-1;
		if ( prev_index>=0)
		{
		    var xold=getx.eval(*(g->get(prev_index)));
		    while(prev_index >= 0 && (ignore::it(xold) || (xaxis->logscale() && !(xold.dbl()>0.0)))) { --prev_index; xold=getx.eval(*(g->get(prev_index))); }
		}

		if(next_index >= g->size() && prev_index < 0)
		{
		    cerr<<"Strange, please report this problem [next_index>=g->size() && prev_index<0"
			" in histo::draw(...)]"<<endl;
		    continue;
		}

		//- modifiactions begin

		/*
		  if(next_index < g->size())
		  {
		  x2 = ( getx.eval(*(g->get(index))).dbl() + getx.eval(*(g->get(next_index))).dbl() ) /2;
		  }
		  else
		  {
		  x2 = ( 1.5*getx.eval(*(g->get(index))).dbl() - 0.5*getx.eval(*(g->get(prev_index))).dbl() );
		  }

		  if(prev_index >= 0)
		  {
		  x1 = ( getx.eval(*(g->get(prev_index))).dbl() + getx.eval(*(g->get(index))).dbl() )/2;
		  }
		  else
		  {
		  x1 = ( 1.5*getx.eval(*(g->get(index))).dbl() - 0.5*getx.eval(*(g->get(next_index))).dbl() );
		  }
		*/

		if (ignore::it(xorig_var) || (xaxis->logscale() && !(xorig_var.dbl()>0.0)))
		{
		    if(next_index >= g->size()) continue;
		    else if(prev_index < 0) continue;
		    else
		    {
			double xother=getx.eval(*(g->get(prev_index))).dbl();
			double xoother=getx.eval(*(g->get(next_index))).dbl();
			x1 = xother + 0.5*( xoother - xother )/(next_index-prev_index);
			x2 = xoother - 0.5*( xoother - xother )/(next_index-prev_index);
			if ( xxx.size()>0 )
			{
			    if ( ::fabs(xxx[xxx.size()-1]-xaxis->map_point(x1))<0.000001*(xaxis->map_point(x2)-xaxis->map_point(x1)) ) continue;
			}
		    }
		}
		else
		{
		    if(next_index < g->size())
		    {
			double xother=getx.eval(*(g->get(next_index))).dbl();
			x2 = 0.5*( getx.eval(*(g->get(index))).dbl() + xother )/(next_index-index);
		    }
		    else
		    {
			double xother=0.5*getx.eval(*(g->get(prev_index))).dbl();
			x2 = ( 1.5*getx.eval(*(g->get(index))).dbl() - xother )/(index-prev_index);
		    }

		    if(prev_index >= 0)
		    {
			double xother=getx.eval(*(g->get(index))).dbl();
			x1 = 0.5*( getx.eval(*(g->get(prev_index))).dbl() + xother )/(index-prev_index);
		    }
		    else
		    {
			double xother=0.5*getx.eval(*(g->get(next_index))).dbl();
			x1 = ( 1.5*getx.eval(*(g->get(index))).dbl() - xother )/(next_index-index);
		    }
		}
	    }

	    if(x2<xmin || xmax<x1) continue;

	    double x1_inrange = x1;
	    if(x1 < xmin) x1_inrange = xmin;

	    double x2_inrange = x2;
	    if(x2 > xmax) x2_inrange = xmax;

	    double x1a = xaxis->map_point(x1_inrange);
	    double x2a = xaxis->map_point(x2_inrange);

	    double ya = 0;
	    const var yorig = gety.eval(*(g->get(index)));
	    if(!ignore::it(xorig_var) &&
	       (!xaxis->logscale() || xorig_var.dbl()>0.0) &&
	       !ignore::it(yorig) &&
	       (!yaxis->logscale() || yorig.dbl()>0))
	    {
		ya = yaxis->map_point(yorig.dbl());
	    }
	    else
	    {
		ya = unset;
	    }

	    // if first bin's low edge fully contained in the xaxis range,
	    // add an y=0 value
	    if(xxx.empty() && (x1 == x1_inrange || g->fill()))
	    {
		xxx.push_back(x1a);
		double map0;
		if(!yaxis->logscale()) map0 = yaxis->map_point(0);
		else map0 = 0;
		yyy.push_back(map0);
	    }

	    if(ya != unset)
	    {
		xxx.push_back(x1a);
		yyy.push_back(ya);
		xxx.push_back(x2a);
		yyy.push_back(ya);
	    }
	    else
	    {
		if(g->fill())
		{
		    xxx.push_back(x1a);
		    yyy.push_back(0);
		}
		xxx.push_back(unset);
		yyy.push_back(unset);
		if(g->fill())
		{
		    xxx.push_back(x2a);
		    yyy.push_back(0);
		}
	    }

	    if( (index == g->size()-1 && x2<=xmax) || // last point fully contained in displayed range
		(xmax <= x2 && g->fill()))            // next point will be outside
	    {
		xxx.push_back(x2a);
		double map0;
		if(!yaxis->logscale()) map0 = yaxis->map_point(0);
		else map0 = 0;
		yyy.push_back(map0);
	    }

	}

	draw_lines_raw(xxx, yyy, xmin_a, xmax_a, ymin_a, ymax_a, term, g->fill());

    }
    CATCH("histo::draw(plottable *,frame *,terminal *)")


    // ---------- bars ----------------------------------------------------------------


    bars::bars(const length &width, const length &offset)
    {
	l1_ = offset - 0.5*width;
	l2_ = offset + 0.5*width;
    }


    bars &bars::operator() (const length &width, const length &offset=0.0)
    {
	l1_ = offset - 0.5*width;
	l2_ = offset + 0.5*width;
	return *this;
    }

    void bars::set_ranges(plottable *g,axis *x,axis *y)
    {
	set_ranges_default(g,x,y,1,2);
	if(!y->logscale()) y->extend_range(0);
    }

    void bars::draw_sample(const length &x,
			   const length &y,
			   const length &size,
			   const plottable *style,
			   terminal *t)
    {
	length x1 = -0.5*!size;
	length x2 =  0.5*!size;
	length y1 = -0.5*EX;
	length y2 = 0.5*EX;

	x1.specialize(t);
	x2.specialize(t);
	y1.specialize(t);
	y2.specialize(t);
	
	vector<terminal::coord> c;

	c.push_back(terminal::coord(x1.termspecific_id(),y1.termspecific_id()));
	c.push_back(terminal::coord(x2.termspecific_id(),y1.termspecific_id()));
	c.push_back(terminal::coord(x2.termspecific_id(),y2.termspecific_id()));
	c.push_back(terminal::coord(x1.termspecific_id(),y2.termspecific_id()));
	c.push_back(terminal::coord(x1.termspecific_id(),y1.termspecific_id()));
	t->translate(x.termspecific_id(),y.termspecific_id());
	t->set_color(style->fillcolor());
	t->fill_polygon(c);
	t->reset_transformation();
    }

    void bars::draw(plottable *g,frame *f,terminal *term)
    {
	if(!g) err("Plottable to draw not set");
	if(g->empty()) return;

	function getx = g->x_hint();
	function gety = g->y_hint();
	if(!getx.initialized()) getx = _1;
	if(!gety.initialized()) gety = _2;

	axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());

	double xmin = xaxis->min();
	if(g->xmin()!=unset && g->xmin()>xmin) xmin = g->xmin();
	double xmax = xaxis->max();
	if(g->xmax()!=unset && g->xmax()<xmax) xmax = g->xmax();

	double ymin = yaxis->min();
	if(g->ymin()!=unset && g->ymin()>ymin) ymin = g->ymin();
	double ymax = yaxis->max();
	if(g->ymax()!=unset && g->ymax()<ymax) ymax = g->ymax();

	term->set_color(g->fillcolor());
	vector<terminal::coord> c;
	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    const var xorig_var = getx.eval(*(g->get(i)));
	    const var yorig_var = gety.eval(*(g->get(i)));
	    if(ignore::it(xorig_var) || ignore::it(yorig_var)) continue;
	    const double xorig = xorig_var.dbl();
	    const double yorig = yorig_var.dbl();
	    if ( xaxis->logscale() && xorig<=0.0 ) continue;
	    if ( yaxis->logscale() && yorig<=0.0 ) continue;
	    if(xorig < xmin) continue;
	    if(xorig > xmax) continue;
	    if(yorig < ymin) continue;
	    if(yorig > ymax) continue;

	    double x = xaxis->map_point(xorig);
	    double y2 = yaxis->map_point(yorig);
	    double y1 = yaxis->map_point(0);

	    term->translate(terminal::id(x,1),0);
	    c.clear();
	    c.push_back(terminal::coord(l1_.termspecific_id(),terminal::id(y1,2)));
	    c.push_back(terminal::coord(l2_.termspecific_id(),terminal::id(y1,2)));
	    c.push_back(terminal::coord(l2_.termspecific_id(),terminal::id(y2,2)));
	    c.push_back(terminal::coord(l1_.termspecific_id(),terminal::id(y2,2)));
	    term->fill_polygon(c);
	    term->reset_transformation();
	}
    }

    graph_drawer *bars::clone() const
    {
	return new bars(*this);
    }

    void bars::prepare_for_draw(plottable *, frame *, int count)
    {
	if(count==1)
	{
	    l1_.register_me();
	    l2_.register_me();
	}
    }

    // ---------- labels --------------------------------------------------------------

    labels &labels::xalign(sym::position a)
    {
	xalign_ = a;
	return *this;
    }
    labels &labels::yalign(sym::position a)
    {
	yalign_ = a;
	return *this;
    }
    labels &labels::align(sym::position xa, sym::position ya)
    {
	xalign_ = xa;
	yalign_ = ya;
	return *this;
    }
    labels &labels::angle(double a)
    {
	angle_ = a;
	return *this;
    }
    labels &labels::xoffset(const length &l)
    {
	x_offset_ = l;
	return *this;
    }
    labels &labels::yoffset(const length &l)
    {
	y_offset_ = l;
	return *this;
    }
    labels &labels::offset(const length &xoff, const length &yoff)
    {
	x_offset_ = xoff;
	y_offset_ = yoff;
	return *this;
    }

    void labels::prepare_for_draw(plottable *g, frame *, int count)
    {
	if(count==1)
	{
	    x_offset_.register_me();
	    y_offset_.register_me();
	}
    }

    void labels::set_ranges(plottable *g,axis *xaxis,axis *yaxis)
	TRY
    {
	set_ranges_default(g,xaxis,yaxis,1,2);
    }
    CATCH("labels::set_ranges(...)")

    graph_drawer *labels::clone() const
    {
	return new labels(*this);
    }

    void labels::draw(plottable *g,frame *f,terminal *term)
	TRY
    {
	if(!g) err("Plottable to draw not set");
	if(g->empty()) return;

	function getx = g->x_hint();
	function gety = g->y_hint();
	function getz = g->z_hint();
	if(!getx.initialized()) getx = _1;
	if(!gety.initialized()) gety = _2;
	if(!getz.initialized()) getz = _3;

	term->set_color(g->linecolor());

	axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());


	double xmin = xaxis->min();
	if(g->xmin()!=unset && g->xmin()>xmin) xmin = g->xmin();
	double xmax = xaxis->max();
	if(g->xmax()!=unset && g->xmax()<xmax) xmax = g->xmax();

	double ymin = yaxis->min();
	if(g->ymin()!=unset && g->ymin()>ymin) ymin = g->ymin();
	double ymax = yaxis->max();
	if(g->ymax()!=unset && g->ymax()<ymax) ymax = g->ymax();

	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    const var xorig_var = getx.eval(*(g->get(i)));
	    const var yorig_var = gety.eval(*(g->get(i)));
	    if(ignore::it(xorig_var) || ignore::it(yorig_var)) continue;
	    const std::string  zorig = getz.eval(*(g->get(i))).str();
	    const double xorig = xorig_var.dbl();
	    const double yorig = yorig_var.dbl();
	    if ( xaxis->logscale() && !(xorig>0.0) ) continue;
	    if ( yaxis->logscale() && !(yorig>0.0) ) continue;
	    if(xorig < xmin || xmax < xorig || yorig < ymin || ymax < yorig) continue;

	    double xm = xaxis->map_point(xorig);
	    double ym = yaxis->map_point(yorig);

	    term->translate(terminal::id(xm,1),terminal::id(ym,2));
	    term->draw_text(terminal::coord(x_offset_.termspecific_id(),
					    y_offset_.termspecific_id()),
			    zorig,
			    xalign_,yalign_,angle_);
	    term->reset_transformation();
	}
    }
    CATCH("labels::draw(terminal *)")

    labels::labels(sym::position xal, sym::position yal,
		   double ang, const length &xoff, const length &yoff)
    {
	x_offset_ = xoff;
	y_offset_ = yoff;
	xalign_ = xal;
	yalign_ = yal;
	angle_ = ang;
    }

    labels &labels::operator() (sym::position xal, sym::position yal, double angle,
				const length &xoff, const length &yoff)
    {
	x_offset_ = xoff;
	y_offset_ = yoff;
	xalign_ = xal;
	yalign_ = yal;
	angle_ = angle;
	return *this;
    }

    void labels::draw_sample(const length &x,
			     const length &y,
			     const length &size,
			     const plottable *s,
			     terminal *t)
    {
    }

    // ---------- graphd_errorbar   ---------------------------------------------------

    bool errorbars::default_clip_points_ = true;
    bool errorbars::default_clip_errorbars_ = false;

    errorbars &errorbars::endmarker_size(const length &l) { endmarker_size_ = l; return *this; }

    length &errorbars::default_endmarker_size_()
    {
        static length &s = PS;
        return s;
    }

    int errorbars::req_components() const
    {
	int n = 2;
	if(x_) {if(symmetric_x_) n += 1; else n += 2;}
	if(y_) {if(symmetric_y_) n += 1; else n += 2;}
	return n;
    }

    void errorbars::get_functions(const plottable *g,
				  function &getx, function &getx1, function &getx2,
				  function &gety, function &gety1, function &gety2)
    {
	getx = g->x_hint();
	if(!getx.initialized()) getx = _1;

	gety = g->y_hint();
	if(!gety.initialized()) gety = _2;

	getx1 = g->x1_hint();
	getx2 = g->x2_hint();

	if(x_ && (!getx1.initialized() || !getx2.initialized()))
	{
	    if(symmetric_x_)
	    {
		function getdx = g->dx_hint();

		// if plottable did not specify a dx_hint, or if it's an fgraph, then
		// use the 3rd column. 
		// This is a dirty code. Think it over!!!!
		if(!getdx.initialized() || dynamic_cast<const fgraph *>(g))
		{
		    getx1 = getx - _3;
		    getx2 = getx + _3;
		}
		else
		{
		    getx1 = getx - getdx/2;
		    getx2 = getx + getdx/2;
		}
	    }
	    else
	    {
		getx1 = _3;
		getx2 = _4;
	    }
	}

	gety1 = g->y1_hint();
	gety2 = g->y2_hint();

	if(y_ && (!gety1.initialized() || !gety2.initialized()))
	{
	    if(symmetric_y_)
	    {
		function getdy = g->dy_hint();
		if(!getdy.initialized())
		{
		    if(!x_) getdy = _3*2;
		    else
		    {
			if(symmetric_x_) getdy = _4*2;
			else             getdy = _5*2;
		    }
		}
		gety1 = gety - getdy/2;
		gety2 = gety + getdy/2;
	    }
	    else
	    {
		if(!x_)
		{
		    gety1 = _3;
		    gety2 = _4;
		}
		else
		{
		    gety1 = _4;
		    gety2 = _5;
		}
	    }
	}
	if(!getx1.initialized()) getx1 = getx;
	if(!getx2.initialized()) getx2 = getx;
	if(!gety1.initialized()) gety1 = gety;
	if(!gety2.initialized()) gety2 = gety;
    }

    void errorbars::set_ranges(plottable *g,axis *xaxis,axis *yaxis)
    {
	if(!g) return;

	if(g->empty())
	{
	    warning::print("Empty plottable","errorbars::set_ranges(...)");
	    return;
	}

	function getx, getx1, getx2, gety, gety1, gety2;
	get_functions(g,getx, getx1, getx2, gety, gety1, gety2);

	if(g->xmin()!=unset) xaxis->extend_range(g->xmin());
	if(g->xmax()!=unset) xaxis->extend_range(g->xmax());
	if(g->ymin()!=unset) yaxis->extend_range(g->ymin());
	if(g->ymax()!=unset) yaxis->extend_range(g->ymax());

	if ( g->xmin()!=unset && g->xmax()!=unset && g->ymin()!=unset && g->ymax()!=unset ) return;

	double xmin=unset;
	double xmax=unset;
	double ymin=unset;
	double ymax=unset;

	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    const var xorig_var = getx.eval(*(g->get(i)));
	    const var x1orig_var = getx1.eval(*(g->get(i)));
	    const var x2orig_var = getx2.eval(*(g->get(i)));
	    if(ignore::it(xorig_var) || ignore::it(x1orig_var) || ignore::it(x2orig_var)) continue;
	    const double xorig = xorig_var.dbl();
	    const double x1orig = x1orig_var.dbl();
	    const double x2orig = x2orig_var.dbl();
	    if(xaxis->logscale() && !(xorig>0.0)) continue;
	    if((g->xmin()!=unset && xorig<g->xmin()) || (g->xmax()!=unset && xorig>g->xmax())) continue;
	    if(xmin==unset && xmax==unset) { xmin = xorig; xmax=xorig; }
	    if(xorig < xmin) xmin = xorig;
	    if(xorig > xmax) xmax = xorig;
	    if(x2orig > xmax) xmax = x2orig;
	    if( !(xaxis->logscale() && !(x1orig>0.0)) ) { if(x1orig < xmin) xmin = x1orig; }

	    const var yorig_var = gety.eval(*(g->get(i)));
	    const var y1orig_var = gety1.eval(*(g->get(i)));
	    const var y2orig_var = gety2.eval(*(g->get(i)));
	    if(ignore::it(yorig_var) || ignore::it(y1orig_var) || ignore::it(y2orig_var)) continue;
	    const double yorig = yorig_var.dbl();
	    const double y1orig = y1orig_var.dbl();
	    const double y2orig = y2orig_var.dbl();
	    if(yaxis->logscale() && !(yorig>0.0)) continue;
	    if((g->ymin()!=unset && yorig<g->ymin()) || (g->ymax()!=unset && yorig>g->ymax())) continue;
	    if(ymin==unset && ymax==unset) { ymin = yorig; ymax=yorig; }
	    if(yorig < ymin) ymin = yorig;
	    if(yorig > ymax) ymax = yorig;
	    if(y2orig > ymax) ymax = y2orig;
	    if( !(yaxis->logscale() && !(y1orig>0.0)) ) { if(y1orig < ymin) ymin = y1orig; }

	}
	if(g->xmin()==unset && xmin!=unset) xaxis->extend_range(xmin);
	if(g->xmax()==unset && xmax!=unset) xaxis->extend_range(xmax);
	if(g->ymin()==unset && ymin!=unset) yaxis->extend_range(ymin);
	if(g->ymax()==unset && ymax!=unset) yaxis->extend_range(ymax);
    }

    void errorbars::draw(plottable *g,frame *f,terminal *t)
    {
	axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());

	double xmin = xaxis->min();
	if(g->xmin() != unset) xmin = ::max(xmin, g->xmin());

	double xmax = xaxis->max();
	if(g->xmax() != unset) xmax = ::min(xmax, g->xmax());

	double ymin = yaxis->min();
	if(g->ymin() != unset) ymin = ::max(ymin, g->ymin());

	double ymax = yaxis->max();
	if(g->ymax() != unset) ymax = ::min(ymax, g->ymax());

	t->set_linewidth(g->linewidth().termspecific_id());
	t->set_linestyle(g->linestyle());

	bool call_noclip = false;
	if(clip_errorbars_)
	{
	    if(t->clip(terminal::coord(terminal::id(xaxis->map_point(xmin),1),
				       terminal::id(yaxis->map_point(ymin),2)),
		       terminal::coord(terminal::id(xaxis->map_point(xmax),1),
				       terminal::id(yaxis->map_point(ymax),2))))
	    {
		call_noclip = true;
	    }
	    else
	    {
		warning::print("This terminal does not support clipping, you might see something different from what you expect",
			       "errorbars::draw(...)");
	    }
	}

	function getx, getx1, getx2, gety, gety1, gety2;
	get_functions(g,getx, getx1, getx2, gety, gety1, gety2);

	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    const var xorig_var = getx.eval(*(g->get(i))); 
	    const var yorig_var = gety.eval(*(g->get(i))); 
	    if(ignore::it(xorig_var) || ignore::it(yorig_var)) continue;

	    const double xorig = xorig_var.dbl();
	    const double yorig = yorig_var.dbl();

	    if(clip_points_)
	    {
		if(xorig < xmin || xmax < xorig) continue;
		if(yorig < ymin || ymax < yorig) continue;
	    }

	    if(xaxis->logscale() && xorig <= 0) continue;
	    if(yaxis->logscale() && yorig <= 0) continue;

	    double x_m = xaxis->map_point(xorig);
	    double x1_m=-0.1,x2_m=1.1;
	    if(x_)
	    {
		var x1_var = getx1.eval(*(g->get(i)));
		var x2_var = getx2.eval(*(g->get(i)));
		if(ignore::it(x1_var) || ignore::it(x2_var)) continue;
		double x1=x1_var.dbl();
		double x2=x2_var.dbl();

		if(x1 > xorig)
		{
		    warning::print("Lower edge of x-errorbar larger than the point itself. Correcting",
				   "errorbars::draw(...)");
		    x1 = xorig;
		}
		if(x2 < xorig)
		{
		    warning::print("Upper edge of x-errorbar smaller than the point itself. Correcting",
				   "errorbars::draw(...)");
		    x2 = xorig;
		}

		if(!xaxis->logscale())
		{
		    x1_m = xaxis->map_point(x1);
		    x2_m = xaxis->map_point(x2);
		}
		else
		{
		    if(x1 > 0) x1_m = xaxis->map_point(x1);
		    else x1_m = 0;
		    if(x2 > 0) x2_m = xaxis->map_point(x2);
		    else x2_m = 0;
		}
	    }

	    double y_m = yaxis->map_point(yorig);
	    double y1_m=-0.1,y2_m=1.1;
	    if(y_)
	    {
		var y1_var = gety1.eval(*(g->get(i)));
		var y2_var = gety2.eval(*(g->get(i)));
		if(ignore::it(y1_var) || ignore::it(y2_var)) continue;
		double y1=y1_var.dbl();
		double y2=y2_var.dbl();

		if(y1 > yorig)
		{
		    warning::print("Lower edge of y-errorbar larger than the point itself. Correcting",
				   "errorbars::draw(...)");
		    y1 = yorig;
		}
		if(y2 < yorig)
		{
		    warning::print("Upper edge of y-errorbar smaller than the point itself. Correcting",
				   "errorbars::draw(...)");
		    y2 = yorig;
		}

		if(!yaxis->logscale())
		{
		    y1_m = yaxis->map_point(y1);
		    y2_m = yaxis->map_point(y2);
		}
		else
		{
		    if(y1 > 0) y1_m = yaxis->map_point(y1);
		    else y1_m = 0;
		    if(y2 > 0) y2_m = yaxis->map_point(y2);
		    else y2_m = 0;
		}
	    }

	    t->set_color(g->linecolor());
            int i1 = t->newlength();
            int i2 = t->newlength();
	    if(x_)
	    {
		t->draw_line(terminal::coord(terminal::id(x1_m,1),terminal::id(y_m,2)),
			     terminal::coord(terminal::id(x2_m,1),terminal::id(y_m,2)));
                t->overwrite(i1,1,terminal::id(y_m,2),-0.5,endmarker_size_.termspecific_id());
                t->overwrite(i2,1,terminal::id(y_m,2),+0.5,endmarker_size_.termspecific_id());
                t->draw_line(terminal::coord(terminal::id(x1_m,1),terminal::id(i1)),
                             terminal::coord(terminal::id(x1_m,1),terminal::id(i2)));
                t->draw_line(terminal::coord(terminal::id(x2_m,1),terminal::id(i1)),
                             terminal::coord(terminal::id(x2_m,1),terminal::id(i2)));
	    }
	    if(y_)
	    {
		t->draw_line(terminal::coord(terminal::id(x_m,1),terminal::id(y1_m,2)),
			     terminal::coord(terminal::id(x_m,1),terminal::id(y2_m,2)));
                t->overwrite(i1,1,terminal::id(x_m,1),-0.5,endmarker_size_.termspecific_id());
                t->overwrite(i2,1,terminal::id(x_m,1),+0.5,endmarker_size_.termspecific_id());
                t->draw_line(terminal::coord(terminal::id(i1),terminal::id(y1_m,2)),
                             terminal::coord(terminal::id(i2),terminal::id(y1_m,2)));
                t->draw_line(terminal::coord(terminal::id(i1),terminal::id(y2_m,2)),
                             terminal::coord(terminal::id(i2),terminal::id(y2_m,2)));
	    }
		
	    if(g->pointtype())
	    {
		t->set_color(g->pointcolor());
		t->translate(terminal::id(x_m,1), terminal::id(y_m,2));
		g->pointtype()->draw(t);
		t->reset_transformation();
	    }
	}
	if(call_noclip) t->noclip();
    }

    void errorbars::draw_sample(const length &x,const length &y,const length &size,
				const plottable *style, terminal *t)
    {
	t->set_color(style->linecolor());
	t->set_linewidth(style->linewidth().termspecific_id());
	t->set_linestyle(style->linestyle());

	length l1 = !x - 0.5*!size;
	length l2 = !x + 0.5*!size;

	l1.specialize(t);
	l2.specialize(t);

	t->draw_line(terminal::coord(l1.termspecific_id(), y.termspecific_id()),
		     terminal::coord(l2.termspecific_id(), y.termspecific_id()));

	if(style->pointtype())
	{
	    t->set_color(style->pointcolor());
	    t->translate(x.termspecific_id(), y.termspecific_id());
	    style->pointtype()->draw(t);
	    t->reset_transformation();
	}
    }

    errorbars &errorbars::clip_points(bool f)
    {
	clip_points_ = f;
	return *this;
    }

    errorbars &errorbars::clip_errorbars(bool f)
    {
	clip_errorbars_ = f;
	return *this;
    }

    

    graph_drawer *errorbars::clone() const
    {
	return new errorbars(*this);
    }

    errorbars::errorbars(bool x,bool sx,bool y,bool sy)
    {
	x_ = x;
	symmetric_x_ = sx;
	y_ = y;
	symmetric_y_ = sy;
	clip_errorbars_ = default_clip_errorbars_;
	clip_points_    = default_clip_points_;
        endmarker_size_ = default_endmarker_size_();
    }

    errorbars::errorbars(const errorbars &rhs)
    {
	x_ = rhs.x_;
	symmetric_x_ = rhs.symmetric_x_;
	y_ = rhs.y_;
	symmetric_y_ = rhs.symmetric_y_;
	clip_errorbars_ = rhs.clip_errorbars_;
	clip_points_    = rhs.clip_points_;
        endmarker_size_ = rhs.endmarker_size_;
    }

    void errorbars::prepare_for_draw(plottable *g, frame *, int count)
    {
	if(count==1)
	{
	    if(g->pointtype()) g->pointtype()->prepare_for_draw(g->pointsize());
	    g->linewidth().register_me();
            endmarker_size_.register_me();
	}
    }
    

    // ---------- graphd_ticlabels  ---------------------------------------------------


    ticlabels::ticlabels(bool x, bool y,const graph_drawer &d) : x_(x), y_(y)
    {
	drawer_ = d.clone();
    }

    ticlabels::ticlabels(const ticlabels &rhs)
    {
	x_ = rhs.x_;
	y_ = rhs.y_;
	drawer_ = rhs.drawer_->clone();
    }

    ticlabels ticlabels::operator() (const graph_drawer &d)
    {
	if(drawer_) cerr<<"drawer_ should be deleted in ticlabels::operator()"<<endl;
	drawer_ = d.clone();
	return *this;
    }

    void ticlabels::set_ranges(plottable *g, axis *xaxis, axis *yaxis)
    {
	if(drawer_ == 0)
	{
	    cerr<<"This should not happen: drawer_ = 0 in ticlabels::set_ranges(...)"<<endl;
	    exit(1);
	}

	drawer_->set_ranges(g,xaxis,yaxis);

	// set the tic labels first

	vector<double> xtic_pos, ytic_pos;
	vector<var>    xtic_lab, ytic_lab;

	function getx = drawer_->get_x(g);
	function gety = drawer_->get_y(g);
	function getxl = 0;
	function getyl = 0;
	if(x_)
	{
	    if(g->columns() < drawer_->req_components()+1)
	    {
		warning::print("Graph has too few columns, using x values as xlabels",
			       "ticlabels::set_scale");
		getxl = getx;
	    }
	    else
	    {
		getxl = ARG(drawer_->req_components()+1);
	    }

	    if(y_)
	    {
		if(g->columns() < drawer_->req_components()+2)
		{
		    warning::print("Graph has too few columns, using y values as ylabels",
				   "ticlabels::set_scale");
		    getyl = gety;
		}
		else
		{
		    getyl = ARG(drawer_->req_components()+2);
		}
	    }
	}
	else
	{
	    if(g->columns() < drawer_->req_components()+1)
	    {
		warning::print("Graph has too few colums, using y values as ylabels",
			       "ticlabels::set_scale");
		getyl = gety;
	    }
	    else
	    {
		getyl = ARG(drawer_->req_components()+1);
	    }
	}

	double xmin = xaxis->min();
	if(g->xmin()!=unset && g->xmin()>xmin) xmin = g->xmin();
	double xmax = xaxis->max();
	if(g->xmax()!=unset && g->xmax()<xmax) xmax = g->xmax();

	double ymin = yaxis->min();
	if(g->ymin()!=unset && g->ymin()>ymin) ymin = g->ymin();
	double ymax = yaxis->max();
	if(g->ymax()!=unset && g->ymax()<ymax) ymax = g->ymax();

	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    const var x_var = getx.eval(*(g->get(i)));
	    if ( ignore::it(x_var) ) continue;
	    const double x = x_var.dbl(); 
	    if ( xaxis->logscale() && !(x>0.0) ) continue;
	    if ( x<xmin || x>xmax ) continue;
	    const var y_var = gety.eval(*(g->get(i)));
	    if ( ignore::it(y_var) ) continue;
	    const double y = y_var.dbl(); 
	    if ( yaxis->logscale() && !(y>0.0) ) continue;
	    if ( y<ymin || y>ymax ) continue;

	    if(x_)
	    {
		var label = getxl.eval(*(g->get(i)));
		bool already_put = false;
		for(unsigned int i=0; i<xtic_pos.size(); ++i)
		{
		    if(xtic_pos[i] == x && xtic_lab[i].str() == label.str())
		    {
			already_put = true;
			break;
		    }
		}
		if(!already_put)
		{
		    xaxis->tics(x,label);
		    xtic_pos.push_back(x);
		    xtic_lab.push_back(label);
		}
	    }

	    if(y_)
	    {
		const double y = gety.eval(*(g->get(i))).dbl();
		var label      = getyl.eval(*(g->get(i)));
		bool already_put = false;
		for(unsigned int i=0; i<ytic_pos.size(); ++i)
		{
		    if(ytic_pos[i] == y && ytic_lab[i].str() == label.str())
		    {
			already_put = true;
			break;
		    }
		}
		if(!already_put)
		{
		    yaxis->tics(y,label);
		    ytic_pos.push_back(y);
		    ytic_lab.push_back(label);
		}
	    }
	}
    }

    void ticlabels::draw(plottable *g, frame *f, terminal *t)
    {
	drawer_->draw(g,f,t);
    }

    void ticlabels::draw_sample(const length &x, const length &y,
				const length &size, const plottable *style, terminal *t)
    {
	drawer_->draw_sample(x,y,size,style,t);
    }

    graph_drawer *ticlabels::clone() const
    {
	return new ticlabels(*this);
    }

    void ticlabels::prepare_for_draw(plottable *g, frame *f, int count)
    {
	if(count==1) drawer_->prepare_for_draw(g,f,count);
    }

    int ticlabels::req_components() const
    {
	// Since the ticlabels drawstyle uses the x/y values,
	// if no extra colums are provided, therefore it return
	// the same number of required columns as the wrapped drawstyle
	return drawer_->req_components();
    }


    // ---------- graphd_colorscale ---------------------------------------------------


    graphd_colorscale &graphd_colorscale::draw_colorlegend(bool f)
    {
	if(f) draw_colorlegend_ = 2;
	else draw_colorlegend_ = 1;
	return *this;
    }

    graphd_colorscale &graphd_colorscale::color_min(double v)
    {
	color_min_ = v;
	if(color_min_ == unset) color_min_fixed_ = false;
	else color_min_fixed_ = true;
	return *this;
    }
    graphd_colorscale &graphd_colorscale::color_max(double v)
    {
	color_max_ = v;
	if(color_max_ == unset) color_max_fixed_ = false;
	else color_max_fixed_ = true;
	return *this; 
    }

    graphd_colorscale &graphd_colorscale::color_range(double min, double max)
    {
	color_min(min);
	color_max(max);
	return *this;
    }

    graphd_colorscale &graphd_colorscale::color_range(const color &startcolor, const color &endcolor)
    {
	cmapping_.set(startcolor,endcolor);
	return *this; 
    }

    void graphd_colorscale::set_colorlegend(color_legend *l)
    {
	// If it has already a colorlegend, unregister myself
	// from that colorlegend's owners, and delete it
	if(colorlegend_)
	{
	    colorlegend_->remove_owner(this);
	    if(colorlegend_->owners_.empty() && colorlegend_->autodel())
		delete colorlegend_;
	    colorlegend_ = 0;
	}

	// and set now the new link
	colorlegend_ = l;
	if(colorlegend_) colorlegend_->owners_.push_back(this);
    }

    void graphd_colorscale::setup_when_added(plottable *p, frame *f)
    {
	// If the colorlegend has already been specified, we have
	// nothing to do here...
	if(colorlegend_ != 0) return;

	if(!f || !p) return;
	
	if(draw_colorlegend_ == 0)  // do not create at all (i.e. delete it)
	{
	    if(colorlegend_)
	    {
		colorlegend_->remove_owner(this);
		if(colorlegend_->owners_.empty() && colorlegend_->autodel())
		    delete colorlegend_;
		colorlegend_ = 0;
	    }
	    return;
	}

	// If the user specified a legendname (i.e. a name identifying a legend
	// to be shared), first of all search for it:
	if(legendname_ != "")
	{
	    // Search for the topmost container, i.e. the canvas
	    container *topcanvas = f;
	    for(; topcanvas->parent(); topcanvas=topcanvas->parent());

	    // If a colorlegend with this name is already present in the canvas,
	    // use that one, and return.
	    if(color_legend *l = dynamic_cast<color_legend*>(topcanvas->find("colorlegend_" + legendname_)))
	    {
		colorlegend_ = l;
		colorlegend_->owners_.push_back(this);
		return;
	    }
	}

	// loop over all the graphs in the given frame
	for(int i=0; i<f->ngraphs(); ++i)
	{
	    plottable *g = f->get_graph(i);

	    // Check if the drawstyle of this graph is a graphd_colorscale type
	    // if so, share the color_legend with that graph
	    if(graphd_colorscale *c = dynamic_cast<graphd_colorscale*>(g->drawstyle()))
	    {
		if(c->colorlegend_)
		{
		    colorlegend_ = c->colorlegend_;
		    colorlegend_->owners_.push_back(this);
		    break;
		}
	    }
	}

	// If no such graph was found, i.e. colorlegend is still 0,
	// create it.
	if(!colorlegend_)
	{
            // If we are within an mframe, check all graphs in the same row of the mframe
            // whether they have a colorlegend. If they do, just share it.
            if(mframe *mf = dynamic_cast<mframe*>(f->parent()))
            {
                int row = mf->row(f);
                if(row==0) warning::print("This should never happen","graphd_colorscale::setup_when_added(...)");
                else
                {
                    for(int col=1; col<=mf->ncols(); ++col)
                    {
                        for(int i=0; i<(*mf)(col,row)->ngraphs(); ++i)
                        {
                            plottable *g = (*mf)(col,row)->get_graph(i);
                            if(graphd_colorscale *c = dynamic_cast<graphd_colorscale*>(g->drawstyle()))
                            {
                                if(c->colorlegend_)
                                {
                                    colorlegend_ = c->colorlegend_;
                                    colorlegend_->owners_.push_back(this);
                                    break;
                                }
                            }
                        }
                    }
                }
            }
        }
        if(!colorlegend_)
        {
            colorlegend_ = new color_legend(this);
            colorlegend_->autodel(true);
            f->rmarginobject(colorlegend_);
            
            // if the legendname (for sharing) is specified, but it did not exist before,
            // then we need to set its name now.
            colorlegend_->name("colorlegend_" + legendname_);
	}
    }

    graphd_colorscale::~graphd_colorscale()
    {
	if(colorlegend_)
	{
	    // Remove myself (this) from the owners of the colorlegend first
	    colorlegend_->remove_owner(this);
	    
	    // If there are no more owners of the colorlegend, delete it.
	    if(colorlegend_->owners_.empty()) delete colorlegend_;
	}
    }

    graphd_colorscale::graphd_colorscale()
    {
	legendname_ = "";
	colorlegend_ = 0;
	color_min_ = color_max_ = unset;
	color_min_fixed_ = color_max_fixed_ = false;
	color_logscale_ = false;
	color_samples_ = 80;
	underflow_color_ = transparent;
	overflow_color_ = transparent;
	draw_colorlegend_ = 2;
    }

    graphd_colorscale::graphd_colorscale(const graphd_colorscale &rhs)
    {
	copy(rhs);
    }
    
    void graphd_colorscale::copy(const graphd_colorscale &rhs)
    {
	color_min_ = rhs.color_min_;
	color_max_ = rhs.color_max_;
	color_logscale_ = rhs.color_logscale_;
	color_min_fixed_ = rhs.color_min_fixed_;
	color_max_fixed_ = rhs.color_max_fixed_;
	color_samples_ = rhs.color_samples_;
	cmapping_ = rhs.cmapping_;
	underflow_color_ = rhs.underflow_color_;
	overflow_color_  = rhs.overflow_color_;
	color_title_ = rhs.color_title_;
	draw_colorlegend_ = rhs.draw_colorlegend_;
	colorlegend_ = rhs.colorlegend_;
	if(colorlegend_) colorlegend_->owners_.push_back(this);
	legendname_ = rhs.legendname_;
    }

    color graphd_colorscale::map_color(double val, double mini, double maxi) const
    {
	if(color_logscale_ && (val<=0 || mini<=0 || maxi<=0))
	{
	    warning::print("Range must be positive for logscale",
			   "graphd_colorscale::map_color(...)");
	    return color(0,0,0);
	}

	if(val < mini && underflow_color_ != transparent) return underflow_color_;
	if(val > maxi && overflow_color_  != transparent) return overflow_color_;

	val  = (color_logscale_?::log10(val):val);
	mini = (color_logscale_?::log10(mini):mini);
	maxi = (color_logscale_?::log10(maxi):maxi);
	if(val<mini) val=mini;
	if(val>maxi) val=maxi;
	return cmapping_.map(val,mini,maxi);
    }



    // ---------- cbox_base  ----------------------------------------------------------

    void cbox_base::init_lengths()
    {
    }

    cbox_base::~cbox_base()
    {
    }

    // for some reason this construction does not work, the default constructor of
    // colorscale is called, I don't know why...
    //cbox_base::cbox_base(const cbox_base &rhs) : colorscale((colorscale)rhs)
    cbox_base::cbox_base(const cbox_base &rhs)
    {
	graphd_colorscale::copy(rhs);
    }

    cbox_base::cbox_base(void *cmapfunc)
    {
	if(cmapfunc != 0) map_function(cmapfunc);
	init_lengths();
    }

    cbox_base::cbox_base(color (*cmapfunc)(double,double,double))
    {
	if(cmapfunc != 0) map_function(cmapfunc);
	init_lengths();
    }

    cbox_base::cbox_base(double mini, double maxi)
    {
	color_range(mini,maxi);
	init_lengths();
    }

    cbox_base::cbox_base(const color &start, const color &end)
    {
	color_range(start,end);
	init_lengths();
    }

    cbox_base::cbox_base()
    {
	color_range(unset,unset);
	init_lengths();
    }
    
    void cbox_base::draw_sample(const length &x,
				const length &y,
				const length &size,
				const plottable *s,
				terminal *t)
    {
	vector<terminal::coord> c;
	length l1 = !x - 0.5*!size;
	length l2 = !x + 0.5*!size;
	l1.specialize(t);
	l2.specialize(t);
	int nsample=(color_samples()>10 ? 10 : color_samples());
	for ( int i=0 ; i<nsample ; ++i )
	{
	    double val = 0.0 + (1.0-0.0)/nsample*(i+1);
	    color the_color = map_color(val,0.0,1.0);
	    length l3 = !y - 0.5*0.8*EM + 0.8*(double)i/nsample*EM;
	    length l4 = !y - 0.5*0.8*EM + 0.8*(double)(i+1)/nsample*EM;
	    l3.specialize(t);
	    l4.specialize(t);
	    c.clear();
	    c.push_back(terminal::coord(l1.termspecific_id(),l3.termspecific_id()));
	    c.push_back(terminal::coord(l2.termspecific_id(),l3.termspecific_id()));
	    c.push_back(terminal::coord(l2.termspecific_id(),l4.termspecific_id()));
	    c.push_back(terminal::coord(l1.termspecific_id(),l4.termspecific_id()));
	    c.push_back(terminal::coord(l1.termspecific_id(),l3.termspecific_id()));
	    t->set_color(the_color);
	    t->fill_polygon(c);
	}
    }

    // ---------- graphd_cbox  --------------------------------------------------------

/*
  bool cboxes::default_frame_foreground_ = true;
  bool cboxes::default_grid_foreground_ = false;

  cboxes &cboxes::legend(color_legend *l)
  {
  draw_clegend_ = 2;
  set_clegend(l);
  return *this; 
  }

  cboxes &cboxes::legend(const var &legendname)
  {
  legendname_ = legendname.str();
  return *this;
  }


  void cboxes::prepare_for_draw(plottable *g, frame *f, int count)
  {
  if(global::debug>0) cout<<"[blop] [cboxes] prepare_for_draw starts... pass="<<count<<endl;
  if(count==1)
  {
  if(frame_foreground_) f->foreground(true);
  if(grid_foreground_) f->grid_foreground(true);
	    
  axis *xaxis =
  (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
  axis *yaxis =
  (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());
	    
  calc(g,xaxis,yaxis);
  if(colorlegend_)
  {
  colorlegend_->clear_tics();
  if(draw_clegend_<2) colorlegend_->off(); // Switch it off
  if(ctitle_.str()!="")
  {
  var new_title = colorlegend_->title();
  if(new_title.str() != "") new_title &= ", ";
  new_title &= ctitle_;
  colorlegend_->title(new_title);
  }
  }
  }

  if(count==2)
  {
  // Let the colorlegend calculate the unified range for all
  // of its graphs
  if(colorlegend_)
  {
  colorlegend_->calculate_tics();
  if(draw_clegend_>=2) colorlegend_->on();  // if any of the graphs have legend on...
  }
  }
  if(global::debug>0) cout<<"[blop] [boxes] prepare_for_draw finished."<<endl;
  }

  cboxes::cboxes(const cboxes &rhs) : cbox_base(rhs)
  {
  dx_=rhs.dx_;
  dy_=rhs.dy_;
  frame_foreground_ = rhs.frame_foreground_;
  grid_foreground_  = rhs.grid_foreground_;
  skip_outrange_ = rhs.skip_outrange_;
  }

  cboxes::cboxes(void *p) : cbox_base(p) {}
  cboxes::cboxes(color (*p)(double,double,double)) : cbox_base(p) {}

  cboxes::cboxes(double mini, double maxi) : cbox_base(mini, maxi)
  {
  dx_ = dy_ = unset;
  frame_foreground_ = default_frame_foreground_;
  grid_foreground_  = default_grid_foreground_;
  skip_outrange_ = true;
  }

  cboxes::cboxes(const color &start, const color &end) : cbox_base(start,end)
  {
  dx_ = dy_ = unset;
  frame_foreground_ = default_frame_foreground_;
  grid_foreground_  = default_grid_foreground_;
  skip_outrange_ = true;
	
  }

  cboxes::cboxes() : cbox_base()
  {
  dx_ = dy_ = unset;
  frame_foreground_ = default_frame_foreground_;
  grid_foreground_  = default_grid_foreground_;
  skip_outrange_ = true;
  }

  void cboxes::set_ranges(plottable *g,axis *x,axis *y)
  {
  // prepare_for_draw is called before this function, so everything
  // is set up already

  x->autoextend_min_soft(false);
  x->autoextend_max_soft(false);
  y->autoextend_min_soft(false);
  y->autoextend_max_soft(false);
  }


  static void get_cellsize(plottable *g, const function &getx, const function &gety, const function &getdx, const function &getdy,
  double &dx, double &dy)
  {
  dx = dy = unset;

  if(!getx.initialized() || !gety.initialized())
  {
  warning::print("getx or gety is uninitialized","get_cellsize(...)");
  return;
  }

  if ( g->ordered() >= 2 )
  {
  double xlast=unset, ylast=unset;
  for(plottable::size_type i1=0; i1<g->size(); ++i1)
  {
  const var x_var = getx.eval(*(g->get(i1)));
  const var y_var = gety.eval(*(g->get(i1)));
  if(ignore::it(x_var) || ignore::it(y_var)) continue;
  double x = x_var.dbl();
  double y = y_var.dbl();
  if(g->xmin() != unset && x < g->xmin()) continue;
  if(g->xmax() != unset && x > g->xmax()) continue;
  if(g->ymin() != unset && y < g->ymin()) continue;
  if(g->ymax() != unset && y > g->ymax()) continue;
		    
  if(xlast != unset && x != xlast)
  {
  if(dx == unset || ::fabs(x-xlast)<dx) dx = ::fabs(x-xlast);
  }
  if(ylast != unset && y != ylast)
  {
  if(dy == unset || ::fabs(y-ylast)<dy) dy = ::fabs(y-ylast);
  }
  xlast = x;
  ylast = y;
  }
  }
	    
  else if( getdx.initialized() && getdy.initialized() )
  {
  for(plottable::size_type i=0; i<g->size(); ++i)
  {
  const datapoint *pp = g->get(i);
  const var x_var=getx.eval(*pp);
  const var y_var=gety.eval(*pp);
  const var dx_tmp_var=getdx.eval(*pp);
  const var dy_tmp_var=getdy.eval(*pp);
  if ( ignore::it(x_var) || ignore::it(y_var) ||
  ignore::it(dx_tmp_var) || ignore::it(dy_tmp_var) ) continue;
  const double dx_tmp = dx_tmp_var.dbl();
  const double dy_tmp = dy_tmp_var.dbl();
  if(dx == unset || dx_tmp < dx) dx = dx_tmp;
  if(dy == unset || dy_tmp < dy) dy = dy_tmp;
  }
  }

  // otherwise, if not ordered, and no dx/dy hints are given, then
  // loop over all data point pairs (very slow for large data!!!)
  else 
  {
  map<double,bool> xvalues,yvalues;
  for(plottable::size_type i=0; i<g->size(); ++i)
  {
  const var xorig_var = getx.eval(*(g->get(i)));
  const var yorig_var = gety.eval(*(g->get(i)));
  if(ignore::it(xorig_var) || ignore::it(yorig_var)) continue;
  const double xorig = xorig_var.dbl();
  const double yorig = yorig_var.dbl();
  xvalues[xorig] = true;
  yvalues[yorig] = true;
  }

  for(map<double,bool>::iterator i=xvalues.begin(); i!=xvalues.end(); ++i)
  {
  if(i != xvalues.begin())
  {
  map<double,bool>::iterator prev = i;
  --prev;
  if(dx == unset || ::fabs((*prev).first - (*i).first) < dx)
  dx = ::fabs((*prev).first - (*i).first);
  }
  }
  for(map<double,bool>::iterator i=yvalues.begin(); i!=yvalues.end(); ++i)
  {
  if(i != yvalues.begin())
  {
  map<double,bool>::iterator prev = i;
  --prev;
  if(dy == unset || ::fabs((*prev).first - (*i).first) < dy)
  dy = ::fabs((*prev).first - (*i).first);
  }
  }
  }
  }

  void cboxes::calc(plottable *g,axis *xaxis,axis *yaxis)
  TRY
  {

  if(global::debug>0) cout<<"[blop] [DEBUG] cboxes::calc starts..."<<endl;
  if(!g) err("Plottable to draw not set");

  function getx = get_x(g);
  function gety = get_y(g);
  function getdx = g->dx_hint();
  function getdy = g->dy_hint();

  function getz = g->z_hint();
  if(!getz.initialized()) getz = _3;

  double dx_new=unset, dy_new = unset;

  // only calculate the cellsize if any of dx_ or dy_ are unset. othersize it's not necessary
  if(dx_ == unset || dy_ == unset) get_cellsize(g,getx,gety,getdx,getdy,dx_new,dy_new);

  if ( dx_==unset ) dx_=dx_new;
  if ( dy_==unset ) dy_=dy_new;
	

  if(dx_ == unset || dx_ <=0)
  {
  warning::print("Couldn't figure out boxwidth, using dx=1.0",
  "cboxes::calc(...)");
  dx_ = 1.0;
  }
  if(dy_ == unset || dy_ <= 0)
  {
  warning::print("Couldn't figure out boxheight, using dy=1.0",
  "cboxes::calc(...)");
  dy_ = 1.0;
  }
	
  // Set ranges:
  if(g->xmin()!=unset) xaxis->extend_range(g->xmin());
  if(g->xmax()!=unset) xaxis->extend_range(g->xmax());
  if(g->ymin()!=unset) yaxis->extend_range(g->ymin());
  if(g->ymax()!=unset) yaxis->extend_range(g->ymax());

  if ( g->xmin()==unset ||
  g->xmax()==unset ||
  g->ymin()==unset ||
  g->ymax()==unset ||
  color_min_fixed_==false ||
  color_max_fixed_==false )
  {
  double xmin=unset;
  double xmax=unset;
  double ymin=unset;
  double ymax=unset;
  double vmin=unset;
  double vmax=unset;

  double xmin_accept = g->xmin();
  if(xaxis->min_fixed() && (xmin_accept==unset || xaxis->min()<xmin_accept))
  xmin_accept = xaxis->min();
  double xmax_accept = g->xmax();
  if(xaxis->max_fixed() && (xmax_accept==unset || xaxis->max()<xmax_accept))
  xmax_accept = xaxis->max();
  double ymin_accept = g->ymin();
  if(yaxis->min_fixed() && (ymin_accept==unset || yaxis->min()<ymin_accept))
  ymin_accept = yaxis->min();
  double ymax_accept = g->ymax();
  if(yaxis->max_fixed() && (ymax_accept==unset || yaxis->max()<ymax_accept))
  ymax_accept = yaxis->max();

  if(xmin_accept!=unset) xmin = xmin_accept;
  if(xmax_accept!=unset) xmax = xmax_accept;
  if(ymin_accept!=unset) ymin = ymin_accept;
  if(ymax_accept!=unset) ymax = ymax_accept;

  for(plottable::size_type i=0; i < g->size(); ++i)
  {
  const var xorig_var = getx.eval(*(g->get(i)));
  const var yorig_var = gety.eval(*(g->get(i)));
  if ( ignore::it(xorig_var) || ignore::it(yorig_var) ) continue;
  const double xorig=xorig_var.dbl();
  const double yorig=yorig_var.dbl();
  if ( (xaxis->logscale() && !(xorig-0.5*dx_>0.0)) ||
  (yaxis->logscale() && !(yorig-0.5*dy_>0.0)) ) continue;

  if ( (xmin_accept!=unset && xorig+0.5*dx_<xmin_accept) ||
  (xmax_accept!=unset && xorig-0.5*dx_>xmax_accept) ||
  (ymin_accept!=unset && yorig+0.5*dy_<ymin_accept) ||
  (ymax_accept!=unset && yorig-0.5*dy_>ymax_accept) ) continue;

  if(xmin==unset || xorig-0.5*dx_<xmin) xmin=xorig-0.5*dx_;
  if(xmax==unset || xorig+0.5*dx_>xmax) xmax=xorig+0.5*dx_;
  if(ymin==unset || yorig-0.5*dy_<ymin) ymin=yorig-0.5*dy_;
  if(ymax==unset || yorig+0.5*dy_>ymax) ymax=yorig+0.5*dy_;
		
  const var v_var = getz.eval(*(g->get(i)));
  if ( ignore::it(v_var) ) continue;
  const double v = v_var.dbl();
  if(clogscale_ && !(v>0.0)) continue;
  if((color_min_fixed_ && v<color_min_) || (color_max_fixed_ && v>color_max_)) continue;
  if(vmin==unset && vmax==unset) { vmin=v; vmax=v; }
  if( v > vmax) vmax = v;
  else if( v < vmin) vmin = v;
  }
  if(xmin!=unset) xaxis->extend_range(xmin);
  if(xmax!=unset) xaxis->extend_range(xmax);
  if(ymin!=unset) yaxis->extend_range(ymin);
  if(ymax!=unset) yaxis->extend_range(ymax);
  if(color_min_fixed_==false && vmin!=unset) color_min_=vmin;
  if(color_max_fixed_==false && vmax!=unset) color_max_=vmax;

  }
  if(global::debug>0) cout<<"[blop] [DEBUG] cboxes::calc finished."<<endl;
  }
  CATCH("cboxes::calc(...)")

  void cboxes::draw(plottable *g,frame *f,terminal *t)
  TRY
  {
  if(!g || g->empty())  return;
	
  t->reset_transformation();
	
  axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
  axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());
	
  double xmin = xaxis->min();
  if(g->xmin()!=unset && g->xmin()>xmin) xmin = g->xmin();
  double xmax = xaxis->max();
  if(g->xmax()!=unset && g->xmax()<xmax) xmax = g->xmax();

  double ymin = yaxis->min();
  if(g->ymin()!=unset && g->ymin()>ymin) ymin = g->ymin();
  double ymax = yaxis->max();
  if(g->ymax()!=unset && g->ymax()<ymax) ymax = g->ymax();
	
  vector<terminal::coord> c;
	
  function getx = get_x(g);
  function gety = get_y(g);
  function getz = g->z_hint();
  if(!getz.initialized()) getz = _3;
	
  for(plottable::size_type i=0; i<g->size(); ++i)
  {
  const var xorig_var = getx.eval(*(g->get(i)));
  const var yorig_var = gety.eval(*(g->get(i)));
  if(ignore::it(xorig_var) || ignore::it(yorig_var)) continue;
  const double xorig = xorig_var.dbl();
  const double yorig = yorig_var.dbl();
	    
  if(xaxis->logscale() && xorig-0.5*dx_<=0.0) continue;
  if(yaxis->logscale() && yorig-0.5*dy_<=0.0) continue;
  if(xorig+0.5*dx_<xmin || xorig-0.5*dx_>xmax ||
  yorig+0.5*dy_<ymin || yorig-0.5*dy_>ymax) continue;
	    
  double x1 = xaxis->map_point(xorig - dx_/2);
  double x2 = xaxis->map_point(xorig + dx_/2);

  double y1 = yaxis->map_point(yorig - dy_/2);
  double y2 = yaxis->map_point(yorig + dy_/2);

  // make the boxes 'slightly' bigger - otherwise bitmap outputs (jpg, png) and
  // antialiased postscript/pdf will appear with white lines
  x1 -= ::fabs(x2-x1)/80;
  x2 += ::fabs(x2-x1)/80;
  y1 -= ::fabs(y2-y1)/80;
  y2 += ::fabs(y2-y1)/80;
	    
  if(0<=x1 && x1<=1 &&
  0<=x2 && x2<=1 &&
  0<=y1 && y1<=1 &&
  0<=y2 && y2<=1)
  {
  const var zorig_var = getz.eval(*(g->get(i)));
  if ( ignore::it(zorig_var) ) continue;
  double zorig = zorig_var.dbl();
  if (clogscale_ && zorig<=0.0) continue;

  if(zorig<color_min_ && underflow_color_ == transparent && skip_outrange_)
  {
  continue;
  }
  if(zorig>color_max_ && overflow_color_ == transparent && skip_outrange_)
  {
  continue;
  }

  double minimum = min();
  double maximum = max();
  if(colorlegend_)
  {
  minimum = colorlegend_->min();
  maximum = colorlegend_->max();
  }
  const color the_color = map_color(zorig,minimum,maximum);
  c.clear();
  c.push_back(terminal::coord(terminal::id(x1,1),terminal::id(y1,2)));
  c.push_back(terminal::coord(terminal::id(x1,1),terminal::id(y2,2)));
  c.push_back(terminal::coord(terminal::id(x2,1),terminal::id(y2,2)));
  c.push_back(terminal::coord(terminal::id(x2,1),terminal::id(y1,2)));
  c.push_back(terminal::coord(terminal::id(x1,1),terminal::id(y1,2)));
  t->set_color(the_color);
  t->fill_polygon(c);
  }
  }
  }
  CATCH("cboxes::draw(...)")

  graph_drawer *cboxes::clone() const
  {
  cboxes *r =  new cboxes(*this);
  return r;
  }

*/

    // ---------- csboxes   -----------------------------------------------------------

    bool csboxes::default_frame_foreground_ = true;
    bool csboxes::default_grid_foreground_ = false;

    csboxes &csboxes::colormapping(color(*p)(double,double,double))
    {
	cmapping_.set(p);
	return *this;
    }

    csboxes &csboxes::independent_xy(bool f)
    {
	if(f)
	{
	    if(color_func_.is_constant())
	    {
		boxsize_x_func_ = _3;
		boxsize_y_func_ = _4;
	    }
	    else
	    {
		boxsize_x_func_ = _4;
		boxsize_y_func_ = _5;
	    }
	}
	else
	{
	    if(color_func_.is_constant())
	    {
		boxsize_x_func_ = _3;
		boxsize_y_func_ = _3;
	    }
	    else
	    {
		boxsize_x_func_ = _4;
		boxsize_y_func_ = _4;
	    }
	}
	return *this;
    }

    csboxes &csboxes::xsize_min(double min)
    {
	xsize_min_ = min;
	if(xsize_min_ != unset) xsize_min_fixed_ = true;
	else xsize_min_fixed_ = false;
	return *this;
    }
    csboxes &csboxes::xsize_max(double max)
    {
	xsize_max_ = max;
	if(xsize_max_ != unset) xsize_max_fixed_ = true;
	else xsize_max_fixed_ = false;
	return *this;
    }
    csboxes &csboxes::xsize_range(double min, double max)
    {
	xsize_min(min);
	xsize_max(max);
	return *this;
    }

    csboxes &csboxes::ysize_min(double min)
    {
	ysize_min_ = min;
	if(ysize_min_ != unset) ysize_min_fixed_ = true;
	else ysize_min_fixed_ = false;
	return *this;
    }
    csboxes &csboxes::ysize_max(double max)
    {
	ysize_max_ = max;
	if(ysize_max_ != unset) ysize_max_fixed_ = true;
	else ysize_max_fixed_ = false;
	return *this;
    }

    csboxes &csboxes::ysize_range(double min, double max)
    {
	ysize_min(min);
	ysize_max(max);
	return *this;
    }

    csboxes &csboxes::xsize_goal_range(double min, double max)
    {
	xsize_goal_min(min);
	xsize_goal_max(max);
	return *this; 
    }
    csboxes &csboxes::xsize_goal_min(double min)
    {
	if(min!=unset) normalize_xsize_ = true;
	if(min<0)
	{
	    warning::print("Negative number in xsize_goal_min. Setting value to 0");
	    min = 0;
	}
	xsize_goal_min_ = min;
	return *this;
    }
    csboxes &csboxes::xsize_goal_max(double max)
    {
	if(max!=unset) normalize_xsize_ = true;
	xsize_goal_max_ = max;
	return *this;
    }

    csboxes &csboxes::ysize_goal_range(double min, double max)
    {
	ysize_goal_min(min);
	ysize_goal_max(max);
	return *this; 
    }
    csboxes &csboxes::ysize_goal_min(double min)
    {
	if(min!=unset) normalize_ysize_ = true;
	if(min<0)
	{
	    warning::print("Negative number in ysize_goal_min. Setting value to 0");
	    min = 0;
	}
	ysize_goal_min_ = min;
	return *this;
    }
    csboxes &csboxes::ysize_goal_max(double max)
    {
	if(max!=unset) normalize_ysize_ = true;
	ysize_goal_max_ = max;
	return *this;
    }

    csboxes &csboxes::dx(double v)
    {
	cell_dx_ = v;
	if(cell_dx_ == unset) cell_dx_fixed_ = false;
	else cell_dx_fixed_ = true;
	return *this;
    }

    csboxes &csboxes::dy(double v)
    {
	cell_dy_ = v;
	if(cell_dy_ == unset) cell_dy_fixed_ = false;
	else cell_dy_fixed_ = true;
	return *this;
    }

    void csboxes::prepare_for_draw(plottable *g, frame *f, int count)
    {
	if(global::debug>0) cout<<"[blop] [csboxes] prepare_for_draw("<<g<<","<<f<<","<<count<<")"<<endl;
	if(count==1)
	{
	    if(frame_foreground_) f->foreground(true);
	    if(grid_foreground_) f->grid_foreground(true);
	    
	    axis *xaxis =
		(g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	    axis *yaxis =
		(g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());
	    
	    calc(g,xaxis,yaxis);

	    if(colorlegend_)
	    {
		colorlegend_->clear_tics();
		if(draw_colorlegend_<2) colorlegend_->off();
		if(color_title_.str()!="")
		{
		    var new_title = colorlegend_->title();
		    if(new_title.str() != "") new_title &= ", ";
		    new_title &= color_title_;
		    colorlegend_->title(new_title);
		}
	    }
	}

	if(count==2)
	{
	    // Let the colorlegend calculate the unified range for all
	    // of its graphs
	    if(colorlegend_)
	    {
		colorlegend_->calculate_tics();
		if(draw_colorlegend_>=2) colorlegend_->on();  // if any of the graphs have legend on...
	    }
	}
	if(global::debug>0) cout<<"[blop] [csboxes] prepare_for_draw finished"<<endl;
    }

    csboxes::csboxes(double min, double max)
    {
	init_csboxes_();
	color_range(min,max);
	xsize_range(min,max);
	ysize_range(min,max);
    }

    csboxes::csboxes(const color &startcolor, const color &endcolor)
    {
	init_csboxes_();
	color_range(startcolor,endcolor);
    }

    csboxes::csboxes(const csboxes &rhs) : cbox_base(rhs)
    {
	boxsize_x_func_ = rhs.boxsize_x_func_;
	boxsize_y_func_ = rhs.boxsize_y_func_;
	color_func_     = rhs.color_func_;
	cell_dx_ = rhs.cell_dx_;
	cell_dy_ = rhs.cell_dy_;
	cell_dx_fixed_ = rhs.cell_dx_fixed_;
	cell_dy_fixed_ = rhs.cell_dy_fixed_;
	normalize_xsize_ = rhs.normalize_xsize_;
	normalize_ysize_ = rhs.normalize_ysize_;
	xsize_min_ = rhs.xsize_min_;
	xsize_max_ = rhs.xsize_max_;
	ysize_min_ = rhs.ysize_min_;
	ysize_max_ = rhs.ysize_max_;
	xsize_min_fixed_ = rhs.xsize_min_fixed_;
	xsize_max_fixed_ = rhs.xsize_max_fixed_;
	ysize_min_fixed_ = rhs.ysize_min_fixed_;
	ysize_max_fixed_ = rhs.ysize_max_fixed_;
	xsize_goal_min_ = rhs.xsize_goal_min_;
	xsize_goal_max_ = rhs.xsize_goal_max_;
	ysize_goal_min_ = rhs.ysize_goal_min_;
	ysize_goal_max_ = rhs.ysize_goal_max_;
	skip_outrange_x_ = rhs.skip_outrange_x_;
	skip_outrange_y_ = rhs.skip_outrange_y_;
	skip_outrange_color_ = rhs.skip_outrange_color_;
	fill_ = true;
	frame_foreground_ = rhs.frame_foreground_;
	grid_foreground_  = rhs.grid_foreground_;
    }

    csboxes::csboxes()
    {
	init_csboxes_();
    }

    void csboxes::init_csboxes_()
    {
	boxsize_x_func_ = _4;
	boxsize_y_func_ = _5;
	color_func_     = _3;
	cell_dx_ = unset;
	cell_dy_ = unset;
	cell_dx_fixed_ = false;
	cell_dy_fixed_ = false;
	normalize_xsize_ = true;
	normalize_ysize_ = true;
	xsize_min_ = unset;
	xsize_max_ = unset;
	ysize_min_ = unset;
	ysize_max_ = unset;
	xsize_min_fixed_ = false;
	xsize_max_fixed_ = false;
	ysize_min_fixed_ = false;
	ysize_max_fixed_ = false;
	skip_outrange_x_ = true;
	skip_outrange_y_ = true;
	skip_outrange_color_ = true;
	xsize_goal_min_ = 0;
	xsize_goal_max_ = 1;
	ysize_goal_min_ = 0;
	ysize_goal_max_ = 1;
	fill_ = true;
	frame_foreground_ = default_frame_foreground_;
	grid_foreground_  = default_grid_foreground_;
    }


    void csboxes::get_functions(plottable *g,
				function &getx, 
				function &gety,
				function &getboxsizex,
				function &getboxsizey,
				function &getcolor)
    {
	getx = get_x(g);
	gety = get_y(g);

	if(color_func_.initialized()) getcolor = color_func_;
	else getcolor = _3;

	if(boxsize_x_func_.initialized()) getboxsizex = boxsize_x_func_;
	else getboxsizex = _4;
	if(getboxsizex.equals(_4) && g->columns() < 4) getboxsizex = _3;

	if(boxsize_y_func_.initialized()) getboxsizey = boxsize_y_func_;
	else getboxsizey = _5;

	if(getboxsizey.nargs() > g->columns()) getboxsizey = getboxsizex;
    }

    void csboxes::set_ranges(plottable *g,axis *x,axis *y)
    {
	// prepare_for_draw is already called, so everything is set up now.
	x->autoextend_min_soft(false);
	x->autoextend_max_soft(false);
	y->autoextend_min_soft(false);
	y->autoextend_max_soft(false);
    }

    void csboxes::calc_box_x_(double x, double dx, double *x1, double *x2)
    {
	const double a = (xsize_goal_max_-xsize_goal_min_)/(xsize_max_-xsize_min_);
	const double b = (xsize_goal_min_*(xsize_max_-xsize_min_)-xsize_min_*(xsize_goal_max_-xsize_goal_min_))/(xsize_max_-xsize_min_);

	*x1 = x - 
	    (normalize_xsize_ ? (a*dx+b)*cell_dx_ : dx )/2;
	*x2 = x + 
	    (normalize_xsize_ ? (a*dx+b)*cell_dx_ : dx )/2;
    }

    void csboxes::calc_box_y_(double y, double dy, double *y1, double *y2)
    {
	const double a = (ysize_goal_max_-ysize_goal_min_)/(ysize_max_-ysize_min_);
	const double b = (ysize_goal_min_*(ysize_max_-ysize_min_)-ysize_min_*(ysize_goal_max_-ysize_goal_min_))/(ysize_max_-ysize_min_);
	*y1 = y - 
	    (normalize_ysize_ ? (a*dy+b)*cell_dy_ : dy )/2;
	*y2 = y + 
	    (normalize_ysize_ ? (a*dy+b)*cell_dy_ : dy )/2;
    }

    void csboxes::calc(plottable *g,axis *xaxis,axis *yaxis)
	TRY
    {
	if(!g) err("Plottable to draw not set");

	function getx, gety, getboxsizex, getboxsizey, getcolor;
	get_functions(g,getx,gety,getboxsizex,getboxsizey,getcolor);

	// Set ranges:
	if(g->xmin()!=unset) xaxis->extend_range(g->xmin());
	if(g->xmax()!=unset) xaxis->extend_range(g->xmax());
	if(g->ymin()!=unset) yaxis->extend_range(g->ymin());
	if(g->ymax()!=unset) yaxis->extend_range(g->ymax());
	if ( g->xmin()==unset || g->xmax()==unset || g->ymin()==unset || g->ymax()==unset ||
	     (!getcolor.is_constant() && (!color_min_fixed_ || !color_max_fixed_))  ||
	     (!getboxsizex.is_constant() && normalize_xsize_ && (!xsize_min_fixed_ || !xsize_max_fixed_)) ||
	     (!getboxsizey.is_constant() && normalize_ysize_ && (!ysize_min_fixed_ || !ysize_max_fixed_)) ||
	     (normalize_xsize_ && !cell_dx_fixed_) ||
	     (normalize_ysize_ && !cell_dy_fixed_) )
	{
	    double xmin=unset;
	    double xmax=unset;
	    double ymin=unset;
	    double ymax=unset;

	    if(!color_min_fixed_) color_min_ = unset;
	    if(!color_max_fixed_) color_max_ = unset;
	    if(!cell_dx_fixed_) cell_dx_ = unset;
	    if(!cell_dy_fixed_) cell_dy_ = unset;
	    if(!xsize_min_fixed_) xsize_min_ = 0;
	    if(!xsize_max_fixed_) xsize_max_ = unset;
	    if(!ysize_min_fixed_) ysize_min_ = 0;
	    if(!ysize_max_fixed_) ysize_max_ = unset;
	    
	    map<double,bool> xvalues,yvalues;
	    double xlast=unset, ylast=unset;

	    for(plottable::size_type i=0; i < g->size(); ++i)
	    {
		const var xorig_var = getx.eval(*(g->get(i)));
		const var yorig_var = gety.eval(*(g->get(i)));
		const var boxsize_x_var = getboxsizex.eval(*(g->get(i)));
		const var boxsize_y_var = getboxsizey.eval(*(g->get(i)));
		const var color_var     = getcolor.eval(*(g->get(i)));

		if ( ignore::it(xorig_var) || ignore::it(yorig_var) ||
		     ignore::it(boxsize_x_var) || ignore::it(boxsize_y_var) ||
		     ignore::it(color_var) )
		    continue;

		if ( xaxis->logscale() && xorig_var.dbl()<=0.0 ) continue;
		if ( yaxis->logscale() && yorig_var.dbl()<=0.0 ) continue;
		if ( color_logscale_    && color_var.dbl()<=0.0 ) continue;
		// check and skip points which are outside of g->xmin(), etc...

//		if ( xmin==unset ) xmin = xorig_var.dbl();
//		if ( xmax==unset ) xmax = xorig_var.dbl();
//		if ( ymin==unset ) ymin = yorig_var.dbl();
//		if ( ymax==unset ) ymax = yorig_var.dbl();

		if(xmin==unset || xorig_var.dbl() < xmin) xmin = xorig_var.dbl();
		if(xmax==unset || xorig_var.dbl() > xmax) xmax = xorig_var.dbl();
		if(ymin==unset || yorig_var.dbl() < ymin) ymin = yorig_var.dbl();
		if(ymax==unset || yorig_var.dbl() > ymax) ymax = yorig_var.dbl();

		if( color_min_fixed_ && color_var.dbl()<color_min_ ) continue;
		if( color_max_fixed_ && color_var.dbl()>color_max_ ) continue;
		
		if(!color_min_fixed_)
		{
		    if(color_min_ == unset) color_min_ = color_var.dbl();
		    if(color_var.dbl() < color_min_) color_min_ = color_var.dbl();
		}
		if(!color_max_fixed_)
		{
		    if(color_max_ == unset) color_max_ = color_var.dbl();
		    if(color_var.dbl() > color_max_) color_max_ = color_var.dbl();
		}

		if(!xsize_min_fixed_)
		{
		    if(xsize_min_ == unset) xsize_min_ = boxsize_x_var.dbl();
		    if(boxsize_x_var.dbl() < xsize_min_) xsize_min_ = boxsize_x_var.dbl();
		}
		if(!xsize_max_fixed_)
		{
		    if(xsize_max_ == unset) xsize_max_ = boxsize_x_var.dbl();
		    if(boxsize_x_var.dbl() > xsize_max_) xsize_max_ = boxsize_x_var.dbl();
		}
		if(!ysize_min_fixed_)
		{
		    if(ysize_min_ == unset) ysize_min_ = boxsize_y_var.dbl();
		    if(boxsize_y_var.dbl() < ysize_min_) ysize_min_ = boxsize_y_var.dbl();
		}
		if(!ysize_max_fixed_)
		{
		    if(ysize_max_ == unset) ysize_max_ = boxsize_y_var.dbl();
		    if(boxsize_y_var.dbl() > ysize_max_) ysize_max_ = boxsize_y_var.dbl();
		}

		if(!cell_dx_fixed_)
		{
		    if(g->ordered()>=2)
		    {
			if(xlast!=unset && xlast!=xorig_var.dbl() && (cell_dx_==unset || ::fabs(xorig_var.dbl()-xlast)<cell_dx_))
			    cell_dx_ = ::fabs(xorig_var.dbl()-xlast);
		    }
		    else xvalues[xorig_var.dbl()] = true;
		}
		if(!cell_dy_fixed_)
		{
		    if(g->ordered()>=2)
		    {
			if(ylast!=unset && ylast!=yorig_var.dbl() && (cell_dy_==unset || ::fabs(yorig_var.dbl()-ylast)<cell_dy_))
			    cell_dy_ = ::fabs(yorig_var.dbl()-ylast);
		    }
		    else yvalues[yorig_var.dbl()] = true;
		}
		xlast = xorig_var.dbl();
		ylast = yorig_var.dbl();
	    }

	    if(g->xmin()==unset && xmin!=unset) xaxis->extend_range(xmin);
	    if(g->xmax()==unset && xmax!=unset) xaxis->extend_range(xmax);
	    if(g->ymin()==unset && ymin!=unset) yaxis->extend_range(ymin);
	    if(g->ymax()==unset && ymax!=unset) yaxis->extend_range(ymax);

	    if(cell_dx_==unset)
	    {
		for(map<double,bool>::iterator i=xvalues.begin(); i!=xvalues.end(); ++i)
		{
		    if(i!=xvalues.begin())
		    {
			map<double,bool>::iterator prev=i;
			--prev;
			const double diff = ::fabs((*prev).first-(*i).first);
			if(cell_dx_==unset || diff < cell_dx_) cell_dx_ = diff;
		    }
		}
	    }
	    if(cell_dy_==unset)
	    {
		for(map<double,bool>::iterator i=yvalues.begin(); i!=yvalues.end(); ++i)
		{
		    if(i!=yvalues.begin())
		    {
			map<double,bool>::iterator prev=i;
			--prev;
			const double diff = ::fabs((*prev).first-(*i).first);
			if(cell_dy_==unset || diff < cell_dy_) cell_dy_ = diff;
		    }
		}
	    }

	    for(plottable::size_type i=0; i < g->size(); ++i)
	    {
		const var xorig_var = getx.eval(*(g->get(i)));
		const var yorig_var = gety.eval(*(g->get(i)));
		const var boxsize_x_var = getboxsizex.eval(*(g->get(i)));
		const var boxsize_y_var = getboxsizey.eval(*(g->get(i)));
		const var color_var     = getcolor.eval(*(g->get(i)));

		if(ignore::it(xorig_var) || ignore::it(yorig_var) ||
		   ignore::it(boxsize_x_var) || ignore::it(boxsize_y_var) ||
		   ignore::it(color_var)) continue;

		double x1=0, x2=0, y1=0, y2=0;
		calc_box_x_(xorig_var.dbl(), boxsize_x_var.dbl(), &x1, &x2);
		calc_box_y_(yorig_var.dbl(), boxsize_y_var.dbl(), &y1, &y2);

		if(g->xmin()==unset) xaxis->extend_range(x1);
		if(g->xmax()==unset) xaxis->extend_range(x2);
		if(g->ymin()==unset) yaxis->extend_range(y1);
		if(g->ymax()==unset) yaxis->extend_range(y2);
	    }
            if(cell_dx_<(xmax-xmin)/300) warning::print("The x cellsize was most probably not determined correctly (not a regular grid, or not exactly equal numbers)");
            if(cell_dy_<(ymax-ymin)/300) warning::print("The y cellsize was most probably not determined correctly (not a regular grid, or not exactly equal numbers)");
	}
    }
    CATCH("csboxes::calc(...)")


    void csboxes::draw(plottable *g,frame *f,terminal *t)
    TRY
    {
	if(!g || g->empty()) return;

	if(!normalize_xsize_ && xsize_min_<0) warning::print("Requested xsize_min is negative, but xsize is not normalized. Negative values will be skipped",
							     "csboxes::draw");
	if(!normalize_ysize_ && ysize_min_<0) warning::print("Requested ysize_min is negative, but ysize is not normalized. Negative values will be skipped",
							     "csboxes::draw");

	t->reset_transformation();

	axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());

	double xmin = xaxis->min();                                         
	if(g->xmin()!=unset && g->xmin()>xmin) xmin = g->xmin();
	double xmax = xaxis->max();
	if(g->xmax()!=unset && g->xmax()<xmax) xmax = g->xmax();

	double ymin = yaxis->min();
	if(g->ymin()!=unset && g->ymin()>ymin) ymin = g->ymin();
	double ymax = yaxis->max();
	if(g->ymax()!=unset && g->ymax()<ymax) ymax = g->ymax();

	vector<terminal::coord> cc;

	function getx, gety, getboxsizex, getboxsizey, getcolor;
	get_functions(g,getx,gety,getboxsizex,getboxsizey,getcolor);

	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    const var xorig_var     = getx.eval(*(g->get(i)));
	    const var yorig_var     = gety.eval(*(g->get(i)));
	    const var boxsize_x_var = getboxsizex.eval(*(g->get(i)));
	    const var boxsize_y_var = getboxsizey.eval(*(g->get(i)));
	    const var color_var     = getcolor.eval(*(g->get(i)));

	    if(ignore::it(xorig_var) || ignore::it(yorig_var) ||
	       ignore::it(boxsize_x_var) || ignore::it(boxsize_y_var) ||
	       ignore::it(color_var)) continue;

	    // skip points outside of the limits.
	    if( (boxsize_x_var.dbl()<xsize_min_ || xsize_max_<boxsize_x_var.dbl()) && skip_outrange_x_) continue;
	    if( (boxsize_y_var.dbl()<ysize_min_ || ysize_max_<boxsize_y_var.dbl()) && skip_outrange_y_) continue;
	    if( (color_var.dbl()<color_min_ || color_max_<color_var.dbl()) && skip_outrange_color_) continue;

	    double x1=0, x2=0, y1=0, y2=0;
	    calc_box_x_(xorig_var.dbl(), boxsize_x_var.dbl(), &x1, &x2);
	    calc_box_y_(yorig_var.dbl(), boxsize_y_var.dbl(), &y1, &y2);

	    // skip points which have 'negative' size, since they would just show up
	    // as normal boxes, confusing the user

	    if(x1 >= xorig_var.dbl() || y1 >= yorig_var.dbl())
	    {
		warning::print("Skipping undersized box","csboxes::draw");
		continue;
	    }

	    if(xaxis->logscale() && x1<=0) continue;
	    if(xaxis->logscale() && x2<=0) continue;
	    if(yaxis->logscale() && y1<=0) continue;
	    if(yaxis->logscale() && y2<=0) continue;
	    if(x2<xmin) continue;
	    if(x1>xmax) continue;
	    if(y2<ymin) continue;
	    if(y1>ymax) continue;

	    if(!color_func_.is_constant() && color_logscale_ && color_var.dbl()<=0.0) continue;

	    double x1_mapped = xaxis->map_point(x1);
	    double x2_mapped = xaxis->map_point(x2);
	    double y1_mapped = yaxis->map_point(y1);
	    double y2_mapped = yaxis->map_point(y2);

            const double xmin_mapped = xaxis->map_point(xmin);
            const double xmax_mapped = xaxis->map_point(xmax);
            const double ymin_mapped = yaxis->map_point(ymin);
            const double ymax_mapped = yaxis->map_point(ymax);

            if(x1_mapped<xmin_mapped) x1_mapped = xmin_mapped;
            if(x2_mapped>xmax_mapped) x2_mapped = xmax_mapped;
            if(y1_mapped<ymin_mapped) y1_mapped = ymin_mapped;
            if(y2_mapped>ymax_mapped) y2_mapped = ymax_mapped;
            if(x1_mapped >= x2_mapped) continue;
            if(y1_mapped >= y2_mapped) continue;

	    //if(-0.01<=x1_mapped && x1_mapped<=1    &&
	    //   0    <=x2_mapped && x2_mapped<=1.01 &&
	    //   -0.01<=y1_mapped && y1_mapped<=1    &&
	    //   0    <=y2_mapped && y2_mapped<=1.01)
	    {
                

                //const double dx = x2_mapped-x1_mapped;
                //x1_mapped -= dx*cell_margin_;
                //x2_mapped += dx*cell_margin_;
                //const double dy = y2_mapped-y1_mapped;
                //y1_mapped -= dy*cell_margin_;
                //y2_mapped += dy*cell_margin_;

		cc.clear();
		cc.push_back(terminal::coord(terminal::id(x1_mapped,1),terminal::id(y1_mapped,2)));
		cc.push_back(terminal::coord(terminal::id(x1_mapped,1),terminal::id(y2_mapped,2)));
		cc.push_back(terminal::coord(terminal::id(x2_mapped,1),terminal::id(y2_mapped,2)));
		cc.push_back(terminal::coord(terminal::id(x2_mapped,1),terminal::id(y1_mapped,2)));
		cc.push_back(terminal::coord(terminal::id(x1_mapped,1),terminal::id(y1_mapped,2)));
		if(!color_func_.is_constant())
		{
		    double minimum = min();
		    double maximum = max();
//                    cerr<<"min1 = "<<minimum<<endl;
//                    cerr<<"max1 = "<<maximum<<endl;
		    if(colorlegend_ && draw_colorlegend_>=2)
		    {
			minimum = colorlegend_->min();
			maximum = colorlegend_->max();
		    }
		    const color the_color = map_color(color_var.dbl(),minimum,maximum);
//                    cerr<<"min = "<<minimum<<endl;
//                    cerr<<"max = "<<maximum<<endl;
//                    cerr<<"val = "<<color_var.dbl()<<endl;
//                    cerr<<"col = "<<the_color<<endl;
//                    cerr<<endl;
		    t->set_color(the_color);
		}
		else
		{
		    if(fill_) t->set_color(g->fillcolor());
		    else      t->set_color(g->linecolor());
		}
		if(fill_) t->fill_polygon(cc);
		else
		{
		    t->set_linewidth(g->linewidth().termspecific_id());
		    t->set_linestyle(g->linestyle());
		    t->draw_lines(cc);
		}
	    }
	}
    }
    CATCH("csboxes::draw(...)")

    graph_drawer *csboxes::clone() const
    {
	return new csboxes(*this);
    }

    // ---------- cboxes  -------------------------------------------------------------

    void cboxes::init_()
    {
	boxsize_x_func_ = 1;
	boxsize_y_func_ = 1;
	color_func_ = _3;
	xsize_range(0,1);
	ysize_range(0,1);
        //cell_margin_ = 0.02;
    }
    cboxes::cboxes() {init_();}
    cboxes::cboxes(const cboxes &rhs) : csboxes(rhs) {init_();}
    cboxes::cboxes(const csboxes &rhs) : csboxes(rhs) {init_();}
    cboxes::cboxes(double min, double max) { init_(); color_range(min,max); }
    cboxes::cboxes(const color &start, const color &end)
    {
	init_();
	color_range(start,end);
    }
    cboxes::cboxes(void *p)
    {
	init_();
	operator()(p);
    }
    cboxes::cboxes(color (*p)(double,double,double))
    {
	init_();
	operator()(p);
    }
    cboxes &cboxes::operator()(double min,double max)
    {
	color_range(min,max);
	return *this;
    }
    cboxes &cboxes::operator()(const color &start, const color &end)
    {
	color_range(start,end);
	return *this;
    }
    cboxes &cboxes::operator()(void *p) {map_function(p); return *this; }
    cboxes &cboxes::operator()(color (*p)(double,double,double)) {map_function(p); return *this; }

    // ---------- sboxes  -------------------------------------------------------------

    void sboxes::init_()
    {
	color_func_ = 1;
	draw_colorlegend_ = 0;
	boxsize_x_func_ = _3;
	boxsize_y_func_ = _3;
	fill_ = false;
    }
    sboxes::sboxes() {init_();}
    sboxes::sboxes(const sboxes &rhs) : csboxes(rhs) {init_();}
    sboxes::sboxes(const csboxes &rhs) : csboxes(rhs) {init_();}
    sboxes::sboxes(double min, double max) {init_(); xsize_range(min,max); ysize_range(min,max);}
    sboxes &sboxes::operator()(double min, double max)
    {
	xsize_range(min,max);
	ysize_range(min,max);
	return *this;
    }
    void sboxes::draw_sample(const length &x,
			     const length &y,
			     const length &size,
			     const plottable *s,
			     terminal *t)
    {
	vector<terminal::coord> cc;
	length l1 = !x - 0.5*!size;
	length l2 = !x + 0.5*!size;
	length l3 = !y - 0.5*0.8*EM;
	length l4 = !y + 0.5*0.8*EM;
	l1.specialize(t);
	l2.specialize(t);
	l3.specialize(t);
	l4.specialize(t);
	cc.push_back(terminal::coord(l1.termspecific_id(),l3.termspecific_id()));
	cc.push_back(terminal::coord(l2.termspecific_id(),l3.termspecific_id()));
	cc.push_back(terminal::coord(l2.termspecific_id(),l4.termspecific_id()));
	cc.push_back(terminal::coord(l1.termspecific_id(),l4.termspecific_id()));
	cc.push_back(terminal::coord(l1.termspecific_id(),l3.termspecific_id()));
	if(s->fill())
	{
	    t->set_color(s->fillcolor());
	    t->fill_polygon(cc);
	}
	else
	{
	    t->set_color(s->linecolor());
	    t->set_linewidth(s->linewidth().termspecific_id());
	    t->set_linestyle(s->linestyle());
	    t->draw_lines(cc);
	}
    }    

/*
// ----------------------  sboxes  ------------------------------------

void sboxes::prepare_for_draw(plottable *g, frame *f, int count)
{
if(count==1)
{
axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());
g->linewidth().register_me();
calc(g,xaxis,yaxis);
}
}

sboxes::sboxes(const sboxes &rhs)
{
dx_ = rhs.dx_;
dy_ = rhs.dy_;
min_ = rhs.min_;
max_ = rhs.max_;
min_fixed_ = rhs.min_fixed_;
max_fixed_ = rhs.max_fixed_;
}

sboxes::sboxes(double minimum,double maximum)
{
dx_ = unset;
dy_ = unset;
min_ = minimum;
if(min_ == unset) min_fixed_ = false;
else min_fixed_ = true;
max_ = maximum;
if(max_ == unset) max_fixed_ = false;
else max_fixed_ = true;
}


void sboxes::set_ranges(plottable *g,axis *x,axis *y)
{
// As prepare_for_draw is already called, everything is set up.
x->autoextend_min_soft(false);
x->autoextend_max_soft(false);
y->autoextend_min_soft(false);
y->autoextend_max_soft(false);
}


void sboxes::calc(plottable *g,axis *xaxis,axis *yaxis)
TRY
{
if(!g) err("Plottable to draw not set");

function getx = get_x(g);
function gety = get_y(g);
function getdx = g->dx_hint();
function getdy = g->dy_hint();
function getz = g->z_hint();
if(!getz.initialized()) getz = _3;

const bool calculate_dx = (dx_ == unset);
const bool calculate_dy = (dy_ == unset);

if ( dx_==unset || dy_==unset )
{
if ( g->ordered() >= 2)
{
double xmin=unset;
double xmax=unset;
double ymin=unset;
double ymax=unset;
int xbins=0;
int ybins=0;
int xdir=false;
int ydir=false;
double xfirst=unset;
double yfirst=unset;
double xsec=unset;
double ysec=unset;
for ( unsigned int i=0 ; i<g->size() ; ++i )
{
const var x_var=getx.eval(*(g->get(i)));
const var y_var=gety.eval(*(g->get(i)));
if( ignore::it(x_var) || ignore::it(y_var) ) continue;
const double x=x_var.dbl();
const double y=y_var.dbl();
if ( xfirst!=unset && xsec==unset && yfirst!=unset && ysec==unset )
{
xsec=x; 
ysec=y;
if ( xsec>xfirst ) { xdir=true; ydir=false; ++xbins; }
else if ( ysec>yfirst ) { ydir=true; xdir=false; ++ybins; }
else
{
warning::print("Cannot determine grid direction. "
"Assuming to be ordered in (y,x).",
"void cboxes::calc");
ydir=true;
xdir=false;
++ybins;
}
}
if ( xfirst==unset && yfirst==unset )
{
xmin=x; xmax=x; xfirst=x;
ymin=y; ymax=y; yfirst=y;
}
if ( x<xmin ) xmin=x;
else if ( xmax<x ) xmax=x;
if ( y<ymin ) ymin=y;
else if ( ymax<y ) ymax=y;
if ( xdir && !(y>yfirst) ) ++xbins;
else if ( ydir && !(x>xfirst) ) ++ybins;
}
if ( xdir ) ybins=g->size()/xbins;
if ( ydir ) xbins=g->size()/ybins;
if ( calculate_dx ) dx_=(xmax-xmin)/(xbins-1);
if ( calculate_dy ) dy_=(ymax-ymin)/(ybins-1);
}
else if ( getdx.initialized() && getdy.initialized() )
{
for(plottable::size_type i=0; i<g->size(); ++i)
{
const var x_var=getx.eval(*(g->get(i)));
const var y_var=gety.eval(*(g->get(i)));
const var dx_var=getdx.eval(*(g->get(i)));
const var dy_var=getdy.eval(*(g->get(i)));
if ( ignore::it(x_var) || ignore::it(y_var) || ignore::it(dx_var) || ignore::it(dy_var) ) continue;
double dx = dx_var.dbl();
double dy = dy_var.dbl();
if(calculate_dx && (dx_==unset || dx < dx_)) dx_ = dx;
if(calculate_dy && (dy_==unset || dy < dy_)) dy_ = dy;
}
}
else
{
map<double,bool> xvalues,yvalues;
for(plottable::size_type i1=0; i1<g->size(); ++i1)
{
const var xorig1_var = getx.eval(*(g->get(i1)));
const var yorig1_var = gety.eval(*(g->get(i1)));
if(ignore::it(xorig1_var) || ignore::it(yorig1_var)) continue;
double xorig1 = xorig1_var.dbl();
double yorig1 = yorig1_var.dbl();
for(plottable::size_type i2=i1+1; i2<g->size(); ++i2)
{
const var xorig2_var = getx.eval(*(g->get(i2)));
const var yorig2_var = gety.eval(*(g->get(i2)));
if(ignore::it(xorig2_var) || ignore::it(yorig2_var)) continue;
double xorig2 = xorig2_var.dbl();
double yorig2 = yorig2_var.dbl();
		    
if(calculate_dx && xorig1!=xorig2)
{
if(xorig1 > xorig2) swap(xorig1,xorig2);
if(dx_ == unset || (xorig2-xorig1) < dx_) dx_ = (xorig2-xorig1);
}
		    
if(calculate_dy && yorig1!=yorig2)
{
if(yorig1 > yorig2) swap(yorig1,yorig2);
if(dy_ == unset || (yorig2-yorig1) < dy_) dy_ = (yorig2-yorig1);
}
}
}
}
}

if(dx_ == unset || dx_ <=0)
{
warning::print("Couldn't figure out boxwidth, using dx=1","sboxes::calc(...)");
dx_ = 1;
}
if(dy_ == unset || dy_ <=0)
{
warning::print("Couldn't figure out boxheight, using dx=1","sboxes::calc(...)");
dy_ = 1;
}

// Set ranges:
if(g->xmin()!=unset) xaxis->extend_range(g->xmin());
if(g->xmax()!=unset) xaxis->extend_range(g->xmax());
if(g->ymin()!=unset) yaxis->extend_range(g->ymin());
if(g->ymax()!=unset) yaxis->extend_range(g->ymax());
if ( g->xmin()==unset || g->xmax()==unset || g->ymin()==unset || g->ymax()==unset || min_fixed_==false || max_fixed_==false )
{
double xmin=unset;
double xmax=unset;
double ymin=unset;
double ymax=unset;
double vmin=unset;
double vmax=unset;
for(plottable::size_type i=0; i < g->size(); ++i)
{
const var xorig_var = getx.eval(*(g->get(i)));
const var yorig_var = gety.eval(*(g->get(i)));
if ( ignore::it(xorig_var) || ignore::it(yorig_var) ) continue;
const double xorig=xorig_var.dbl();
const double yorig=yorig_var.dbl();
if ( (xaxis->logscale() && !(xorig-0.5*dx_>0.0)) || (yaxis->logscale() && !(yorig-0.5*dy_>0.0)) ) continue;
if ( (g->xmin()!=unset && xorig+0.5*dx_<g->xmin()) || (g->xmax()!=unset && xorig-0.5*dx_>g->xmax()) || (g->ymin()!=unset && yorig+0.5*dy_<g->ymin()) || (g->ymax()!=unset && yorig-0.5*dy_>g->ymax()) ) continue;
if ( xmin==unset && xmax==unset ) { xmin=xorig-0.5*dx_; xmax=xorig+0.5*dx_; }
if ( ymin==unset && ymax==unset ) { ymin=yorig-0.5*dy_; ymax=yorig+0.5*dy_; }
if ( xorig-0.5*dx_<xmin ) xmin=xorig-0.5*dx_;
else if ( xorig+0.5*dx_>xmax ) xmax=xorig+0.5*dx_;
if ( yorig-0.5*dy_<ymin ) ymin=yorig-0.5*dy_;
else if ( yorig+0.5*dy_>ymax ) ymax=yorig+0.5*dy_;

const var v_var = getz.eval(*(g->get(i)));
if ( ignore::it(v_var) ) continue;
const double v = v_var.dbl();
if((min_fixed_ && v<min_) || (max_fixed_ && v>max_)) continue;
if(vmin==unset && vmax==unset) { vmin=v; vmax=v; }
if( v > vmax) vmax = v;
else if( v < vmin) vmin = v;
}
if(g->xmin()==unset && xmin!=unset) xaxis->extend_range(xmin);
if(g->xmax()==unset && xmax!=unset) xaxis->extend_range(xmax);
if(g->ymin()==unset && ymin!=unset) yaxis->extend_range(ymin);
if(g->ymax()==unset && ymax!=unset) yaxis->extend_range(ymax);
if(min_fixed_==false && vmin!=unset) min_=vmin;
if(max_fixed_==false && vmax!=unset) max_=vmax;
}
}
CATCH("sboxes::calc(...)")

void sboxes::draw(plottable *g,frame *f,terminal *t)
TRY
{
if(!g || g->empty())  return;

function getx = get_x(g);
function gety = get_y(g);
function getz = g->z_hint();
if(!getz.initialized()) getz = _3;

if(g->fill())
{
t->set_color(g->fillcolor());
}
else
{
t->set_linewidth(g->linewidth().termspecific_id());
t->set_linestyle(g->linestyle());
t->set_color(g->linecolor());
}

t->reset_transformation();

axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());

double xmin = xaxis->min();
if(g->xmin()!=unset && g->xmin()>xmin) xmin = g->xmin();
double xmax = xaxis->max();
if(g->xmax()!=unset && g->xmax()<xmax) xmax = g->xmax();

double ymin = yaxis->min();
if(g->ymin()!=unset && g->ymin()>ymin) ymin = g->ymin();
double ymax = yaxis->max();
if(g->ymax()!=unset && g->ymax()<ymax) ymax = g->ymax();

vector<terminal::coord> cc;

double maxi;
double mini;

{
mini =  min_ - 0.1*(max_-min_);
maxi =  max_ + 0.1*(max_-min_);
}
if(min_fixed_) mini = min_;
if(max_fixed_) maxi = max_;

for(plottable::size_type i=0; i<g->size(); ++i)
{	
const var xorig_var = getx.eval(*(g->get(i)));
const var yorig_var = gety.eval(*(g->get(i)));
const var zorig_var = getz.eval(*(g->get(i)));
if(ignore::it(xorig_var) || ignore::it(yorig_var) || ignore::it(zorig_var)) continue;
const double xorig = xorig_var.dbl();
const double yorig = yorig_var.dbl();
const double zorig = zorig_var.dbl();

double dx = (zorig-mini)/(maxi-mini);
double dy = (zorig-mini)/(maxi-mini);
if ( dx<0.0 ) dx=0.0;
if ( dx>1.0 ) dx=1.0;
if ( dy<0.0 ) dy=0.0;
if ( dy>1.0 ) dy=1.0;
dx*=dx_;
dy*=dy_;

if((xaxis->logscale() && !(xorig-0.5*dx>0.0)) || (yaxis->logscale() && !(yorig-0.5*dy>0.0))) continue;
if(xorig+0.5*dx<xmin || xorig-0.5*dx>xmax || yorig+0.5*dy<ymin || yorig-0.5*dy>ymax) continue;

double x1 = xaxis->map_point(xorig - dx/2);
double x2 = xaxis->map_point(xorig + dx/2);
double y1 = yaxis->map_point(yorig - dy/2);
double y2 = yaxis->map_point(yorig + dy/2);

if(0<=x1 && x1<=1 &&
0<=x2 && x2<=1 &&
0<=y1 && y1<=1 &&
0<=y2 && y2<=1)
{
cc.clear();
cc.push_back(terminal::coord(terminal::id(x1,1),terminal::id(y1,2)));
cc.push_back(terminal::coord(terminal::id(x1,1),terminal::id(y2,2)));
cc.push_back(terminal::coord(terminal::id(x2,1),terminal::id(y2,2)));
cc.push_back(terminal::coord(terminal::id(x2,1),terminal::id(y1,2)));
cc.push_back(terminal::coord(terminal::id(x1,1),terminal::id(y1,2)));
if(g->fill()) t->fill_polygon(cc);
else t->draw_lines(cc);
}
}
}
CATCH("sboxes::draw(...)")

graph_drawer *sboxes::clone() const
{
return new sboxes(*this);
}

void sboxes::set_range_(double mini, double maxi)
{
min_ = mini;
max_ = maxi;
if(mini == unset) min_fixed_ = false;
else min_fixed_ = true;
if(maxi == unset) max_fixed_ = false;
else max_fixed_ = true;
}



*/

// ---------- band   --------------------------------------------------------------

function bands::get_x1_(plottable *g) const
{
    function result = g->x1_hint();
    if(!result.initialized()) result = _1;
    return result;
}
    function bands::get_x2_(plottable *g) const
    {
	function result = g->x2_hint();
	if(!result.initialized())
	{
	    if(x_) result = _2;
	    else return _1;
	}
	return result;
    }
    function bands::get_y1_(plottable *g) const
    {
	function result = g->y1_hint();
	if(!result.initialized())
	{
	    if(x_) result = _3;
	    else result = _2;
	}
	return result;
    }
    function bands::get_y2_(plottable *g) const
    {
	function result = g->y2_hint();
	if(!result.initialized())
	{
	    if(x_)
	    {
		if(y_) result = _4;
		else result = _3;
	    }
	    else
	    {
		if(y_) 	result = _3;
		else result = _2;
	    }
	}
	return result;
    }


    void bands::set_ranges(plottable *g,axis *xaxis,axis *yaxis)
	TRY
    {
	if(!g) err("Plottable to draw not set");
	if(!xaxis) err("Xaxis not set");
	if(!yaxis) err("Yaxis not set");
	if(g->empty()) return;

	if(g->xmin()!=unset) xaxis->extend_range(g->xmin());
	if(g->xmax()!=unset) xaxis->extend_range(g->xmax());
	if(g->ymin()!=unset) yaxis->extend_range(g->ymin());
	if(g->ymax()!=unset) yaxis->extend_range(g->ymax());

	if(g->xmin()!=unset && g->xmax()!=unset && g->ymin()!=unset && g->ymax()!=unset) return;

	function gety1 = get_y1_(g);
	function gety2 = get_y2_(g);
	function getx1 = get_x1_(g);
	function getx2 = get_x2_(g);

	double xmin=unset;
	double xmax=unset;
	double ymin=unset;
	double ymax=unset;

	for ( unsigned int i=0 ; i<g->size() ; ++i )
	{

	    const var x1orig_var=getx1.eval(*g->get(i));
	    if ( !ignore::it(x1orig_var) )
	    {
		const double x1orig=x1orig_var.dbl();
		if ( !(xaxis->logscale() && !(x1orig>0.0)) )
		{
		    if ( !((g->xmin()!=unset && x1orig<g->xmin()) || (g->xmax()!=unset && x1orig>g->xmax())) )
		    {
			if ( xmin==unset && xmax==unset ) { xmin=x1orig; xmax=x1orig; }
			if ( x1orig<xmin ) xmin=x1orig;
			else if ( x1orig>xmax ) xmax=x1orig;

			const var y1orig_var=gety1.eval(*g->get(i));
			if ( !ignore::it(y1orig_var) )
			{
			    const double y1orig=y1orig_var.dbl();
			    if ( !(yaxis->logscale() && !(y1orig>0.0)) )
			    {
				if ( !((g->ymin()!=unset && y1orig<g->ymin()) || (g->ymax()!=unset && y1orig>g->ymax())) )
				{
				    if ( ymin==unset && ymax==unset ) { ymin=y1orig; ymax=y1orig; }
				    if ( y1orig<ymin ) ymin=y1orig;
				    else if ( y1orig>ymax ) ymax=y1orig;
				}
			    }
			}
		    }
		}
	    }

	    const var x2orig_var=getx2.eval(*g->get(i));
	    if ( !ignore::it(x2orig_var) )
	    {
		const double x2orig=x2orig_var.dbl();
		if ( !(xaxis->logscale() && !(x2orig>0.0)) )
		{
		    if ( !((g->xmin()!=unset && x2orig<g->xmin()) || (g->xmax()!=unset && x2orig>g->xmax())) )
		    {
			if ( xmin==unset && xmax==unset ) { xmin=x2orig; xmax=x2orig; }
			if ( x2orig<xmin ) xmin=x2orig;
			else if ( x2orig>xmax ) xmax=x2orig;

			const var y2orig_var=gety2.eval(*g->get(i));
			if ( !ignore::it(y2orig_var) )
			{
			    const double y2orig=y2orig_var.dbl();
			    if ( !(yaxis->logscale() && !(y2orig>0.0)) )
			    {
				if ( !((g->ymin()!=unset && y2orig<g->ymin()) || (g->ymax()!=unset && y2orig>g->ymax())) )
				{
				    if ( ymin==unset && ymax==unset ) { ymin=y2orig; ymax=y2orig; }
				    if ( y2orig<ymin ) ymin=y2orig;
				    else if ( y2orig>ymax ) ymax=y2orig;
				}
			    }
			}
		    }
		}
	    }

	}

	if (g->xmin()==unset && xmin!=unset) xaxis->extend_range(xmin);
	if (g->xmax()==unset && xmax!=unset) xaxis->extend_range(xmax);
	if (g->ymin()==unset && ymin!=unset) yaxis->extend_range(ymin);
	if (g->ymax()==unset && ymax!=unset) yaxis->extend_range(ymax);
    }
    CATCH("bands::set_ranges()")

    void bands::draw(plottable *g,frame *f,terminal *t)
    TRY
    {
	if(!g) err("Plottable to draw not set");
	if(g->empty())  return;
	axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());
	if(xaxis == 0 || yaxis == 0) return;

	double xmin = xaxis->min();
	if(g->xmin() != unset && g->xmin() > xmin) xmin = g->xmin();
	double xmax = xaxis->max();
	if(g->xmax() != unset && g->xmax() < xmax) xmax = g->xmax();
	    
	double ymin = yaxis->min();
	if(g->ymin() != unset && g->ymin() > ymin) ymin = g->ymin();
	double ymax = yaxis->max();
	if(g->ymax() != unset && g->ymax() < ymax) ymax = g->ymax();

	function gety1 = get_y1_(g);
	function gety2 = get_y2_(g);
	function getx1 = get_x1_(g);
	function getx2 = get_x2_(g);

	vector<terminal::coord> cc;

	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    const var x1orig = getx1.eval(*(g->get(i)));
	    const var y1orig = gety1.eval(*(g->get(i)));
	    const var x2orig = getx2.eval(*(g->get(i)));
	    const var y2orig = gety2.eval(*(g->get(i)));
	    if(ignore::it(x1orig) || ignore::it(x2orig) ||
	       ignore::it(y1orig) || ignore::it(y2orig)) continue;
	    if((xaxis->logscale() && !(x1orig.dbl()>0.0)) ||
	       (xaxis->logscale() && !(x2orig.dbl()>0.0)) ||
	       (yaxis->logscale() && !(y1orig.dbl()>0.0)) ||
	       (yaxis->logscale() && !(y2orig.dbl()>0.0))) continue;

	    const double x = xaxis->map_point(x1orig.dbl());
	    const double y = yaxis->map_point(y1orig.dbl());

	    cc.push_back(terminal::coord(terminal::id(x,1),terminal::id(y,2)));
	}

	for(int i=g->size()-1; i >= 0; --i)
	{
	    const var x1orig = getx1.eval(*(g->get(i)));
	    const var y1orig = gety1.eval(*(g->get(i)));
	    const var x2orig = getx2.eval(*(g->get(i)));
	    const var y2orig = gety2.eval(*(g->get(i)));
	    if(ignore::it(x1orig) || ignore::it(x2orig) ||
	       ignore::it(y1orig) || ignore::it(y2orig)) continue;
	    if((xaxis->logscale() && !(x1orig.dbl()>0.0)) ||
	       (xaxis->logscale() && !(x2orig.dbl()>0.0)) ||
	       (yaxis->logscale() && !(y1orig.dbl()>0.0)) ||
	       (yaxis->logscale() && !(y2orig.dbl()>0.0))) continue;

	    double x = xaxis->map_point(x2orig.dbl());
	    double y = yaxis->map_point(y2orig.dbl());

	    cc.push_back(terminal::coord(terminal::id(x,1),terminal::id(y,2)));
	}

	t->set_color(g->fillcolor());
	bool termclip = t->clip(terminal::coord(terminal::id(xaxis->map_point(xmin),1),
						terminal::id(yaxis->map_point(ymin),2)),
				terminal::coord(terminal::id(xaxis->map_point(xmax),1),
						terminal::id(yaxis->map_point(ymax),2)));
	if(!termclip) warning::print("This terminal does not support clipping. "
				     "Graph drawn with band style might reach out of frame",
				     "bands::draw(...)");
	t->fill_polygon(cc);
	if(termclip) t->noclip();
    }
    CATCH("bands::draw(plottable *,frame *,terminal *)")

    void bands::draw_sample(const length &x,const length &y,
			    const length &size,
			    const plottable *g,terminal *t)
    {
	length x1 = -0.5*!size;
	length x2 =  0.5*!size;
	length y1 = -0.5*EX;
	length y2 = 0.5*EX;

	x1.specialize(t);
	x2.specialize(t);
	y1.specialize(t);
	y2.specialize(t);
	
	vector<terminal::coord> cc;

	cc.push_back(terminal::coord(x1.termspecific_id(),y1.termspecific_id()));
	cc.push_back(terminal::coord(x2.termspecific_id(),y1.termspecific_id()));
	cc.push_back(terminal::coord(x2.termspecific_id(),y2.termspecific_id()));
	cc.push_back(terminal::coord(x1.termspecific_id(),y2.termspecific_id()));
	cc.push_back(terminal::coord(x1.termspecific_id(),y1.termspecific_id()));
	t->translate(x.termspecific_id(),y.termspecific_id());
	t->set_color(g->fillcolor());
	t->fill_polygon(cc);
	t->reset_transformation();

    }

    graph_drawer *bands::clone() const
    {
	return new bands(*this);
    }


    // ---------- mosaic  -------------------------------------------------------------

    bool mosaic::default_frame_foreground_ = true;
    bool mosaic::default_grid_foreground_ = false;

    mosaic::mosaic(const mosaic &rhs)
	: cbox_base(rhs), params_to_xy_(rhs.params_to_xy_), fixdp_(rhs.fixdp_),
	  frame_foreground_(rhs.frame_foreground_), grid_foreground_(rhs.grid_foreground_)
    {
    }


    mosaic::mosaic(const function &f1)
	: params_to_xy_(f1), fixdp_(true),
	  frame_foreground_(mosaic::default_frame_foreground_),
	  grid_foreground_(mosaic::default_grid_foreground_)
    {
    }

    mosaic::mosaic(const function &f1, const function &f2)
	: params_to_xy_(f1,f2), fixdp_(true),
	  frame_foreground_(mosaic::default_frame_foreground_),
	  grid_foreground_(mosaic::default_grid_foreground_)
    {
    }

    void mosaic::set_ranges(plottable *g, axis *x, axis *y)
    {
	vector<var> params(2), xy(2);
	for(unsigned int i=0; i<g->size(); ++i)
	{
	    params[0] = (*g->get(i))[0];
	    params[1] = (*g->get(i))[1];
	    params_to_xy_.meval(params,xy);
	    x->extend_range(xy[0].dbl());
	    y->extend_range(xy[1].dbl());

	}
    }

    void mosaic::prepare_for_draw(plottable *g, frame *f, int count)
    {
	if(count==1)
	{
	    if(frame_foreground_) f->foreground(true);
	    double zmin=unset, zmax=unset;
	    for(unsigned int i=0; i<g->size(); ++i)
	    {
		const double z = (*g->get(i))[2];
		if(zmin == unset || z<zmin) zmin = z;
		if(zmax == unset || z>zmax) zmax = z;
	    }
	    if(!color_min_fixed_ && zmin!=unset) color_min_ = zmin;
	    if(!color_max_fixed_ && zmax!=unset) color_max_ = zmax;
	    if(colorlegend_)
	    {
		colorlegend_->clear_tics();
		if(draw_colorlegend_<2) colorlegend_->off();
		if(color_title_.str()!="")
		{
		    var new_title = colorlegend_->title();
		    if(new_title.str() != "") new_title &= ", ";
		    new_title &= color_title_;
		    colorlegend_->title(new_title);
		}
	    }
	}

	if(count==2)
	{
	    // Let the colorlegend calculate the unified range for all
	    // of its graphs
	    if(colorlegend_)
	    {
		colorlegend_->calculate_tics();
		if(draw_colorlegend_>=2) colorlegend_->on();  // if any of the graphs have legend on...
	    }		
	}
    }

    void mosaic::draw(plottable *g, frame *f, terminal *term)
    {
	if(!g || g->empty()) return;
	term->reset_transformation();
	axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());

	std::map<double,bool> p1map, p2map;
	double zmin = unset, zmax = unset;

	// first loop over the whole graph, and collect the values of the parameters
	for(unsigned int i=0; i<g->size(); ++i)
	{
	    p1map[(*g->get(i))[0].dbl()];
	    p2map[(*g->get(i))[1].dbl()];

	    const double z = (*g->get(i))[2];
	    if(zmin == unset || z < zmin) zmin = z;
	    if(zmax == unset || z > zmax) zmax = z;
	}

	if(color_min_fixed_) zmin = color_min_;
	if(color_max_fixed_) zmax = color_max_;

	// if we assume a fix stepsize in both parameters, then loop over
	// all neighbouring parameter values, and calculate the smallest
	// difference between them
	double dp1=unset, dp2=unset;
	if(fixdp_)
	{
	    {
		std::map<double,bool>::iterator i=p1map.begin();
		for(++i; i != p1map.end(); ++i)
		{
		    std::map<double,bool>::iterator prev = i; --prev;
		    const double diff = ::fabs((*i).first-(*prev).first);
		    if(dp1 == unset || diff<dp1) dp1=diff;
		}
	    }
	    {
		std::map<double,bool>::iterator i=p2map.begin();
		for(++i; i != p2map.end(); ++i)
		{
		    std::map<double,bool>::iterator prev = i; --prev;
		    const double diff = ::fabs((*i).first-(*prev).first);
		    if(dp2 == unset || diff<dp2) dp2=diff;
		}
	    }
	}

	std::vector<blop::var> param(2);
	std::vector<blop::var> xy   (2);
	for(unsigned int i=0; i<g->size(); ++i)
	{
	    const double z = (*g->get(i))[2];
	    if(z<zmin || zmax<z) continue;
	    param[0] = (*g->get(i))[0];
	    param[1] = (*g->get(i))[1];
	    
	    double p1_1, p1_2, p2_1, p2_2;

	    if(fixdp_)
	    {
		p1_1 = param[0].dbl()-0.5*dp1;
		p1_2 = param[0].dbl()+0.5*dp1;
		p2_1 = param[1].dbl()-0.5*dp2;
		p2_2 = param[1].dbl()+0.5*dp2;
	    }
	    else
	    {
		const std::map<double,bool>::iterator p1i = p1map.find(param[0].dbl());
		const std::map<double,bool>::iterator p2i = p2map.find(param[1].dbl());
		
		map<double,bool>::iterator other;
		
		if(p1i == p1map.begin())
		{
		    other = p1i; ++other; // get the next one
		    p1_1 = 1.5*(*p1i).first-0.5*(*other).first;
		}
		else
		{
		    other = p1i; --other;  // get the previous one
		    p1_1 = 0.5*( (*other).first + (*p1i).first );
		}
		
		// set upper limit of p1
		other = p1i; ++other;
		if(other == p1map.end())
		{
		    other = p1i; --other;  // get previous one
		    p1_2 = 1.5*(*p1i).first - 0.5*(*other).first;
		}
		else
		{
		    p1_2 = 0.5*( (*p1i).first + (*other).first );
		}
		
		if(p2i == p2map.begin())
		{
		    other = p2i; ++other; // get the next one
		    p2_1 = 1.5*(*p2i).first-0.5*(*other).first;
		}
		else
		{
		    other = p2i; --other;  // get the previous one
		    p2_1 = 0.5*( (*other).first + (*p2i).first );
		}
		
		// set upper limit of p2
		other = p2i; ++other;
		if(other == p2map.end())
		{
		    other = p2i; --other;  // get previous one
		    p2_2 = 1.5*(*p2i).first - 0.5*(*other).first;
		}
		else
		{
		    p2_2 = 0.5*( (*p2i).first + (*other).first );
		}
	    }

	    vector<double> x, y;

	    const int N = 5;

	    // fix the first param to p1_1, and scan the second param
	    param[0] = p1_1;  
	    for(int n=0; n<N; ++n)
	    {
		param[1] = p2_1 + n*(p2_2-p2_1)/N;
		params_to_xy_.meval(param,xy);
		x.push_back(xaxis->map_point(xy[0].dbl()));
		y.push_back(yaxis->map_point(xy[1].dbl()));
	    }

	    // now fix second param to p2_2, and scan the first param
	    param[1] = p2_2;
	    for(int n=0; n<N; ++n)
	    {
		param[0] = p1_1 + n*(p1_2-p1_1)/N;
		params_to_xy_.meval(param,xy);
		x.push_back(xaxis->map_point(xy[0].dbl()));
		y.push_back(yaxis->map_point(xy[1].dbl()));
	    }

	    // now fix the first param to p1_2, and scan the second param
	    param[0] = p1_2;
	    for(int n=0; n<N; ++n)
	    {
		param[1] = p2_2 - n*(p2_2-p2_1)/N;
		params_to_xy_.meval(param,xy);
		x.push_back(xaxis->map_point(xy[0].dbl()));
		y.push_back(yaxis->map_point(xy[1].dbl()));
	    }

	    // and finally fix second param to p2_1, and scan the first param
	    param[1] = p2_1;
	    for(int n=0; n<N; ++n)
	    {
		param[0] = p1_2 - n*(p1_2-p1_1)/N;
		params_to_xy_.meval(param,xy);
		x.push_back(xaxis->map_point(xy[0].dbl()));
		y.push_back(yaxis->map_point(xy[1].dbl()));
	    }

	    term->set_color(map_color(z,zmin,zmax));
	    draw_lines_raw(x,y,0,1,0,1,term,true);
	}
	//draw_legend(term,g);
    }


    // --------------------------------------------------------------------------------

    bool isolines::default_center_labels_ = false;

    isolines::isolines()
	: isovalues_specified_(false),
	  min_(unset), max_(unset), step_(unset),
	  minfixed_(false), maxfixed_(false), stepfixed_(false),
	  logscale_(false), draw_labels_(true), turn_labels_(true), labelformat_("%g"),
	  center_labels_(default_center_labels_),
	  x1repulsion_(1), x2repulsion_(1), y1repulsion_(1), y2repulsion_(1),
	  label_xalign_(sym::center), label_yalign_(sym::center)
    {
	labelrepulsion_[0] = 1;
	labelrepulsion_[1] = 0.7;
	labelrepulsion_[2] = 0.4;
	colors_.push_back(black);
	colors_.push_back(red);
	colors_.push_back(blue);
	colors_.push_back(green);
	colors_.push_back(magenta);
	colors_.push_back(cyan);
	colors_.push_back(orange);
    }

    isolines &isolines::label_align(sym::position xalign, sym::position yalign)
    {
	label_xalign_ = xalign;
	label_yalign_ = yalign;
	return *this;
    }

    isolines &isolines::label_repulsion_frame(double w)
    {
	x1repulsion_ = x2repulsion_ = y1repulsion_ = y2repulsion_ = w;
	return *this;
    }

    isolines &isolines::label_repulsion_label(double w1, double w2, double w3)
    {
	labelrepulsion_[0] = w1;
	labelrepulsion_[1] = w2;
	labelrepulsion_[2] = w3;
	return *this;
    }

    void isolines::set_ranges(plottable *g, axis *x, axis *y)
    {
	set_ranges_default(g,x,y,1,2);
    }

    bool isolines::skip_(const var &xorig_var,
			 const var &yorig_var,
			 const var &zorig_var,
			 axis *xaxis,
			 axis *yaxis,
			 double xmin, double xmax,
			 double ymin, double ymax )
    {
	if(ignore::it(xorig_var) ||
	   ignore::it(yorig_var) ||
	   ignore::it(zorig_var)) return true;
	
	const double xorig = xorig_var.dbl();
	const double yorig = yorig_var.dbl();
	
	if(xaxis->logscale() && xorig<=0.0) return true;
	if(yaxis->logscale() && yorig<=0.0) return true;

	if(xorig<xmin || xmax<xorig || yorig<ymin || ymax<yorig) return true;

	const double xx = xaxis->map_point(xorig);
	const double yy = yaxis->map_point(yorig);

	if(xx == unset || yy == unset) return true;
	if(xx<0 || 1<xx || yy<0 || 1<yy) return true;

	return false;
    }


    void isolines::draw(plottable *g, frame *f, terminal *term)
    {
#ifdef HAVE_GTS_H

	if(!g || g->empty()) return;
	axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());

	double xmin = xaxis->min();
	if(g->xmin()!=unset && g->xmin()>xmin) xmin = g->xmin();
	double xmax = xaxis->max();
	if(g->xmax()!=unset && g->xmax()<xmax) xmax = g->xmax();
	
	double ymin = yaxis->min();
	if(g->ymin()!=unset && g->ymin()>ymin) ymin = g->ymin();
	double ymax = yaxis->max();
	if(g->ymax()!=unset && g->ymax()<ymax) ymax = g->ymax();
	
	function getx = get_x(g);
	function gety = get_y(g);
	function getz = g->z_hint();
	if(!getz.initialized()) getz = _3;

	// first loop over all points to calculate number of vertices
	double zmin = unset, zmax = unset;
	if(minfixed_) zmin = min_;
	if(maxfixed_) zmax = max_;
	for(plottable::size_type i=0; i<g->size(); ++i)
	{
	    const var xorig_var = getx.eval(*(g->get(i)));
	    const var yorig_var = gety.eval(*(g->get(i)));
	    const var zorig_var = getz.eval(*(g->get(i)));
	    if(skip_(xorig_var,yorig_var,zorig_var,xaxis,yaxis,xmin,xmax,ymin,ymax))
		continue;
	    if(!minfixed_ && (zmin == unset || zorig_var.dbl()<zmin)) zmin=zorig_var.dbl();
	    if(!maxfixed_ && (zmax == unset || zorig_var.dbl()>zmax)) zmax=zorig_var.dbl();
	}

	if(!isovalues_specified_)
	{
	    isovalues_.clear();
	    std::vector<blop::tic> tics;
	    std::vector<std::pair<double,double> > cuts;
	    tic scale;
	    scale.value(1).label("1");
	    double zmin_local = zmin;
	    if(minfixed_) zmin_local = min_;
	    double zmax_local = zmax;
	    if(maxfixed_) zmax_local = max_;
	    double step = step_;

	    if(stepfixed_ && spec_val_ != unset)
	    {
		for(double v = spec_val_; v<zmax; v += step_)
		{
		    isovalues_.push_back(v);
		}
		for(double v = spec_val_-step_; v>zmin; v -= step_)
		{
		    isovalues_.push_back(v);
		}
	    }
	    else
	    {
		calculate_tics(zmin_local, minfixed_,  
			       zmax_local, maxfixed_,   
			       step, stepfixed_,  
			       unset, unset,
			       scale,
			       1.0,
			       cuts,
			       logscale_,
			       false, // normalform_tics
			       false, // normalform_scale
			       tics);
		
		for(unsigned int i=0; i<tics.size(); ++i)
		{
		    const double val = tics[i].value();
		    if(zmin<val && val<zmax) isovalues_.push_back(val);
		}
	    }
	}
	if(isovalues_.empty())
	{
	    warning::print("No isovalues could be calculated","isolines::draw(...)");
	    return;
	}

	delaunay_surface surface;

	// Now collect the vertices. The x and y coordinates already the
	// mapped values (that is, between 0 and 1)
	for(plottable::size_type aa=0; aa<g->size(); ++aa)
	{
	    int i=aa-1;
	    if(i<0) i = g->size()-1;
	    const var xorig_var = getx.eval(*(g->get(i)));
	    const var yorig_var = gety.eval(*(g->get(i)));
	    const var zorig_var = getz.eval(*(g->get(i)));
	    if(skip_(xorig_var,yorig_var,zorig_var,xaxis,yaxis,xmin,xmax,ymin,ymax)) continue;

	    const double x = xaxis->map_point(xorig_var.dbl());
	    const double y = yaxis->map_point(yorig_var.dbl());
	    const double z = zorig_var.dbl();

	    surface.add_vertex(x,y,z);
	}

	double last_label[3][2] = {{unset,unset},{unset,unset},{unset,unset}};

	term->reset_transformation();

	// now loop over the isovalues, create a cutting plane at this
	// value, and let gts calculate the intersection of the surface
	// with this cutting plane
	for(unsigned int izcut=0; izcut<isovalues_.size(); ++izcut)
	{
	    const double zcut = isovalues_[izcut];

	    if(isovalues_specified_ && (zcut<zmin || zmax<zcut))
	    {
		warning::print(var("Specified isovalue ") & zcut & var(" is out of range"));
		continue;
	    }

	    // take the intersection of the surface with a z=zcut plane surface,
	    // the (x,y) coordinates are returned in the xx and yy arrays, unset values
	    // indicating a discontinuity.
	    vector<double> xx, yy;
	    intersection_z(surface,zcut,xx,yy);

	    double largest_distance = 0;
	    double labelp1[2] = {unset,unset};
	    double labelp2[2] = {unset,unset};

	    for(unsigned int i=0; i<xx.size()-1; ++i)
	    {
		if(xx[i]==unset || xx[i+1]==unset) continue;

		// skip point-pairs which are too close
		const double dx = xx[i+1]-xx[i];
		const double dy = yy[i+1]-yy[i];
		if(std::sqrt(dx*dx+dy*dy)<1e-3) continue;

		// the coordinates of the isovalue label are set to the midpoint between this
		// and the next point
		const double label_x = 0.5*(xx[i]+xx[i+1]);
		const double label_y = 0.5*(yy[i]+yy[i+1]);

		double d =
		    x1repulsion_*::sqrt(::fabs(label_x)) +
		    x2repulsion_*::sqrt(::fabs(1-label_x)) +
		    y1repulsion_*::sqrt(::fabs(label_y)) +
		    y2repulsion_*::sqrt(::fabs(1-label_y));
		for(int a=0; a<3; ++a)
		{
		    if(last_label[a][0] != unset)
		    {
			d += labelrepulsion_[a] *
			    (::sqrt(::fabs(last_label[a][0]-label_x))+
			     ::sqrt(::fabs(last_label[a][1]-label_y)));
		    }
		}
		if(d >= largest_distance)
		{
		    labelp1[0] = xx[i];
		    labelp1[1] = yy[i];
		    labelp2[0] = xx[i+1];
		    labelp2[1] = yy[i+1];
		    largest_distance = d;
		}
	    }

	    term->set_linewidth(g->linewidth().termspecific_id());
	    term->set_linestyle(g->linestyle());
	    term->set_color(colors_[izcut%colors_.size()]);
	    draw_lines_raw(xx,yy,
			   xaxis->map_point(xmin), xaxis->map_point(xmax),
			   yaxis->map_point(ymin), yaxis->map_point(ymax),
			   term,false);

	    if(draw_labels_)
	    {
		if(center_labels_)
		{
		    for(unsigned int i=xx.size()/2; i<xx.size()-1; ++i)
		    {
			if(xx[i]!=unset && xx[i+1]!=unset)
			{
			    labelp1[0] = xx[i];
			    labelp1[1] = yy[i];
			    labelp2[0] = xx[i+1];
			    labelp2[1] = yy[i+1];
			    break;
			}
		    }
		}

		if(labelp1[0] != unset && labelp1[1] != unset &&
		   labelp2[0] != unset && labelp2[1] != unset)
		{
		    double dx = labelp2[0] - labelp1[0];
		    double dy = labelp2[1] - labelp1[1];

		    if(dx<0) { dx *= -1; dy *= -1; }

		    if(::fabs(dx) < numeric_limits<double>::epsilon()*100 && 
		       ::fabs(dy) < numeric_limits<double>::epsilon()*100)
		    {
			dx = 1;
			dy = 0;
		    }
		    else if(::fabs(dx)<numeric_limits<double>::epsilon()*100)
		    {
			dx = 0;
			dy = 1;
		    }
		    else if(::fabs(dy)<numeric_limits<double>::epsilon()*100)
		    {
			dy = 0;
			dx = 1;
		    }
		    else if(::fabs(dy/dx) < 0.001) dy = 0;
		    else if(::fabs(dx/dy) < 0.001) dx = 0;
		    else
		    {
			dy /= dx;
			dx = 1;
		    }

		    const double x  = 0.5*(labelp1[0]+labelp2[0]);
		    const double y  = 0.5*(labelp1[1]+labelp2[1]);

		    string lab;
		    if(!labels_.empty()) lab = labels_[izcut%labels_.size()].str();
		    else
		    {
			char s[256];
			sprintf(s,labelformat_.c_str(),zcut);
			lab = s;
		    }

		    term->translate(terminal::id(x ,1),terminal::id(y ,2));
		    if(turn_labels_)
		    {
			term->draw_text(terminal::coord(terminal::id(0.0,1), terminal::id(0.0,2)),
					lab, label_xalign_, label_yalign_, terminal::id(dx,1),terminal::id(dy,2));
		    }
		    else
		    {
			term->draw_text(terminal::coord(terminal::id(0.0,1), terminal::id(0.0,2)),
					lab, label_xalign_, label_yalign_, 0);
		    }
		    term->reset_transformation();

		    last_label[2][0] = last_label[1][0];
		    last_label[2][1] = last_label[1][1];
		    last_label[1][0] = last_label[0][0];
		    last_label[1][1] = last_label[0][1];
		    last_label[0][0] = x;
		    last_label[0][1] = y;
		}
		else
		{
		    warning::print(var("For some reason did not find the optimum "
				       "point to print the isovalue (") & zcut & ") label",
				   "isolines::draw(...)");
		}
	    }

	} // end of looping over isovalues (cutvalues)

// HAVE_GTS_H
#else
	warning::print("Install libgts-dev package to have isolines");
#endif
    }

    isolines &isolines::at(double val1,
			   double val2,
			   double val3,
			   double val4,
			   double val5,
			   double val6,
			   double val7,
			   double val8,
			   double val9,
			   double val10,
			   double val11,
			   double val12,
			   double val13,
			   double val14,
			   double val15,
			   double val16,
			   double val17,
			   double val18,
			   double val19,
			   double val20)
    {

	isovalues_.push_back(val1);
	if(val2!=unset) isovalues_.push_back(val2);
	if(val3!=unset) isovalues_.push_back(val3);
	if(val4!=unset) isovalues_.push_back(val4);
	if(val5!=unset) isovalues_.push_back(val5);
	if(val6!=unset) isovalues_.push_back(val6);
	if(val7!=unset) isovalues_.push_back(val7);
	if(val8!=unset) isovalues_.push_back(val8);
	if(val9!=unset) isovalues_.push_back(val9);
	if(val10!=unset) isovalues_.push_back(val10);
	if(val11!=unset) isovalues_.push_back(val11);
	if(val12!=unset) isovalues_.push_back(val12);
	if(val13!=unset) isovalues_.push_back(val13);
	if(val14!=unset) isovalues_.push_back(val14);
	if(val15!=unset) isovalues_.push_back(val15);
	if(val16!=unset) isovalues_.push_back(val16);
	if(val17!=unset) isovalues_.push_back(val17);
	if(val18!=unset) isovalues_.push_back(val18);
	if(val19!=unset) isovalues_.push_back(val19);
	if(val20!=unset) isovalues_.push_back(val20);
	isovalues_specified_ = true;
	return *this; 
    }


    isolines &isolines::colors(const color &c1,
			       const color &c2,
			       const color &c3,
			       const color &c4,
			       const color &c5,
			       const color &c6,
			       const color &c7,
			       const color &c8,
			       const color &c9,
			       const color &c10)
    {
	colors_.clear();
	colors_.push_back(c1);
	if(c2.red()>=0) colors_.push_back(c2);
	if(c3.red()>=0) colors_.push_back(c3);
	if(c4.red()>=0) colors_.push_back(c4);
	if(c5.red()>=0) colors_.push_back(c5);
	if(c6.red()>=0) colors_.push_back(c6);
	if(c7.red()>=0) colors_.push_back(c7);
	if(c8.red()>=0) colors_.push_back(c8);
	if(c9.red()>=0) colors_.push_back(c9);
	if(c10.red()>=0) colors_.push_back(c10);
	return *this;
    }

    isolines &isolines::labels(const var &c1,
			       const var &c2,
			       const var &c3,
			       const var &c4,
			       const var &c5,
			       const var &c6,
			       const var &c7,
			       const var &c8,
			       const var &c9,
			       const var &c10)
    {
	labels_.clear();
	labels_.push_back(c1);
	if(c2.dbl() != unset) labels_.push_back(c2);
	if(c3.dbl() != unset) labels_.push_back(c3);
	if(c4.dbl() != unset) labels_.push_back(c4);
	if(c5.dbl() != unset) labels_.push_back(c5);
	if(c6.dbl() != unset) labels_.push_back(c6);
	if(c7.dbl() != unset) labels_.push_back(c7);
	if(c8.dbl() != unset) labels_.push_back(c8);
	if(c9.dbl() != unset) labels_.push_back(c9);
	if(c10.dbl() != unset) labels_.push_back(c10);
	return *this;
    }

    graph_drawer *isolines::clone() const
    {
	return new isolines(*this);
    }

    void isolines::prepare_for_draw(plottable *g, frame *, int count)
    {
	if(count==1) g->linewidth().register_me();
    }

    void isolines::draw_sample(const length &x,
			       const length &y,
			       const length &size,
			       const plottable *style,
			       terminal *t)
    {
	t->set_linewidth(style->linewidth().termspecific_id());
	t->set_linestyle(style->linestyle());
	length l1 = !x - 0.5*!size;
	length l2 = !x + 0.5*!size;
	length l3 = height(style->legend()); //3*style->linewidth();
	l1.specialize(t);
	l2.specialize(t);
	l3.specialize(t);
	for ( unsigned int l=0 ; l<isovalues_.size() ; ++l )
	{
	    double n = isovalues_.size();
	    length l5 = !y + (l-n/2)/n*l3;
	    l5.specialize(t);
	    t->set_color(colors_[l%colors_.size()]);
	    t->draw_line(terminal::coord(l1.termspecific_id(), l5.termspecific_id()),
			 terminal::coord(l2.termspecific_id(), l5.termspecific_id()));
	}
    }

    // ---------- vectors  ------------------------------------------------------------

    vectors::~vectors()
    {
	delete arrow_;
    }

    void vectors::setup_when_added(plottable *g,frame *f)
    {
	if(g && f && use_color_) graphd_colorscale::setup_when_added(g,f);
	else colorlegend_ = 0;
    }

    vectors::vectors(const vectors &rhs) : graphd_colorscale(rhs)
    {
	colorfunc_ = rhs.colorfunc_;
	colorlegend_ = 0;
	scale_arrow_ = rhs.scale_arrow_;
	use_color_ = rhs.use_color_;
	xunit_ = rhs.xunit_;
	yunit_ = rhs.yunit_;
	arrow_ = rhs.arrow_->clone();
	dx_ = rhs.dx_;
	dy_ = rhs.dy_;
	norm_ = rhs.norm_;
	norm_fixed_ = rhs.norm_fixed_;
	min_length_cut_ = rhs.min_length_cut_;
	min_length_cut_fixed_ = rhs.min_length_cut_fixed_;
	pos_ = rhs.pos_;
	clip_ = rhs.clip_;
    }

    vectors::vectors()
    {
	arrow_ = new arrowhead::simple(EX);
	dx_ = dy_ = norm_ = unset;
	min_length_cut_ = unset;
	norm_fixed_ = false;
	min_length_cut_fixed_ = false;
	pos_ = sym::center;
	min_length_ = max_length_ = unset;
	scale_arrow_ = true;
	use_color_ = false;
	colorfunc_ = sqrt(_3*_3+_4*_4);
	colorlegend_ = 0;
	clip_ = true;
    }

    graph_drawer *vectors::clone() const
    {
	return new vectors(*this);
    }

    int vectors::req_components() const
    {
	int result = 4;
	if(use_color_ && colorfunc_.nargs()>4) result = colorfunc_.nargs();
	return result;
    }

    void vectors::prepare_for_draw(plottable *g, frame *f, int count)
    {
	if(g->empty()) return;

	if(count==1)
	{
	    g->linewidth().register_me();
	    
	    axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	    axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());
	    
	    function getx=g->x_hint();
	    if (!getx.initialized()) getx=ARG(1);
	    function gety=g->y_hint();
	    if (!gety.initialized()) gety=ARG(2);
	    function getvx = ARG(3);
	    function getvy = ARG(4);
	    
	    map<double,bool> xvalues,yvalues;

	    double xmin = unset, ymin = unset, xmax = unset, ymax = unset;
	    if(xaxis->min_fixed()) xmin = xaxis->min();
	    if(xaxis->max_fixed()) xmax = xaxis->max();
	    if(yaxis->min_fixed()) ymin = yaxis->min();
	    if(yaxis->max_fixed()) ymax = yaxis->max();
	    
	    if(g->xmin() != unset)
	    {
		if(xmin == unset || g->xmin() > xmin) xmin = g->xmin();
	    }
	    if(g->xmax() != unset)
	    {
		if(xmax == unset || g->xmax() < xmax) xmax = g->xmax();
	    }
	    if(g->ymin() != unset)
	    {
		if(ymin == unset || g->ymin() > ymin) ymin = g->ymin();
	    }
	    if(g->ymax() != unset)
	    {
		if(ymax == unset || g->ymax() < ymax) ymax = g->ymax();
	    }
	    
	    if(!norm_fixed_) norm_ = 0;

	    max_length_ = min_length_ = unset;
	    if(!color_min_fixed_) color_min_ = unset;
	    if(!color_max_fixed_) color_max_ = unset;
	    for(plottable::size_type i=0; i<g->size(); ++i)
	    {
		const var x = getx.eval(*(g->get(i)));
		const var y = gety.eval(*(g->get(i)));
		const var vx = getvx.eval(*(g->get(i)));
		const var vy = getvy.eval(*(g->get(i)));
		
		if(ignore::it(x) || ignore::it(y) || ignore::it(vx) || ignore::it(vy)) continue;
		if(xmin != unset && x.dbl() < xmin) continue;
		if(xmax != unset && x.dbl() > xmax) continue;
		if(ymin != unset && y.dbl() < ymin) continue;
		if(ymax != unset && y.dbl() > ymax) continue;
		
		xvalues[x.dbl()] = true;
		yvalues[y.dbl()] = true;
		
		const double l = ::sqrt(vx.dbl()*vx.dbl()+vy.dbl()*vy.dbl());
		if(max_length_ == unset || l > max_length_) max_length_ = l;
		if(min_length_ == unset || l < min_length_) min_length_ = l;
		
		if(use_color_)
		{
		    double c = colorfunc_.eval(*(g->get(i)));
		    if(!color_min_fixed_) if(color_min_ == unset || c < color_min_) color_min_ = c;
		    if(!color_max_fixed_) if(color_max_ == unset || c > color_max_) color_max_ = c;
		}
	    }

	    if(!min_length_cut_fixed_) min_length_cut_ = 0.05*max_length_;
	    min_length_ = min_length_cut_;
	    if(!norm_fixed_) norm_   = max_length_;
	    
	    dx_ = unset;
	    dy_ = unset;
	    double xaver = 0, yaver = 0;
	    for(map<double,bool>::iterator i=xvalues.begin(); i!=xvalues.end(); ++i)
	    {
		xaver += (*i).first;
		if(i != xvalues.begin())
		{
		    map<double,bool>::iterator prev = i;
		    --prev;
		    if(dx_ == unset || ::fabs((*prev).first - (*i).first)<dx_)
		    {
			dx_ = ::fabs((*prev).first - (*i).first);
		    }
		}
	    }
	    for(map<double,bool>::iterator i=yvalues.begin(); i!=yvalues.end(); ++i)
	    {
		yaver += (*i).first;
		if(i != yvalues.begin())
		{
		    map<double,bool>::iterator prev = i;
		    --prev;
		    if(dy_ == unset || ::fabs((*prev).first - (*i).first)<dy_)
		    {
			dy_ = ::fabs((*prev).first - (*i).first);
		    }
		}
	    }
	    xaver /= xvalues.size();
	    yaver /= yvalues.size();
	    
	    if(dx_ == unset) dx_ = 1;
	    if(dy_ == unset) dy_ = 1;
	    double delta = ::min(dx_,dy_);
	    
	    xunit_ = length::base_axis_t(xaxis, xaver, xaver+delta);
	    yunit_ = length::base_axis_t(yaxis, yaver, yaver+delta);
	    xunit_.register_me();
	    yunit_.register_me();
	    
	    arrow_->prepare_for_draw();

	    if(colorlegend_) colorlegend_->clear_tics();
	} // end if(count==1)

	if(count==2)
	{
	    // Let the colorlegend calculate the unified range for all
	    // of its graphs
	    if(colorlegend_) colorlegend_->calculate_tics();
	}
    }

    vectors &vectors::pos(sym::position p)
    {
	if(p != sym::center && p != sym::begin && p != sym::end)
	{
	    warning::print("Argument should be either of sym::begin, sym::end or sym::center. Using sym::center!",
			   "vectors::pos(symbolic p)");
	    p = sym::center;
	    
	}
	pos_ = p;
	return *this; 
    }


    vectors &vectors::arrow(const arrowhead &a)
    {
	delete arrow_;
	arrow_ = a.clone();
	return *this; 
    }

    vectors &vectors::norm(double value)
    {
	norm_ = value;
	if(value == unset) norm_fixed_ = false;
	else norm_fixed_ = true;
	return *this;
    }

    vectors &vectors::min(double value)
    {
	min_length_cut_ = value;
	if(value == unset) min_length_cut_fixed_ = false;
	else min_length_cut_fixed_ = true;
	return *this;
    }


    void vectors::set_ranges(plottable *g, axis *xaxis, axis *yaxis)
    {
	if(!g) err("Plottagle to draw is zero");
	if(!xaxis) err("Xaxis not set");
	if(!yaxis) err("yaxis not set");

	if(g->empty())
	{                                                                       
	    warning::print("Warning: empty plottable!",
			   "void set_ranges_default(plottable*, axis*, axis*, int, int)");
	    return;                                                             
	}

	if(g->xmin()!=unset) xaxis->extend_range(g->xmin());
	if(g->xmax()!=unset) xaxis->extend_range(g->xmax());
	if(g->ymin()!=unset) yaxis->extend_range(g->ymin());
	if(g->ymax()!=unset) yaxis->extend_range(g->ymax());

	// if all ranges are defined, then simply return at this point
	if(g->xmin()!=unset && g->xmax()!=unset &&
	   g->ymin()!=unset && g->ymax()!=unset) return;

	double xmin=unset;
	double xmax=unset;
	double ymin=unset;
	double ymax=unset;
	function getx=g->x_hint();
	if (!getx.initialized()) getx=_1;
	function gety=g->y_hint();
	if (!gety.initialized()) gety=_2;

	for ( unsigned int i=0 ; i<g->size() ; ++i )
	{
	    const var xorig_var=getx.eval(*g->get(i));
	    if( ignore::it(xorig_var) ) continue;
	    const double xorig=xorig_var.dbl();
	    if( xaxis->logscale() && xorig<=0.0 ) continue;
	    if( g->xmin() != unset && xorig<g->xmin() ) continue;
	    if( g->xmax() != unset && xorig>g->xmax() ) continue;
	    if(xmin == unset || xorig<xmin) xmin = xorig;
	    if(xmax == unset || xorig>xmax) xmax = xorig;

	    const var yorig_var=gety.eval(*g->get(i));
	    if ( ignore::it(yorig_var) ) continue;
	    const double yorig=yorig_var.dbl();
	    if( yaxis->logscale() && yorig<=0.0 ) continue;
	    if( g->ymin()!=unset && yorig<g->ymin()) continue;
	    if( g->ymax()!=unset && yorig>g->ymax()) continue;
	    if(ymin == unset || yorig < ymin) ymin = yorig;
	    if(ymax == unset || yorig > ymax) ymax = yorig;

	}

	if(xmin != unset) xmin -= dx_/2;
	if(xmax != unset) xmax += dx_/2;
	if(ymin != unset) ymin -= dy_/2;
	if(ymax != unset) ymax += dy_/2;
	if (g->xmin()==unset && xmin!=unset) xaxis->extend_range(xmin);
	if (g->xmax()==unset && xmax!=unset) xaxis->extend_range(xmax);
	if (g->ymin()==unset && ymin!=unset) yaxis->extend_range(ymin);
	if (g->ymax()==unset && ymax!=unset) yaxis->extend_range(ymax);	
    }

    void vectors::draw(plottable *g, frame *f, terminal *t)
    {
	if(g->empty()) return;

	bool call_noclip = false;
	if(clip_) call_noclip = t->clip(terminal::coord(terminal::id(0,1),terminal::id(0,2)),
					terminal::coord(terminal::id(1,1),terminal::id(1,2)));

	axis *xaxis = (g->xaxis() == axis::x1 ? f->x1axis() : f->x2axis());
	axis *yaxis = (g->yaxis() == axis::y1 ? f->y1axis() : f->y2axis());

	function getx=g->x_hint();
	if (!getx.initialized()) getx=ARG(1);
	function gety=g->y_hint();
	if (!gety.initialized()) gety=ARG(2);
	function getvx = ARG(3);
	function getvy = ARG(4);

	int x1id = t->newlength();
	int y1id = t->newlength();
	int x2id = t->newlength();
	int y2id = t->newlength();

	double xmin = unset, ymin = unset, xmax = unset, ymax = unset;
	if(xaxis->min_fixed()) xmin = xaxis->min();
	if(xaxis->max_fixed()) xmax = xaxis->max();
	if(yaxis->min_fixed()) ymin = yaxis->min();
	if(yaxis->max_fixed()) ymax = yaxis->max();

	if(g->xmin() != unset)
	{
	    if(xmin == unset || g->xmin() > xmin) xmin = g->xmin();
	}
	if(g->xmax() != unset)
	{
	    if(xmax == unset || g->xmax() < xmax) xmax = g->xmax();
	}
	if(g->ymin() != unset)
	{
	    if(ymin == unset || g->ymin() > ymin) ymin = g->ymin();
	}
	if(g->ymax() != unset)
	{
	    if(ymax == unset || g->ymax() < ymax) ymax = g->ymax();
	}

	t->set_color(g->linecolor());
	t->set_linestyle(g->linestyle());
	t->set_linewidth(g->linewidth().termspecific_id());

        int plotted = 0;

	for(plottable::size_type i = 0; i<g->size(); ++i)
	{
	    var x = getx.eval(*(g->get(i)));
	    var y = gety.eval(*(g->get(i)));
	    var vx = getvx.eval(*(g->get(i)));
	    var vy = getvy.eval(*(g->get(i)));

	    if(ignore::it(x) || ignore::it(y) || ignore::it(vx) || ignore::it(vy)) continue;
	    if(xmin != unset && x.dbl() < xmin) continue;
	    if(xmax != unset && x.dbl() > xmax) continue;
	    if(ymin != unset && y.dbl() < ymin) continue;
	    if(ymax != unset && y.dbl() > ymax) continue;

	    const double length = ::sqrt(vx.dbl()*vx.dbl() + vy.dbl()*vy.dbl());
	    if(min_length_cut_ != unset && length<min_length_cut_) continue;

	    double xx = xaxis->map_point(x.dbl());
	    double yy = yaxis->map_point(y.dbl());

	    const double vxdbl = vx.dbl()/norm_;
	    const double vydbl = vy.dbl()/norm_;

            if(!finite(vxdbl) || !finite(vydbl)) continue;

	    if(use_color_)
	    {
		double minimum = color_min();
		double maximum = color_max();
		if(colorlegend_)
		{
		    minimum = colorlegend_->min();
		    maximum = colorlegend_->max();
		}
		const double colorkey = colorfunc_.eval(*(g->get(i)));
		t->set_color(cmapping_.map(colorkey, minimum, maximum));
	    }
            else
            {
                t->set_color(g->linecolor());
            }

	    if(pos_ == sym::center)
	    {
		t->overwrite(x1id, 1, terminal::id(xx,1), -0.5*vxdbl, xunit_.termspecific_id());
		t->overwrite(x2id, 1, terminal::id(xx,1),  0.5*vxdbl, xunit_.termspecific_id());
		t->overwrite(y1id, 1, terminal::id(yy,2), -0.5*vydbl, yunit_.termspecific_id());
		t->overwrite(y2id, 1, terminal::id(yy,2),  0.5*vydbl, yunit_.termspecific_id());
	    }
	    else if(pos_ == sym::begin)
	    {
		t->overwrite(x1id, 1, terminal::id(xx,1),  0*vxdbl, xunit_.termspecific_id());
		t->overwrite(x2id, 1, terminal::id(xx,1),  1*vxdbl, xunit_.termspecific_id());
		t->overwrite(y1id, 1, terminal::id(yy,2),  0*vydbl, yunit_.termspecific_id());
		t->overwrite(y2id, 1, terminal::id(yy,2),  1*vydbl, yunit_.termspecific_id());
	    }
	    else
	    {
		t->overwrite(x1id, 1, terminal::id(xx,1), -1*vxdbl, xunit_.termspecific_id());
		t->overwrite(x2id, 1, terminal::id(xx,1),  0*vxdbl, xunit_.termspecific_id());
		t->overwrite(y1id, 1, terminal::id(yy,2), -1*vydbl, yunit_.termspecific_id());
		t->overwrite(y2id, 1, terminal::id(yy,2),  0*vydbl, yunit_.termspecific_id());
	    }

	    t->draw_line(terminal::coord(terminal::id(x1id),terminal::id(y1id)),
			 terminal::coord(terminal::id(x2id),terminal::id(y2id)));

	    t->overwrite(x1id, -1, x1id, 1, x2id);
	    t->overwrite(y1id, -1, y1id, 1, y2id);

	    if(scale_arrow_) t->scale(length/norm_);
	    t->rotate(x1id,y1id);
	    t->translate(terminal::id(x2id),terminal::id(y2id));
	    arrow_->print(t);
	    t->reset_transformation();
            ++plotted;
	}
	if(call_noclip) t->noclip();

        if(plotted == 0) warning::print("No vector could be plotted (graph empty or all vectors have 0 length)","vectors::draw(...)");
    }



    void vectors::draw_sample(const length &x, const length &y, const length &size,
			      const plottable *g, terminal *t)
    {
	t->set_linewidth(g->linewidth().termspecific_id());
	t->set_color(g->linecolor());
	t->set_linestyle(g->linestyle());
	length l1 = !x - 0.5*!size;
	length l2 = !x + 0.5*!size;
	l1.specialize(t);
	l2.specialize(t);
	t->draw_line(terminal::coord(l1.termspecific_id(), y.termspecific_id()),
		     terminal::coord(l2.termspecific_id(), y.termspecific_id()));

	t->translate(l2.termspecific_id(),y.termspecific_id());
	arrow_->print(t);
	t->reset_transformation();
    }



}